# Fluctuations in the canonical ensemble of an Abelian charge

Bengt Friman  GSI Helmholtzzentrum für Schwerionenforschung, 64291 Darmstadt,
Germany  Krzysztof Redlich  Institute of Theoretical Physics, University of
Wrocław, 50204 Wrocław, Poland

###### Abstract

We study fluctuations in the canonical ensemble, where the net baryon number
is exactly conserved. The focus is on cumulants and factorial cumulants linked
to the baryon and antibaryon multiplicities and their sum or difference in
full phase-space as well as in subsystems. In particular, we connect the
fluctuations of the net baryon number in a subsystem, relevant for fluctuation
studies in nucleus-nucleus collisions, with fluctuations of the baryon and
antibaryon numbers of the total system. We derive analytic expressions for
factorial cumulants of arbitrary order. Compact results are obtained in terms
of cumulants of the baryon number of the total system. Moreover, we derive the
asymptotic forms of the factorial cumulants of baryon and antibaryon
multiplicities in the high- and low-temperature limits and discuss the results
in the context of heavy-ion collision experiments.

##  I Introduction

One of the goals of current experimental and theoretical studies of chiral
symmetry restoration in QCD is to unravel the phase diagram of strongly
interacting matter and to verify whether a chiral critical endpoint exists. A
dedicated research program at RHIC, the beam energy scan (BES), has been
established to explore these issues in collisions of heavy ions at
relativistic energies [1]. Fluctuations and correlations of conserved charges
are considered as possible probes of the QCD phase diagram [2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14]. These are experimentally accessible observables and
reflect the criticality of the chiral transition. Fluctuations of the net
baryon number are particularly interesting, owing to a direct connection to
critical scaling near the chiral phase boundary [5, 10]. Other useful
quantities characterizing multiparticle-correlations, and possibly also the
chiral phase transition, are factorial cumulants of the baryon and antibaryon
multiplicity distributions [13, 15, 16, 17, 18]. Data on proton-, antiproton-
and net-proton-number fluctuations in heavy-ion collisions have been obtained
by the STAR Collaboration at RHIC [19, 20, 21, 22] and the ALICE Collaboration
[23, 24, 25] at LHC energies. These are utilized as proxies for fluctuations
of the baryon, antibaryon and net baryon numbers, respectively. The
corresponding results on the variance, skewness, and kurtosis of the net
proton number are intriguing and have triggered discussions on their physics
origin and interpretation, focused mainly on the connection to chiral
criticality.

However, an important aspect of the analysis of such fluctuations in heavy-ion
collisions, is that the resulting cumulants are affected also by fluctuations
unrelated to criticality. Two such effects are of particular relevance: volume
fluctuations [27, 28, 29], which are linked to event-by-event fluctuations of
the number of participating nucleons, and constraints imposed on fluctuation
observables by exact conservation of the net baryon number in full phase space
[13, 15, 16, 17, 18, 28, 31, 30, 32, 33, 34, 35, 36, 37]

Recently, the constraints on fluctuations set by the global conservation of
additive charges, in particular of the net baryon number, were explored for
current experimental setups. In these studies cumulants of net-baryon-number
fluctuations [13, 28, 31, 30, 32, 34, 35, 36], as well as factorial cumulants
of the baryon and antibaryon multiplicity distributions [13, 15, 18] were
considered. It was shown that the magnitude and energy dependence of the
suppression of net-proton number cumulants relative to the Skellam baseline,
observed by the STAR and ALICE collaborations in ultra-relativistic nucleus-
nucleus collisions, is, given the present experimental uncertainties,
consistent with exact conservation of the total net baryon number [32, 34].

In the following, we extend the so far published results on baryon
conservation effects on fluctuation observables. In particular, we connect the
fluctuations of the net baryon number in a subsystem, relevant for fluctuation
studies in nucleus-nucleus collisions, with fluctuations of the baryon and
antibaryon numbers in the total system. Furthermore, taking into account exact
conservation of the baryon number, we compute factorial cumulants of the
baryon, antibaryon, mixed baryon-antibaryon and net baryon number multiplicity
distributions in a subsystem. We derive general expressions for these
factorial cumulants up to arbitrary order. Compact results are obtained in
terms of the cumulants of the baryon number of the total system. These results
are an extension of the previous finding by Barej and Bzdak [15, 18], where
analytic results for the corresponding factorial cumulants up to sixth order
were derived.

Considering possible applications of our results in heavy-ion collisions we
also derive asymptotic forms for the factorial cumulants of baryon and
antibaryon multiplicities in two limits, which correspond to the conditions in
heavy-ion collision experiments at high and low beam energies, respectively.
The analytic expressions for cumulants and factorial cumulants of arbitrary
order, as well as their relations, allow a systematic approach to fluctuations
in the canonical ensemble, which in turn provide an important baseline for the
analysis of fluctuation data in heavy-ion collisions.

The paper is organized as follows: In the next Section we introduce the
canonical partition function for uncorrelated baryons, where the net baryon
number is conserved and discuss pertinent fluctuation observables. In Section
3 we formulate the cumulant and factorial cumulant generating functions.
Analytic expressions for factorial cumulants of the baryon, antibaryon, and
the net baryon numbers in a subsystem are obtained in Section 4, while the
high- and low-energy limits of the cumulants and factorial cumulants are
presented in Section 5. We summarize our results in Section 6. Mathematical
details and general results, which are independent of the assumed partition
function, are presented in three appendices.

##  II Canonical ensemble of the net baryon number

To compute the influence of exact charge conservation on fluctuation
observables, we adopt a thermal model for the net baryon number conservation
in the Boltzmann approximation and neglect baryon-baryon interactions. The
statistical operator is formulated following the S-matrix approach [38, 39,
40, 41, 42, 43, 44, 45], where to leading order in the fugacity expansion it
has the form of an ideal gas, albeit with the thermal phase-space of free
baryons modified by meson-baryon interactions.

The S-matrix thermodynamic potential reproduces particle production yields in
heavy-ion collisions [46, 47, 48, 49], and describes basic properties of the
net proton number cumulants and their energy dependence, as obtained by the
STAR Collaboration [34]. Moreover, the S-matrix thermodynamic potential of a
hadron gas formulated in the grand canonical ensemble is consistent with the
lattice QCD equation of state and some second-order cumulants and correlations
of conserved charges in the confined phase [50, 51, 52, 53, 54].

In the grand canonical partition function

$\displaystyle\mathcal{Z}(\mu_{B},T)$ $\displaystyle=$
$\displaystyle\sum_{N_{B}=0}^{\infty}\sum_{N_{\bar{B}}=0}^{\infty}\frac{(e^{\mu_{B}/T}\,z_{B})^{N_{B}}}{N_{B}!}\frac{(e^{-\mu_{B}/T}\,z_{\bar{B}})^{N_{\bar{B}}}}{N_{\bar{B}}!}$
(1) $\displaystyle=$
$\displaystyle\exp\big{(}e^{\mu_{B}/T}\,z_{B}+e^{-\mu_{B}/T}\,z_{\bar{B}}\big{)}.$

the fluctuations of baryons and antibaryons are described by two independent
Poisson distributions. In the spirit of the S-matrix approach, the effect of
meson-baryon interactions is subsumed in the baryon and antibaryon single-
particle partition functions $z_{B}$ and $z_{\bar{B}}$.

The canonical partition function in the ensemble, where the net baryon number
is conserved, is given by [34, 46] 111As noted in [34], this form of the
partition function applies also to systems where the composition is non-
uniform, provided the baryon and antibaryon multiplicities are locally Poisson
distributed. Consequently, it is more general than the global thermal
statistical model employed as a motivation.

$\displaystyle Z_{B}$ $\displaystyle=$
$\displaystyle\sum_{N_{B}=0}^{\infty}\sum_{N_{\bar{B}}=0}^{\infty}\frac{(\lambda_{B}\,z_{B})^{N_{B}}}{N_{B}!}\frac{(\lambda_{\bar{B}}\,z_{\bar{B}})^{N_{\bar{B}}}}{N_{\bar{B}}!}\delta({N_{B}-N_{\bar{B}}-B})$
$\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{d\phi}{2\pi}\,e^{-iB\phi}\exp\left({\lambda_{B}\,z_{B}\,e^{i\phi}+\lambda_{\bar{B}}\,z_{\bar{B}}\,e^{-i\phi}}\right)$
$\displaystyle=$
$\displaystyle\left(\frac{\lambda_{B}\,z_{B}}{\lambda_{\bar{B}}\,z_{\bar{B}}}\right)^{\frac{B}{2}}\,I_{B}(2\,z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}}),$

where the auxiliary parameters $\lambda_{B,\bar{B}}$ are introduced for the
calculation of the mean number of baryons and antibaryons and the
corresponding cumulants,

$\displaystyle\langle N_{B}\rangle$ $\displaystyle=$
$\displaystyle\lambda_{B}{\frac{\partial\ln
Z_{B}}{\partial\lambda_{B}}}|_{\lambda_{B},\lambda_{\bar{B}}=1}=z\,\frac{I_{B-1}(2\,z)}{I_{B}(2\,z)},$
(3) $\displaystyle\langle N_{\bar{B}}\rangle$ $\displaystyle=$
$\displaystyle\lambda_{\bar{B}}{\frac{\partial\ln
Z_{B}}{\partial\lambda_{\bar{B}}}}|_{\lambda_{B},\lambda_{\bar{B}}=1}=z\,\frac{I_{B+1}(2\,z)}{I_{B}(2\,z)},$
(4) $\displaystyle\langle(\delta N_{\rm B})^{k}\rangle_{c}$ $\displaystyle=$
$\displaystyle\left[\left(\lambda_{B}\,\frac{\partial}{\partial\lambda_{B}}\right)^{k}\log
Z_{B}\right]_{\lambda_{B},\lambda_{\bar{B}}=1}\qquad(k>1).$ (5)

Here (and below) $\langle\dots\rangle$ denotes an expectation value in the
canonical ensemble, while $\langle\dots\rangle_{c}$ is the connected part
thereof. Thus, fluctuations of the baryon and antibaryon numbers in the
canonical ensemble are given by $\delta N_{B}=N_{B}-\langle N_{B}\rangle$ and
$\delta N_{\bar{B}}=N_{\bar{B}}-\langle N_{\bar{B}}\rangle$, and (5) yields
the cumulants of the baryon number in the full system.

We can rewrite (5) in terms of the mean baryon and antibaryon mulitpicities,
keeping the dependence on the auxiliary parameters $\lambda_{\rm
B},\lambda_{\bar{\rm B}}$

$\langle(\delta N_{\rm
B})^{k}\rangle_{c}=\frac{1}{2}\,\left[\left(\lambda_{B}\,\frac{\partial}{\partial\lambda_{B}}\right)^{k-1}\left(\langle
N_{B}\rangle_{\lambda}+\langle
N_{\bar{B}}\rangle_{\lambda}\right)\right]_{\lambda_{B},\lambda_{\bar{B}}=1}\quad(k>1),$
(6)

where

$\displaystyle\langle N_{B}\rangle_{\lambda}$ $\displaystyle=$ $\displaystyle
z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}}\,\frac{I_{B-1}(2\,z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}})}{I_{B}(2\,z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}})},$
(7) $\displaystyle\langle N_{\bar{B}}\rangle_{\lambda}$ $\displaystyle=$
$\displaystyle
z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}}\,\frac{I_{B+1}(2\,z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}})}{I_{B}(2\,z\,\sqrt{\lambda_{B}\,\lambda_{\bar{B}}})}.$
(8)

Hence, the baryon and antibaryon multiplicities, (7) and (8), are functions of
$z\sqrt{\lambda_{B}\,\lambda_{\bar{B}}}$. This reflects the fact that in the
canonical ensemble $N_{B}-N_{\bar{B}}$ is conserved and thus does not
fluctuate 222This is easily seen by noting that the recurrence relation
$2\,\nu\,I_{\nu}(x)=x\,\big{[}I_{\nu-1}(x)-I_{\nu+1}(x)\big{]}$ implies that
$\langle N_{B}\rangle_{\lambda}-\langle N_{\bar{B}}\rangle_{\lambda}=B$ is
independent of $\lambda_{B}$ and $\lambda_{\bar{B}}$.. It follows that the
fluctuations of $N_{B}$ and $N_{\bar{B}}$ are equal, i.e., $\delta N_{\rm
B}=\delta N_{\bar{\rm B}}$ and that cumulants of the form

$\langle(\delta N_{\rm B})^{n-m}\,(\delta N_{\bar{\rm B}})^{m}\rangle_{c}$ (9)

are independent of $m$, and thus all equal to $\langle(\delta N_{\rm
B})^{n}\rangle_{c}$. In the full system, it is therefore sufficient to
consider only cumulants of the baryon number.

Moreover, it follows that the derivatives of $\langle
N_{B}\rangle_{\lambda}+\langle N_{\bar{B}}\rangle_{\lambda}$ with respect to
$\lambda_{B}$ can be replaced by derivatives with respect to $z$, setting
$\lambda_{B}=\lambda_{\bar{B}}=1$. Using the notation

$\displaystyle c_{1}=$ $\displaystyle\frac{1}{2}\langle N_{\rm B}+N_{\bar{\rm
B}}\rangle,$ (10) $\displaystyle c_{k}=$ $\displaystyle\langle(\delta N_{\rm
B})^{k}\rangle_{c}\qquad(k>1),$ (11)

we can then rewrite (6) as

$c_{k}=\left(\frac{z}{2}\,\frac{d}{dz}\right)^{k-1}c_{1}\qquad(k>1).$ (12)

Clearly, the cumulants satisfy the recurrence relation

$c_{k+1}=\frac{z}{2}\,\frac{d}{dz}\,c_{k}.$ (13)

It is useful to define $N^{(tot)}=N_{B}+N_{\bar{B}}$ and to consider the
fluctuations thereof, $\delta N^{(tot)}=\delta N_{B}+\delta N_{\bar{B}}$. We
introduce the notation

$\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle\langle
N^{(tot)}\rangle,$ (14) $\displaystyle C_{k}$ $\displaystyle=$
$\displaystyle\langle(\delta N^{(tot)})^{k}\rangle_{c}\qquad(k>1).$

Using the fact that $\delta N_{\rm B}=\delta N_{\bar{\rm B}}$ and (10,11,13),
it follows that

$C_{k}=2^{k}\,c_{k}$ (15)

and

$C_{k+1}=z\,\frac{d}{dz}\,C_{k}.$ (16)

It is straighforward to compute the baryon number cumulants $c_{k}$ using the
recurrence relation for $c_{k}$, (13). Given that
$c_{1}=\frac{1}{2}(N_{B}+N_{\bar{B}})$, and using the notation [34]

$\displaystyle S$ $\displaystyle=$ $\displaystyle\langle
N_{B}+N_{\bar{B}}\rangle,$ $\displaystyle P$ $\displaystyle=$
$\displaystyle\langle N_{B}\rangle\,\langle N_{\bar{B}}\rangle,$ (17)
$\displaystyle Q$ $\displaystyle=$ $\displaystyle z^{2}-P,$ $\displaystyle W$
$\displaystyle=$ $\displaystyle Q\,S-P.$

with the derivatives

$S^{\prime}=\frac{4}{z}\,Q,\qquad P^{\prime}=\frac{2}{z}\,Q\,S,\qquad
Q^{\prime}=\frac{2}{z}\,\big{(}Q-W\big{)},$ (18)

one finds for the first few

$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle S/2,$ $\displaystyle
c_{2}$ $\displaystyle=$ $\displaystyle Q,$ (19) $\displaystyle c_{3}$
$\displaystyle=$ $\displaystyle Q-W,$ $\displaystyle c_{4}$ $\displaystyle=$
$\displaystyle Q-W\,+W\,S-2Q^{2}.$

The corresponding expressions for the cumulants $C_{k}$ are trivially obtained
using (15).

##  III Cumulant generating functions

The generating function for the cumulants of $N^{(tot)}=N_{B}+N_{\bar{B}}$,
(14), in the canonical ensemble with net baryon number $B$ is obtained by
evaluating

$\displaystyle G_{c}(t)$ $\displaystyle=$
$\displaystyle\ln\Big{[}\sum_{N_{B}=0}^{\infty}\sum_{N_{\bar{B}}=0}^{\infty}\frac{(z_{B})^{N_{B}}}{N_{B}!}\frac{(z_{\bar{B}})^{N_{\bar{B}}}}{N_{\bar{B}}!}\delta({N_{B}-N_{\bar{B}}-B})\,e^{(N_{B}+N_{\bar{B}})\,t}\Big{]}$
$\displaystyle=$
$\displaystyle\ln\Big{[}\int_{0}^{2\pi}\frac{d\phi}{2\,\pi}\,\exp\big{(}z_{B}\,e^{t+i\phi}\big{)}\,\exp\big{(}z_{\bar{B}}\,e^{t-i\phi}\big{)}\,e^{-i\,\phi\,B}\Big{]}$
$\displaystyle=$
$\displaystyle\ln\Big{[}\left(\frac{z_{B}}{z_{\bar{B}}}\right)^{B/2}\,I_{B}\big{(}2\,z\,e^{t}\big{)}\Big{]},$

where we set the auxiliary parameters $\lambda_{B},\lambda_{\bar{B}}$ to
unity. The $k$:th cumulant is then given by

$C_{k}=\frac{d^{k}G_{c}(t)}{d\,t^{k}}\mid_{t=0}.$ (21)

The generating functions for the corresponding factorial cumulants is obtained
by the replacement [55] $e^{t}\to x$,

$H_{c}(x)=\ln\Big{[}\left(\frac{z_{B}}{z_{\bar{B}}}\right)^{B/2}\,I_{B}\big{(}2\,z\,x\big{)}\Big{]}.$
(22)

The factorial cumulants of $\delta(N_{B}+N_{\bar{B}})$ are then given by

$F_{n}=\frac{d^{n}H_{c}(x)}{d\,x^{n}}\mid_{x=1}.$ (23)

We note that the factorial cumulants are identical to the functions
$f^{(n)}(z)$ defined in [34] and thus satisfy the recurrence relation

$F_{n+1}=z^{n+1}\frac{d}{d\,z}\left(F_{n}/z^{n}\right).$ (24)

The factorial cumulant $F_{n}$ of any order $n$ can be obtained by using (24)
and (II,18), starting from $F_{1}=C_{1}=S$.

The cumulants and factorial cumulants in the full system are connected via the
general relations

$F_{n}=\sum_{k=1}^{k}\,s(n,k)\,C_{k},$ (25)

and

$C_{n}=\sum_{k=1}^{n}\,S(n,k)\,F_{k},$ (26)

derived in C, where $s(n,k)$ and $S(n,k)$ are Stirling numbers of the first
and second kind, respectively.

##  IV Fluctuations in a subsystem

In this section we explore the fluctuations in a subsystem, $A$, of the
canonical system discussed so far. The part of the system that is not in $A$
we denote by $R$. While the net baryon number is conserved in the full system,
$A+R$, it is not in the subsystem, because net baryon number can be
transferred between $A$ and $R$. We obtain analytic results for the cumulants
and factorial cumulants of the baryon and antibaryon numbers in $A$. Moreover,
we derive general relations between these cumulants and those of the net
baryon number.

The cumulants of fluctuations in a subsystem are computed by differentiating
the corresponding generating functions. Thereby we employ the Faà di Bruno
formula [56, 57]

$\frac{d^{n}}{dx^{n}}f(h(x))=\sum_{k=1}^{n}f^{(k)}(h(x))\,B_{n,k}\left(h^{(1)}(x),h^{(2)}(x),\dots,h^{(n-k+1)}(x)\right),$
(27)

and generalizations thereof [56, 58]. In (27), $f^{(k)}$ and $h^{(k)}$ denote
the $k$th derivatives while $B_{n,k}(y_{1},y_{2},\dots)$ are partial Bell
polynomials [59].

###  IV.1 Factorial cumulants of $\delta N_{B}$ and $\delta N_{\bar{B}}$

Analytic expressions for the factorial cumulants of the baryon and antibaryon
numbers up to sixth order in a subsystem, taking the exact conservation of the
net baryon number in the total system into account, were presented in [15]. In
this section we derive closed-form expressions for these factorial cumulants.
Compact expressions are obtained in terms of the cumulants of the baryon
number of the total system, $c_{k}$.

The generating function for factorial cumulants is in this case [15]

$\displaystyle g_{fc}(x,\bar{x})$ $\displaystyle=$
$\displaystyle\frac{B}{2}\big{[}\ln\big{(}p(x)\big{)}-\ln\big{(}\bar{p}(\bar{x})\big{)}\big{]}$
$\displaystyle+$
$\displaystyle\ln\big{[}I_{B}\big{(}2\,z\,\sqrt{p(x)\,\bar{p}(\bar{x})}\big{)}\big{]},$

where $B$ is the total net baryon number, $p(x)=1-\alpha_{B}+\alpha_{B}\,x$
and $\bar{p}(\bar{x})=1-\alpha_{\bar{B}}+\alpha_{\bar{B}}\,\bar{x}$. Here
$\alpha_{B}$ is the probability for finding a baryon in the subsystem $A$ and
$\alpha_{\bar{B}}$ that for an antibaryon, while $z=\sqrt{z_{B}\,z_{\bar{B}}}$
is the geometric mean of the single-particle partition functions $z_{B}$ and
$z_{\bar{B}}$ (cf. Ref. [34]).

The factorial cumulants in the subsystem $A$ are obtained by evaluating the
derivatives with respect to $x$ and $\bar{x}$

$F^{(n,m)}=\frac{\partial^{n}}{\partial\,x^{n}}\,\frac{\partial^{m}}{\partial\bar{x}^{m}}\,g_{fc}(x,\bar{x})|_{x=\bar{x}=1}.$
(29)

The derviatives of the first two terms in (IV.1), $g_{fc,1}(x)$ and
$g_{fc,2}(\bar{x})$, are obtained by employing the Faà di Bruno formula for
functions of the form $\log(f(y))$, where $f(y)$ is a first order polynomial
in y [57]. One finds

$\frac{\partial^{n}}{\partial
x^{n}}\,g_{fc,1}(x)=-\frac{B}{2}\,(-\alpha_{B})^{n}\,(n-1)!,$ (30)

and

$\frac{\partial^{n}}{\partial\bar{x}^{n}}\,g_{fc,2}(\bar{x})=\frac{B}{2}\,(-\alpha_{\bar{B}})^{n}\,(n-1)!.$
(31)

For the last term in (IV.1),

$g_{fc,3}(x,\bar{x})=\ln\big{[}I_{B}\big{(}2\,z\,\sqrt{p(x)\,\bar{p}(\bar{x})}\big{)}\big{]},$
(32)

we consider the composite function $f_{3}(g_{3}(x,\bar{x}))$, where
$f_{3}(y)=\ln\big{[}I_{B}(2\,z\,e^{y})\big{]}$ and
$g_{3}(x,\bar{x})=\frac{1}{2}\,\big{(}\ln\big{[}p(x)\big{]}+\ln\big{[}\bar{p}(\bar{x})\big{]}\big{)}$.
The derivatives of $f_{3}(y)$ are obtained using (III) and (21), which yield
$f^{(n)}_{3}(0)=C_{n}$. Moreover, the non-zero derivatives of
$g_{3}(x,\bar{x})$ are given by

$\displaystyle g_{3}^{(n,0)}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(-\alpha_{B})^{n}\,(n-1)!\quad(n>0)$ (33)
$\displaystyle g_{3}^{(0,m)}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(-\alpha_{\bar{B}})^{m}\,(m-1)!\quad(m>0).$

Now, using the multivariate Faà di Bruno formula to compute the derivatives
with respect to $x$ and $\bar{x}$ of $g_{fc,3}(x,\bar{x})$ we find

$\frac{\partial^{n}}{\partial\,x^{n}}\,\frac{\partial^{m}}{\partial\bar{x}^{m}}\,g_{fc,3}(x,\bar{x})=\sum_{k=1}^{n+m}\,C_{k}\,B_{n,m;k}(\\{g_{3}^{(i,j)}\\})$
(34)

where the $B_{n,m;k}(\\{g_{3}^{(i,j)}\\})$ are multivariate Bell polynomials
(see A and Ref. [58]). As shown in B, the multivariate Bell polynomials are,
for the derivatives (33), given by

$B_{n,m;k}(\\{g_{3}^{(i,j)}\\})=\frac{(\alpha_{B})^{n}\,(\alpha_{\bar{B}})^{m}}{2^{k}}\,\,s(n,m;k),$
(35)

where $s(n,m;k)$ are generalized Stirling numbers of the first kind. They can
be expressed in terms of the standard Stirling numbers of the first kind,

$s(n,m;k)=\sum_{l=0}^{k}\,s(n,k-l)\,s(m,l).$ (36)

Now, collecting all terms and using (15), we find a closed-form expression for
the factorial cumulants of a subsystem, which accounts for baryon-number
conservation in full phase space,

$\displaystyle F^{(n,m)}$ $\displaystyle=$
$\displaystyle-\frac{B}{2}\Big{(}(-\alpha_{B})^{n}\,(n-1)!\,\delta_{m,0}-(-\alpha_{\bar{B}})^{m}\,(m-1)!\,\delta_{n,0}\Big{)}$
(37) $\displaystyle+$
$\displaystyle(\alpha_{B})^{n}\,(\alpha_{\bar{B}})^{m}\sum_{k=1}^{n+m}c_{k}\sum_{l=0}^{k}\,s(n,k-l)\,s(m,l).$

The first few factorial cumulants are

$\displaystyle F^{(1,0)}$ $\displaystyle=$ $\displaystyle\alpha_{B}\,\langle
N_{B}\rangle,$ $\displaystyle F^{(0,1)}$ $\displaystyle=$
$\displaystyle\alpha_{\bar{B}}\,\langle N_{\bar{B}}\rangle,$ $\displaystyle
F^{(2,0)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{2}\,\big{(}c_{2}-\langle N_{B}\rangle\big{)},$
$\displaystyle F^{(1,1)}$ $\displaystyle=$
$\displaystyle\alpha_{B}\,\alpha_{\bar{B}}\,c_{2},$ $\displaystyle F^{(0,2)}$
$\displaystyle=$ $\displaystyle(\alpha_{\bar{B}})^{2}\,\big{(}c_{2}-\langle
N_{\bar{B}}\rangle\big{)},$ (38) $\displaystyle F^{(3,0)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{3}\,\big{(}c_{3}-3\,c_{2}+2\,\langle
N_{B}\rangle\big{)},$ $\displaystyle F^{(2,1)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{2}\,\alpha_{\bar{B}}\,\big{(}c_{3}-c_{2}\big{)},$
$\displaystyle F^{(1,2)}$ $\displaystyle=$
$\displaystyle\alpha_{B}\,(\alpha_{\bar{B}})^{2}\,\big{(}c_{3}-c_{2}\big{)},$
$\displaystyle F^{(0,3)}$ $\displaystyle=$
$\displaystyle(\alpha_{\bar{B}})^{3}\,\big{(}c_{3}-3\,c_{2}+2\,\langle
N_{\bar{B}}\rangle\big{)}.$

###  IV.2 Cumulants of $\delta N_{B}$ and $\delta N_{\bar{B}}$

The generating function for cumulants of $\delta N_{B}$ and $\delta
N_{\bar{B}}$ in the subsystem is obtained from (IV.1) by the replacements
$x\to e^{t}$ and $\bar{x}\to e^{s}$,

$\displaystyle g_{c}(t,s)$ $\displaystyle=$
$\displaystyle\frac{B}{2}\big{[}\ln\big{(}q(t)\big{)}-\ln\big{(}\bar{q}(s)\big{)}\big{]}$
$\displaystyle+$
$\displaystyle\ln\big{[}I_{B}\big{(}2\,z\,\sqrt{q(t)\,\bar{q}(s)}\big{)}\big{]},$

where $q(t)=1-\alpha_{B}+\alpha_{B}\,e^{t}$ and
$\bar{q}(s)=1-\alpha_{\bar{B}}+\alpha_{\bar{B}}\,e^{s}$.

The cumulants in the subsystem $A$, defined by the acceptance probabilities
$\alpha_{B}$ and $\alpha_{\bar{B}}$, are given by

$C^{(n,m)}=\frac{\partial^{n}}{\partial t^{n}}\,\frac{\partial^{m}}{\partial
s^{m}}\,g_{c}(t,s)|_{t=s=0}.$ (40)

As shown in C, it follows from the relation between the generating functions,
$g_{c}(t,s)=g_{fc}(e^{t},e^{s})$, that the cumulants can be obtained from the
corresponding factorial cumulants using

$C^{(n,m)}=\sum_{k_{1}=0}^{n}\,\sum_{k_{2}=0}^{m}\,F^{(k1,k2)}\,S(n,k_{1})\,S(m,k_{2}).$
(41)

Now, inserting the expression for the factorial cumulants in the canonical
ensemble (37) in (41), we find

$\displaystyle C^{(n,m)}$ $\displaystyle=$
$\displaystyle\frac{B}{2}\,\Big{(}\kappa_{\rm
Ber}^{(n)}(\alpha_{B})\,\delta_{m,0}-\kappa_{\rm
Ber}^{(m)}(\alpha_{\bar{B}})\,\delta_{n,0}\Big{)}$ $\displaystyle+$
$\displaystyle\sum_{k=1}^{n+m}c_{k}\sum_{i=0}^{k}B_{n,k-i}\big{(}\kappa_{\rm
Ber}^{(1)}(\alpha_{B}),\kappa_{\rm Ber}^{(2)}(\alpha_{B}),\dots\big{)}$
$\displaystyle\times$ $\displaystyle B_{m,i}\big{(}\kappa_{\rm
Ber}^{(1)}(\alpha_{\bar{B}}),\kappa_{\rm
Ber}^{(2)}(\alpha_{\bar{B}}),\dots\big{)},$

where [34]

$\kappa_{\rm Ber}^{(n)}(\alpha)=\delta_{n,1}+(-1)^{1+n}\,{\rm
Li}_{1-n}(1-1/\alpha),$ (43)

are the cumulants of the Bernoulli distribution with success probability $p$
and ${\rm Li}_{n}(x)$ is the polylogarithm. In obtaining (IV.2) we used the
relations

$\displaystyle\sum_{k=1}^{n}S(n,k)\,(-\alpha)^{k}\,(k-1)!=-\kappa^{(n)}_{\rm
Ber}(\alpha),$ (44)
$\displaystyle\sum_{k=1}^{n}S(n,k)\,(\alpha)^{k}\,s(k,i)=B_{n,i}\big{(}\kappa_{\rm
Ber}^{(1)}(\alpha),\kappa_{\rm Ber}^{(2)}(\alpha),\dots\big{)}.$

We provide explicit expressions for the first few cumulants in the subsystem,

$\displaystyle C^{(1,0)}$ $\displaystyle=$ $\displaystyle\alpha_{B}\,\langle
N_{B}\rangle,$ $\displaystyle C^{(0,1)}$ $\displaystyle=$
$\displaystyle\alpha_{\bar{B}}\,\langle N_{\bar{B}}\rangle,$ $\displaystyle
C^{(2,0)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{2}\,c_{2}+\alpha_{B}(1-\alpha_{B})\,\langle
N_{B}\rangle,$ $\displaystyle C^{(1,1)}$ $\displaystyle=$
$\displaystyle\alpha_{B}\,\alpha_{\bar{B}}\,c_{2},$ $\displaystyle C^{(0,2)}$
$\displaystyle=$
$\displaystyle(\alpha_{\bar{B}})^{2}\,c_{2}+\alpha_{\bar{B}}(1-\alpha_{\bar{B}})\,\langle
N_{\bar{B}}\rangle,$ (45) $\displaystyle C^{(3,0)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{3}\,c_{3}+3(\alpha_{B})^{2}\,(1-\alpha_{B})\,c_{2}$
$\displaystyle+$
$\displaystyle\alpha_{B}(1-3\,\alpha_{B}+2\,(\alpha_{B})^{2})\,\langle
N_{B}\rangle,$ $\displaystyle C^{(2,1)}$ $\displaystyle=$
$\displaystyle(\alpha_{B})^{2}\,\alpha_{\bar{B}}\,c_{3}+\alpha_{B}\,\alpha_{\bar{B}}\,(1-\alpha_{B})\,c_{2},$
$\displaystyle C^{(1,2)}$ $\displaystyle=$
$\displaystyle\alpha_{B}\,(\alpha_{\bar{B}})^{2}\,c_{3}+\alpha_{B}\,\alpha_{\bar{B}}\,(1-\alpha_{\bar{B}})\,c_{2},$
$\displaystyle C^{(0,3)}$ $\displaystyle=$
$\displaystyle(\alpha_{\bar{B}})^{3}\,c_{3}+3(\alpha_{\bar{B}})^{2}\,(1-\alpha_{\bar{B}})\,c_{2}$
$\displaystyle+$
$\displaystyle\alpha_{\bar{B}}(1-3\,\alpha_{\bar{B}}+2\,(\alpha_{\bar{B}})^{2})\,\langle
N_{\bar{B}}\rangle.$

###  IV.3 Cumulants of the net baryon number

The generating function for net baryon number cumulants in a subsystem, is
obtained by making the substitutions $x\to e^{t}$ and $\bar{x}\to e^{-t}$ in
(IV.1). The generating function can also be derived starting from the
probability distribution $P_{A}(B_{A})$ of the net baryon number in the
subsystem [31, 34]. One finds333The generating function for cumulants of
$\delta(N_{B}-N_{\bar{B}})$ in a finite acceptance, $g_{net}$, reduces to the
generating function for cumulants of $\delta(N_{B}+N_{\bar{B}})$ in the full
system, Eq. (III), in the limit $\alpha_{B}\to 1$ and
$\alpha_{\bar{B}}\to-e^{t}$.

$\displaystyle g_{net}(t)$ $\displaystyle=$
$\displaystyle\ln\left(\sum_{B_{A}}P_{A}(B_{A})e^{B_{A}t}\right)$ (46)
$\displaystyle=$
$\displaystyle\frac{B}{2}\left[\ln(q_{1}(t))-\ln(q_{2}(t))\right]+\ln\left\\{I_{B}[2\,z\,\sqrt{q_{1}(t)\,q_{2}(t)}]\right\\},$

where $B$ is the total net baryon number,
$q_{1}(t)=1-\alpha_{B}+\alpha_{B}\,e^{t}$ and
$q_{2}(t)=1-\alpha_{\bar{B}}+\alpha_{\bar{B}}\,e^{-t}$.

In order to compute the net baryon number cumulants, we employ the formula of
Faà di Bruno (27) to compute derivatives of $g_{net}(t)$. For the first two
terms in (46), which we denote by $(a)$ and $(b)$, we choose $f(y)=\ln(y)$,
$h(t)=q_{1}(t)$ and $f(y)=-\ln(y)$, $h(t)=q_{2}(t)$, as in [34]. Employing
again the Faà di Bruno formula for functions of the form $\log(f(y))$ [57], we
find a closed-form expression for the corresponding contribution to the net
baryon number cumulants [34]

$\displaystyle\kappa^{(a+b)}_{n}$ $\displaystyle=$
$\displaystyle\frac{B}{2}\Big{(}\kappa_{\rm
Ber}^{(n)}(\alpha_{B})+(-1)^{(n+1)}\kappa_{\rm
Ber}^{(n)}(\alpha_{\bar{B}})\Big{)}$ (47) $\displaystyle\equiv$ $\displaystyle
B\,k^{(n)}_{+},$

where $\kappa_{\rm Ber}^{(n)}(\alpha_{B})$ is the $n$:th cumulant of the
Bernoulli distribution (43) and the second line defines $k^{(n)}_{+}$.

Now, in the evaluation of the derivatives of the last term in (46), denoted by
${(c)}$, we deviate slightly from [34], and choose
$f(y)=\ln\big{(}I_{B}[2\,z\,e^{y}]\big{)}$ and
$h(t)=\frac{1}{2}\big{(}\ln(q_{1}(t))+\ln(q_{2}(t))\big{)}$ and evaluate the
derivatives at $y=1$ and $t=0$, respectively. Using the fact that the
derivatives of $f(y)$ are given by the cumulants $C_{k}$ (cf. (III,21)) and
those of $h(t)$ by [34]

$h^{(n)}\equiv k_{-}^{(n)}=\frac{1}{2}\Big{(}\kappa_{\rm
Ber}^{(n)}(\alpha_{B})-(-1)^{(n+1)}\kappa_{\rm
Ber}^{(n)}(\alpha_{\bar{B}})\Big{)},$ (48)

we find that the corresponding contributions to the cumulants are given by

$\kappa^{(c)}_{n}=\sum_{k=1}^{n}\,C_{k}\,B_{n,k}(k_{-}^{(1)},\dots,k_{-}^{(n-k+1)}).$
(49)

Thus, the net baryon cumulants in a subsystem are closely related to the
cumulants of $\delta N^{tot}=\delta(N_{B}+N_{\bar{B}})$ in the full canonical
system, where the net baryon number is strictly conserved. Collecting all
terms, we find a compact form for the net baryon cumulants in a subsystem of a
canonical system,

$\kappa_{n}=B\,k_{+}^{(n)}+\sum_{k=1}^{n}\,C_{k}\,B_{n,k}(k_{-}^{(1)},\dots,k_{-}^{(n-k+1)}).$
(50)

We note that, using (26), the cumulants $\kappa_{n}$ can be expressed in terms
of the factorial cumulants $F_{k}$ rather than cumulants $C_{k}$. The
resulting analytic expression for $\kappa_{n}$ is identical to the one
obtained in [34] in terms of the functions $f^{(k)}(z)$, which as noted above,
are equal to the factorial cumulants (23).

In C, we obtain general relations of the net baryon cumulants to the baryon
and antibaryon cumulants in the subsystem,

$\kappa_{n}=\sum_{i=0}^{n}\,\binom{n}{i}\,(-1)^{n-i}\,C^{(i,n-i)},$ (51)

where $\binom{n}{i}$ is a binomial coefficient, and to the factorial ones,

$\kappa_{n}=\sum_{\begin{subarray}{c}k_{1},k_{2}\\\ k_{1}+k_{2}\geq
0\end{subarray}}^{n}F^{(k_{1},k_{2})}\,\sum_{i=0}^{n}\binom{n}{i}\,(-1)^{n-i}\,S(i,k_{1})\,S(n-i,k_{2}).$
(52)

One thus finds,

$\displaystyle\kappa_{1}$ $\displaystyle=$ $\displaystyle
F^{(1,0)}-F^{(0,1)},$ $\displaystyle\kappa_{2}$ $\displaystyle=$
$\displaystyle F^{(1,0)}+F^{(0,1)}+F^{(2,0)}+F^{(0,2)}-2\,F^{(1,1)},$
$\displaystyle\kappa_{3}$ $\displaystyle=$ $\displaystyle
F^{(1,0)}-F^{(0,1)}+3(F^{(2,0)}-F^{(0,2)})$ $\displaystyle-$ $\displaystyle
3(F^{(2,1)}-F^{(1,2)})+F^{(3,0)}-F^{(0,3)},$ $\displaystyle\kappa_{4}$
$\displaystyle=$ $\displaystyle
F^{(1,0)}+F^{(0,1)}+7(F^{(2,0)}+F^{(0,2)})+6(F^{(3,0)}+F^{(0,3)})$ (53)
$\displaystyle-$ $\displaystyle
6(F^{(2,1)}+F^{(1,2)})-4(F^{(3,1)}+F^{(1,3)})+6\,F^{(2,2)}-2\,F^{(1,1)}$
$\displaystyle+$ $\displaystyle F^{(4,0)}+F^{(0,4)},$
$\displaystyle\kappa_{5}$ $\displaystyle=$ $\displaystyle
F^{(1,0)}-F^{(0,1)}+15(F^{(2,0)}-F^{(0,2)})+25(F^{(3,0)}-F^{(0,3)})$
$\displaystyle+$ $\displaystyle
10(F^{(4,0)}-F^{(0,4)})-15(F^{(2,1)}-F^{(1,2)})-20(F^{(3,1)}-F^{(1,3)})$
$\displaystyle-$ $\displaystyle
5(F^{(4,1)}-F^{(1,4)})+10(F^{(3,2)}-F^{(2,3)})+F^{(5,0)}-F^{(0,5)}$
$\displaystyle\kappa_{6}$ $\displaystyle=$ $\displaystyle
F^{(1,0)}+F^{(0,1)}+31(F^{(2,0)}+F^{(0,2)})+90(F^{(3,0)}+F^{(0,3)})$
$\displaystyle+$ $\displaystyle
65(F^{(4,0)}+F^{(0,4)})+15(F^{(5,0)}+F^{(0,5)})-30(F^{(2,1)}+F^{(1,2)})$
$\displaystyle-$ $\displaystyle
80(F^{(3,1)}+F^{(1,3)})-45(F^{(4,1)}+F^{(1,4)})-6(F^{(5,1)}+F^{(1,5)})$
$\displaystyle+$ $\displaystyle
30(F^{(3,2)}+F^{(2,3)})+15(F^{(4,2)}+F^{(2,4)})-20F^{(3,3)}+30F^{(2,2)}$
$\displaystyle-$ $\displaystyle 2F^{(1,1)}+F^{(6,0)}+F^{(0,6)},$

When the factorial cumulants (37, IV.1) are plugged into (52,IV.3), one
recovers the explicit expressions given in [34]

##  V High- and low-energy limits

In the high-energy limit, ${\rm max}(B,1/2)<<2\,z$, one finds [34]

$F_{n}=2\,z\,\delta_{n,1}-(-1)^{n-1}\,(n-1)!\Big{(}\frac{1}{2}-n\frac{4\,B^{2}-1}{16\,z}\Big{)}+\mathcal{O}(1/z^{2}).$
(54)

Inserting this into (26), one finds the high-energy limit of the cumulants of
$\delta(N_{B}+N_{\bar{B}})$,

$C_{m}^{(HE)}=\sum_{n=1}^{m}\,S(m,n)\,\left(2\,z\,\delta_{n,1}-(-1)^{n-1}\,(n-1)!\Big{(}\frac{1}{2}-n\frac{4\,B^{2}-1}{16\,z}\Big{)}\right),$
(55)

where terms of order $\mathcal{O}(1/z^{2})$ and higher have been dropped.
Using the relations

$\displaystyle\sum_{n=1}^{m}\,S(m,n)\,(-1)^{n-1}\,(n-1)!$ $\displaystyle=$
$\displaystyle\delta_{m,1},$ (56)
$\displaystyle\sum_{n=1}^{m}\,S(m,n)\,(-1)^{n-1}\,n!$ $\displaystyle=$
$\displaystyle(-1)^{m-1},$ (57)

we find

$\displaystyle C_{1}^{(HE)}$ $\displaystyle=$ $\displaystyle
2\,z-\frac{1}{2}+\frac{4\,B^{2}-1}{16\,z},$ (58) $\displaystyle C_{n}^{(HE)}$
$\displaystyle=$ $\displaystyle
2\,z+(-1)^{n-1}\,\frac{4\,B^{2}-1}{16\,z}\qquad({\rm for}\,\,n>1).$

When the cumulants (58) are plugged into (50), one recovers the high-energy
limits of $\kappa_{n}$ derived in [34].

We note that for $n>1$ and $z\to\infty$, the cumulants

$C_{n}^{(HE)}=2\,z=C_{1}^{(HE)}+\frac{1}{2}$ (59)

are independent of $n$, like the cumulants of a Poisson distribution. Thus, in
the canonical ensemble, the high-energy limit of the cumulants of $N^{(tot)}$
are, except for the first one, equal to those of a Poisson distribution with
the cumulants $\langle N_{B}+N_{\bar{B}}\rangle+1/2$. Consequently, in the
high-energy limit ($B/z\to 0$), the fluctuations of
$\delta(N_{B}+N_{\bar{B}})$ in the canonical ensemble are, apart from a small
correction, the same as those in the grand-canonical ensemble.

In the low-energy limit, we expand the particle numbers about $z=0$. For
$B>0$, one finds444Compared to the result given in [34], we have here included
the $z^{4}$ terms.

$\displaystyle\langle N_{B}\rangle$ $\displaystyle=$ $\displaystyle
B+\frac{z^{2}}{B+1}-\frac{z^{4}}{(B+1)^{2}(B+2)}+\mathcal{O}(z^{6}/B^{5}),$
(60) $\displaystyle\langle N_{\bar{B}}\rangle$ $\displaystyle=$
$\displaystyle\frac{z^{2}}{B+1}-\frac{z^{4}}{(B+1)^{2}(B+2)}+\mathcal{O}(z^{6}/B^{5}).$
(61)

Thus, retaining terms of order $z^{4}$,

$C_{1}=F_{1}=\langle N_{B}+N_{\bar{B}}\rangle\approx
B+\frac{2\,z^{2}}{B+1}-\frac{2\,z^{4}}{(B+1)^{2}(B+2)}.$ (62)

Using the recurrence relation for $F_{n}$ (24), one finds for the low-energy
factorial cumulants

$\displaystyle F_{n}$ $\displaystyle=$
$\displaystyle(-1)^{n-1}\,(n-1)!\,B+(\delta_{n,1}+\delta_{n,2})\,\frac{2\,z^{2}}{B+1}$
$\displaystyle-$
$\displaystyle(\delta_{n,1}+3\,\delta_{n,2}+6\,\delta_{n,3}+6\,\delta_{n,4})\,\frac{2\,z^{4}}{(B+1)^{2}(B+2)}.$

Similarly, using the recurrence relation for $C_{k}$ (16), one finds 555Note
that both the high- and low-energy limits of the cumulants and factorial
cumulants satisfy the relations (25) and (26).

$C_{k}^{(LE)}=\frac{2^{k}\,z^{2}}{B+1}-\frac{2^{2\,k-1}\,z^{4}}{(B+1)^{2}(B+2)},\quad(k>1).$
(64)

For the cumulants of the baryon number one then finds

$\displaystyle c_{1}^{(LE)}$ $\displaystyle=$ $\displaystyle
B/2+\frac{z^{2}}{B+1}-\frac{z^{4}}{(B+1)^{2}(B+2)}=\frac{1}{2}\langle
N_{B}+N_{\bar{B}}\rangle,$ (65) $\displaystyle c_{k}^{(LE)}$ $\displaystyle=$
$\displaystyle\frac{z^{2}}{B+1}-\frac{2^{k-1}\,z^{4}}{(B+1)^{2}(B+2)}=\langle
N_{\bar{B}}\rangle+\mathcal{O}(z^{4}),\quad(k>1).$ (66)

Thus, in the low-energy limit $(z\to 0)$, the fluctuations of the baryon
number, $c_{k}=\langle\left(\delta N_{B}\right)^{k}\rangle_{c}\,\,\,(k>1)$,
approach the cumulants of a Poisson distribution with expectation value
$\langle N_{\bar{B}}\rangle$. In other words, in the canonical ensemble,
fluctuations of the baryon number are at low energies equal to the
fluctuations of the antibaryon number in the grand canonical ensemble.

The cumulants of $\delta(N_{B}+N_{\bar{B}})$, (14), are shown in Fig. 1 as
functions of $z$ for net baryon number $B=350$.

Figure 1: The cumulants of $\delta(N_{B}+N_{\bar{B}})$, $C_{n}$, normalized to
$N_{B}+N_{\bar{B}}=C_{1}$. The corresponding beam energy increases with
increasing $z$.

In Fig. 2 the corresponding cumulants of $\delta N_{B}$, (11), are shown. The
behaviour in the high- and low-energy limits (large and small $z$,
respectively) of the cumulants shown in Figs. 1 and 2 is consistent with the
analytic results (58) and (64).

Figure 2: The cumulants of $\delta(N_{B})$, $c_{n}$, normalized to
$N_{\bar{B}}$.

##  VI Summary and conclusions

The main objective of these studies was to calculate fluctuations of the
multiplicity distributions of baryons, antibaryons, as well as their sum and
difference in a canonical ensemble. We have focused on cumulants and factorial
cumulants, assuming that baryons are correlated by exact conservation of the
net baryon number in the full system. To this end we have adopted a
statistical model of baryons and antibaryons in the spirit of the S-matrix
approach, where to leading order in the fugacity expansion, the partition
function takes the form of an ideal gas.

We have derived analytic expressions for the cumulants and factorial cumulants
of the baryon and antibaryon numbers in the canonical ensemble. Moreover,
based on the generalized Faà di Bruno formula for computing the derivatives of
composite functions, we obtained closed-form results for cumulants and
factorial cumulants of arbitrary order in a subsystem constrained by exact
conservation of the net baryon number in the full system. We also presented
general relations, independent of assumptions on the partition function,
between cumulants and factorial cumulants as well as between net-baryon-number
fluctuations in a subsystem and those of the baryon and antibaryon numbers of
the full system. Finally, we deduced the asymptotic forms of the cumulants and
factorial cumulants of baryon and antibaryon multiplicities in the high- and
low-energy limits.

The analytic results for fluctuation observables of the baryon and antibaryon
numbers in a subsystem, given a canonical partition function for the full
system, as well as the general relations between observables, provide a firm
baseline for the phenomenological interpretation of fluctuation data obtained
in heavy-ion collisions. Applications of these results will be presented in an
upcoming publication.

##  VII Acknowledgements

We acknowledge stimulating discussions with Peter Braun-Munzinger, Anar
Rustamov and Johanna Stachel. K.R. also acknowledges the supported by the
Polish National Science Center (NCN) under OPUS Grant No. 2018/31/B/ST2/01663,
and by the Polish Ministry of Science.

##  Appendix A Multivariate Bell polynomials

In order to express the derivatives of a composite function of several
variables, like $g_{fc}(x,\bar{x})$ in (IV.1), in a form analogous to the Faà
di Bruno formula one needs a generalization of the Bell polynomials [59].

The multivariate partial Bell polynomials can be obtained from a
generalization of the generating function for the standard partial Bell
polynomials [58]

$\Phi(t,u)=\exp\left(u\sum_{j}x^{(j)}\frac{t^{j}}{j!}\right)=\sum_{n\geq k\geq
0}B_{n,k}(x^{(1)},x^{(2)},\dots,x^{(n-k+1)})\,\frac{t^{n}}{n!}\,u^{k}.$ (67)

The regular partial Bell polynomials are given by

$B_{n,k}(x^{(1)},x^{(2)},\dots,x^{(n-k+1)})=\frac{1}{k!}\frac{\partial^{n}}{\partial
t^{n}}\,\frac{\partial^{k}}{\partial u^{k}}\,\Phi(t,u)|_{t=u=0}.$ (68)

The Bell polynomials enter the Faà di Bruno formula (27) for the derivatives
of a composite function $f(g(z))$, Below we present the generalizations of
(27) to composite functions of the form $f(g(z_{1},z_{2}))$ and
$f(g_{1}(z),g_{2}(z))$.

In Ref. [58] the generating function for the general case of arbitrary
dimensions of $z$ and $g$ is given. In order to keep the notation transparent,
we give the generating function for two-dimensional $\vec{z}=(z_{1},z_{2})$
and $\vec{g}=(g_{1},g_{2})$,

$\displaystyle\Phi(t_{1},t_{2},u_{1},u_{2})$ $\displaystyle=$
$\displaystyle\exp\left(\sum_{\begin{subarray}{c}j_{1},j_{2}\\\
j_{1}+j_{2}>0\end{subarray}}^{\infty}\big{[}u_{1}*x_{1}^{(j_{1},j_{2})}+u_{2}*x_{2}^{(j_{1},j_{2})}\big{]}\frac{t_{1}^{j_{1}}\,t_{2}^{j_{2}}}{j_{1}!\,j_{2}!}\right)$
$\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}n_{1},n_{2},k_{1},k_{2}\\\
n_{1}+n_{2}\geq k_{1}+k_{2}\geq
0\end{subarray}}B_{\\{n_{1},n_{2}\\},\\{k_{1},k_{2}\\}}\big{(}\\{\vec{x}^{(i,j)}\\}\big{)}\,\frac{t_{1}^{n_{1}}\,t_{2}^{n_{2}}}{n_{1}!\,n_{2}!}\,(u_{1}^{k_{1}}+u_{2}^{k_{2}}).$

The generalized Bell polynomials are then given by

$\displaystyle B_{n_{1},n_{2};k_{1},k_{2}}(\\{\vec{x}^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle\frac{1}{k_{1}!\,k_{2}!}$ $\displaystyle\times$
$\displaystyle\frac{\partial^{n_{1}}}{\partial
t_{1}^{n_{1}}}\frac{\partial^{n_{2}}}{\partial
t_{2}^{n_{2}}}\frac{\partial^{k_{1}}}{\partial
u_{1}^{k_{1}}}\frac{\partial^{k_{2}}}{\partial
u_{2}^{k_{2}}}\,\Phi(t_{1},t_{2},u_{1},u_{2})|_{t_{1}=t_{2}=u_{1}=u_{2}=0}.$

Now, for computing the factorial cumulants, we need the case with vector
argument $\vec{z}=(z_{1},z_{2})$ and a scalar function $g(z_{1},z_{2})$,

$\displaystyle\Phi(t_{1},t_{2},u)$ $\displaystyle=$
$\displaystyle\exp\left(u\,\sum_{\begin{subarray}{c}j_{1},j_{2}\\\
j_{1}+j_{2}>0\end{subarray}}^{\infty}x^{(j_{1},j_{2})}\frac{t_{1}^{j_{1}}\,t_{2}^{j_{2}}}{j_{1}!\,j_{2}!}\right)$
$\displaystyle=$ $\displaystyle\sum_{\begin{subarray}{c}n_{1},n_{2},k\\\
n_{1}+n_{2}\geq k\geq
0\end{subarray}}B_{n_{1},n_{2};k}\big{(}\\{x^{(i,j)}\\}\big{)}\,\frac{t_{1}^{n_{1}}\,t_{2}^{n_{2}}}{n_{1}!\,n_{2}!}\,u^{k}.$

The corresponding partial Bell polynomials are given by,

$\displaystyle B_{n_{1},n_{2};k}(\\{x^{(i,j)}\\})$ $\displaystyle=$
$\displaystyle\frac{1}{k!}\,\frac{\partial^{n_{1}}}{\partial
t_{1}^{n_{1}}}\frac{\partial^{n_{2}}}{\partial
t_{2}^{n_{2}}}\frac{\partial^{k}}{\partial
u^{k}}\,\Phi(t_{1},t_{2},u)|_{t_{1}=t_{2}=u=0}.$ (72)

The first few generalized Bell polynomials are

$\displaystyle B_{1,0;1}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
x^{(1,0)}\,,$ $\displaystyle B_{0,1;1}(\\{x^{(i,j)}\\})$ $\displaystyle=$
$\displaystyle x^{(0,1)}\,,$ $\displaystyle B_{2,0;1}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle x^{(2,0)}\,,$ $\displaystyle
B_{2,0;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle(x^{(1,0)})^{2}\,,$
(73) $\displaystyle B_{0,2;1}(\\{x^{(i,j)}\\})$ $\displaystyle=$
$\displaystyle x^{(0,2)}\,,$ $\displaystyle B_{0,2;2}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle(x^{(0,1)})^{2}\,,$ $\displaystyle
B_{1,1;1}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle x^{(1,1)}\,,$
$\displaystyle B_{1,1;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
x^{(1,0)}\,x^{(0,1)}\,,$ $\displaystyle B_{3,0;1}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle x^{(3,0)}\,,$ $\displaystyle
B_{3,0;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
3\,x^{(1,0)}\,x^{(2,0)}\,,$ $\displaystyle B_{3,0;3}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle(x^{(1,0)})^{3}\,,$ $\displaystyle
B_{2,1;1}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle x^{(2,1)}\,,$
$\displaystyle B_{2,1;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
2\,x^{(1,0)}\,x^{(1,1)}+x^{(0,1)}\,x^{(2,0)},$ $\displaystyle
B_{2,1;3}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
x^{(0,1)}\,(x^{(1,0)})^{2}.$ $\displaystyle B_{1,2;1}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle x^{(1,2)}\,,$ $\displaystyle
B_{1,2;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
2\,x^{(0,1)}\,x^{(1,1)}+x^{(1,0)}\,x^{(0,2)},$ $\displaystyle
B_{1,2;3}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
x^{(1,0)}\,(x^{(0,1)})^{2}.$ $\displaystyle B_{0,3;1}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle x^{(0,3)}\,,$ $\displaystyle
B_{0,3;2}(\\{x^{(i,j)}\\})$ $\displaystyle=$ $\displaystyle
3\,x^{(0,1)}\,x^{(0,2)}\,,$ $\displaystyle B_{0,3;3}(\\{x^{(i,j)}\\})$
$\displaystyle=$ $\displaystyle(x^{(0,1)})^{3}\,,$

A combinatorial interpretation of the Bell polynomial
$B_{n,m;k}({x^{\\{i,j\\}}})$ goes as follows. Consider a collection of $n$
blue beads and $m$ red ones. How can these be split into $k$ groups is encoded
in the Bell polynomials. A group consisting of $i$ blue beads and $j$ red ones
is denoted by $x^{\\{i,j\\}}$. For instance, a system consisting of two blue
beads and one red one can be split into two groups in three ways. Two with one
blue and one red bead in one group and the remaining blue bead in the other
group and one with the two blue beads in one group and the red one in the
other. This case corresponds to the Bell polynomial
$B_{2,1;2}(\\{x^{(i,j)}\\})$ in (A).

The generalized Faà di Bruno formula for computing the derivatives of a
composite function of two variables, $f(g(z_{1},z_{2}))$, needed for the
calculation of the factorials cumulants is then,

$\displaystyle\frac{\partial^{n_{1}}}{\partial
z_{1}^{n_{1}}}\,\frac{\partial^{n_{2}}}{\partial
z_{2}^{n_{2}}}f(g(z_{1},z_{2}))=\sum_{k=1}^{n_{1}+n_{2}}f^{(k)}\,B_{n_{1},n_{2};k}(\\{g^{(i,j)}\\}).$
(74)

The first few terms are

$\displaystyle\frac{\partial}{\partial z_{1}}\,f(g(z_{1},z_{2}))$
$\displaystyle=$ $\displaystyle f^{(1)}g^{(1,0)}\,,$
$\displaystyle\frac{\partial}{\partial z_{2}}\,f(g(z_{1},z_{2}))$
$\displaystyle=$ $\displaystyle f^{(1)}g^{(0,1)}\,,$
$\displaystyle\frac{\partial^{2}}{\partial z_{1}^{2}}\,f(g(z_{1},z_{2}))$
$\displaystyle=$ $\displaystyle
f^{(1)}g^{(2,0)}+f^{(2)}\big{(}g^{(1,0)}\big{)}^{2}\,,$ (75)
$\displaystyle\frac{\partial^{2}}{\partial z_{2}^{2}}\,f(g(z_{1},z_{2}))$
$\displaystyle=$ $\displaystyle
f^{(1)}g^{(0,2)}+f^{(2)}\big{(}g^{(0,1)}\big{)}^{2}\,,$
$\displaystyle\frac{\partial}{\partial z_{1}}\frac{\partial}{\partial
z_{2}}\,f(g(z_{1},z_{2}))$ $\displaystyle=$ $\displaystyle
f^{(1)}g^{(1,1)}+f^{(2)}\,g^{(1,0)}\,g^{(0,1)}\,.$

We also need the Bell polynomials for computing the derivatives of a function
of the form $f(g_{1}(z),g_{2}(z))$. The generating function is:

$\displaystyle\Phi(t,u_{1},u_{2})$ $\displaystyle=$
$\displaystyle\exp\left(\sum_{j=1}^{\infty}\big{[}u_{1}*x_{1}^{(j)}+u_{2}*x_{2}^{(j)}\big{]}\frac{t^{j}}{j!}\right)$
$\displaystyle=$ $\displaystyle\sum_{\begin{subarray}{c}n,k_{1},k_{2}\\\ n\geq
k_{1}+k_{2}\geq
0\end{subarray}}B_{n;k_{1},k_{2}}\big{(}\\{\vec{x}^{(i)}\\}\big{)}\,\frac{t^{n}}{n!}\,(u_{1}^{k_{1}}+u_{2}^{k_{2}}).$

and the corresponding generalized Bell polynomials are given by,

$B_{n;k_{1},k_{2}}(\\{\vec{x}^{(i)}\\})=\frac{1}{k_{1}!\,k_{2}!}\frac{\partial^{n}}{\partial
t^{n}}\frac{\partial^{k_{1}}}{\partial
u_{1}^{k_{1}}}\frac{\partial^{k_{2}}}{\partial
u_{2}^{k_{2}}}\,\Phi(t_{1},t_{2},u)|_{t_{1}=t_{2}=u=0}.$ (77)

We note that the multivariate Bell polynomials of this type can be constructed
from the standard (univariate) Bell polynomials, [56]

$\displaystyle
B_{n;k_{1},k_{2}}(x_{1}^{(1)},x_{1}^{(2)},\dots;x_{2}^{(1)},x_{2}^{(2)},\dots)$
(78)
$\displaystyle=\sum_{i=0}^{n}\binom{n}{i}\,B_{i,k_{1}}(x_{1}^{(1)},x_{1}^{(2)},\dots)\,B_{n-i,k2}(x_{2}^{(1)},x_{2}^{(2)},\dots).$

The first few Bell polynomials are

$\displaystyle B_{1;1,0}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$ $\displaystyle
x_{1}^{(1)},$ $\displaystyle B_{1;0,1}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$
$\displaystyle x_{2}^{(1)},$ $\displaystyle B_{2;1,0}(\\{\vec{x}^{(i)}\\})$
$\displaystyle=$ $\displaystyle x_{1}^{(2)},$ $\displaystyle
B_{2;0,1}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$ $\displaystyle x_{2}^{(2)},$
(79) $\displaystyle B_{2;2,0}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$
$\displaystyle(x_{1}^{(1)})^{2},$ $\displaystyle
B_{2;0,2}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$
$\displaystyle(x_{2}^{(1)})^{2},$ $\displaystyle
B_{2;1,1}(\\{\vec{x}^{(i)}\\})$ $\displaystyle=$ $\displaystyle
2\,x_{1}^{(1)}\,x_{2}^{(1)},$

and the corresponding Faà di Bruno formula is [56]

$\displaystyle\frac{\partial^{n}}{\partial
z^{n}}f(g_{1}(z),g_{2}(z))=\sum_{\begin{subarray}{c}k_{1},k_{2}\\\
k_{1}+k_{2}\geq
0\end{subarray}}^{n}f^{(k_{1},k_{2})}\,B_{n;k_{1},k_{2}}(\\{g_{1}^{(i)};g_{2}^{(i)}\\}).$
(80)

Here the first few terms are given by

$\displaystyle\frac{\partial}{\partial z}f(g_{1}(z),g_{2}(z))$
$\displaystyle=$ $\displaystyle f^{(1,0)}g_{1}^{(1)}+f^{(0,1)}g_{2}^{(1)}\,,$
$\displaystyle\frac{\partial^{2}}{\partial z^{2}}f(g_{1}(z),g_{2}(z))$
$\displaystyle=$ $\displaystyle
f^{(1,0)}g_{1}^{(2)}+f^{(0,1)}g_{2}^{(2)}+f^{(2,0)}\big{(}g_{1}^{(1)}\big{)}^{2}$
$\displaystyle+$ $\displaystyle
f^{(0,2)}\big{(}g_{2}^{(1)}\big{)}^{2}+2\,f^{(1,1)}\,g_{1}^{(1)}\,g_{2}^{(1)}\,.$

##  Appendix B Generalized Stirling numbers

Consider multivariate Bell polynomials of the type (72), with the derivatives
$x^{(i,j)}$ given by $g^{(i,j)}$ in (33). In this case, the generating
function (A) reduces to

$\Phi(t_{1},t_{2},u)=\big{[}(1+\alpha_{B}\,t_{1})\,(1+\alpha_{\bar{B}}\,t_{2})\big{]}^{u/2},$
(82)

and the corresponding Bell polynomials are

$B_{n,m;k}(\\{g_{3}^{(i,j)}\\})=\frac{(\alpha_{B})^{n}\,(\alpha_{\bar{B}})^{m}}{2^{k}}\,\,s(n,m;k),$
(83)

where

$s(n,m;k)=\frac{1}{k!}\,\frac{\partial^{n}}{\partial
z^{n}}\,\frac{\partial^{m}}{\partial w^{m}}\,\frac{\partial^{k}}{\partial
u^{k}}\,\Psi(z,w,u)|_{z=w=u=0},$ (84)

are generalized Stirling numbers and

$\Psi(z,w,u)=\big{[}(1+z)\,(1+w)\big{]}^{u}$ (85)

is the corresponding generating function.

The generating function for the regular Stirling numbers of the first kind

$s(n,k)=\frac{1}{k!}\,\frac{\partial^{n}}{\partial
z^{n}}\,\frac{\partial^{k}}{\partial u^{k}}\,\psi(z,w,u)|_{z=u=0}$ (86)

is of the form [57]

$\psi(z,u)=\big{[}1+z\big{]}^{u}.$ (87)

It follows that

$s(n,0;k)=s(0,m;k)=s(n,k),$ (88)

that $s(n,m;0)=0$, except for $n=m=0$, and that $s(n,m;k)=0$ for $k>n+m$ and
that

$s(n,m;k)=\sum_{l=0}^{k}\,s(n,k-l)\,s(m,l).$ (89)

Moreover, the generalized Stirling numbers satisfy the recurrence relations

$\displaystyle s(n,m;k)$ $\displaystyle=$ $\displaystyle
s(n-1,m;k-1)-(n-1)\,s(n-1,m;k),$ (90) $\displaystyle s(n,m;k)$
$\displaystyle=$ $\displaystyle s(n,m-1;k-1)-(m-1)\,s(n,m-1;k),$ (91)

in close analogy to the one obeyed by the regular Stirling numbers of the
first kind [57]. The relation (90) holds for $n,k\geq 1$ and (91) for $m,k\geq
1$. The two relations (90) and (91) combined yield a recurrence relation at
fixed $k$,

$s(n,m;k)=s(n+1,m-1;k)+(n-m+1)\,s(n,m-1;k).$ (92)

The recurrence relation for factorial cumulants,

$F^{(n+1,m)}=\frac{\alpha_{B}}{\alpha_{\bar{B}}}\,F^{(n,m+1)}-(n-m)\,\alpha_{B}\,F^{(n,m)},$
(93)

which was empirically deduced in [15], follows from (92).

The generalized Stirling numbers also have a combinatorial interpretation. The
absolute value, $|s(n,m;k)|$, equals the number of permutations of $n$ blue
beads and $m$ red beads in k disjoint cycles, where each cycle consists of
only blue or only red beads. The phase of $s(n,m;k)$ reproduces the phase
stemming from the derivatives of $g_{3}(x,\bar{x})$, given in (33).

##  Appendix C General relations between cumulants and factorial cumulants

In this appendix we derive general relations between cumulants and factorial
cumulants, which are independent of the assumed partition function. We start
by considering the fluctuations of $N^{(tot)}=N_{B}+N_{\bar{B}}$. The
cumulants $C_{n}$ and factorial cumulants $F_{n}$ are obtained by
differentiating the corresponding generating functions,

$C_{n}=\frac{d^{n}G(t)}{d\,t^{n}}\mid_{t=0}.$ (94)

and

$F_{n}=\frac{d^{n}H(x)}{d\,x^{n}}\mid_{x=1},$ (95)

which are closely related through [55]

$G(t)=H(e^{t}),\qquad H(x)=G(\ln(x)).$ (96)

By applying the Faà di Bruno formula to composite functions of the form $f(\ln
x)$ and $f(e^{t})$ one finds the following general relations between cumulants
and factorial cumulants [57],

$F_{n}=\sum_{k=1}^{k}\,s(n,k)\,C_{k},\qquad
C_{n}=\sum_{k=1}^{n}\,S(n,k)\,F_{k},$ (97)

where $s(n,k)$ and $S(n,k)$ are Stirling numbers of the first and second kind,
respectively. The Stirling numbers satisfy the orthogonality relations [57]

$\displaystyle\sum_{k=0}^{n}\,S(n,k)\,s(k,m)$ $\displaystyle=$
$\displaystyle\delta_{n,m},$ (98)
$\displaystyle\sum_{k=0}^{n}\,s(n,k)\,S(k,m)$ $\displaystyle=$
$\displaystyle\delta_{n,m}.$ (99)

Analogously, one finds the relation between the cumulants and factorial
cumulants of a subsystem by using the multivariant Faá di Bruno formula (74)
for a function of the form $G(e^{t},e^{s})$,

$C^{(n,m)}=\sum_{k_{1}=0}^{n}\,\sum_{k_{2}=0}^{m}\,F^{(k1,k2)}\,S(n,k_{1})\,S(m,k_{2}).$
(100)

For completeness we note that the relation (100) can be inverted, using the
orthogonality property of the Stirling numbers, (99),

$F^{(p,q)}=\sum_{n=0}^{p}\,\sum_{m=0}^{q}\,C^{(n,m)}\,s(p,n)\,s(q,m).$ (101)

Utilizing the fact that cumulants of the net baryon number correspond to
cumulants of $\delta N_{B}-\delta N_{\bar{B}}$, one finds, using the binomial
theorem, the general relation

$\kappa_{n}=\sum_{i=0}^{n}\,\binom{n}{i}\,(-1)^{n-i}\,C^{(i,n-i)},$ (102)

where $\binom{n}{i}$ is a binomial coefficient. Similarly, one obtains the
cumulants $\kappa^{s}_{n}$ of $\delta(N_{B}+N_{\bar{B}})$ in the subsystem
using

$\kappa^{s}_{n}=\sum_{i=0}^{n}\,\binom{n}{i}\,C^{(i,n-i)}.$ (103)

In the limit $\alpha_{B},\alpha_{\bar{B}}\to 1$, $\kappa_{n}\to 0$, while
$\kappa^{s}_{n}\to C_{n}$.

Finally, by inserting (100) in (102) one finds a general relation between the
factorial cumulants and the net-baryon-number cumulants [60],

$\kappa_{n}=\sum_{\begin{subarray}{c}k_{1},k_{2}\\\ k_{1}+k_{2}\geq
0\end{subarray}}^{n}F^{(k_{1},k_{2})}\,\sum_{i=0}^{n}\binom{n}{i}\,(-1)^{n-i}\,S(i,k_{1})\,S(n-i,k_{2}).$
(104)

## References

  * [1] M. M. Aggarwal et al. [STAR Collaboration], arXiv:1007.2613 [nucl-ex].
  * [2] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998).
  * [3] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999).
  * [4] M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989).
  * [5] B. Friman, F. Karsch, K. Redlich and V. Skokov, Eur. Phys. J. C 71, 1694 (2011).
  * [6] S. Ejiri, C. R. Allton, M. Doring, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Nucl. Phys. Proc. Suppl. 140, 505 (2005).
  * [7] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Phys. Rev. D 71, 054508 (2005).
  * [8] S. Ejiri, F. Karsch and K. Redlich, Phys. Lett. B 633, 275 (2006).
  * [9] F. Karsch, S. Ejiri and K. Redlich, Nucl. Phys. A 774, 619 (2006).
  * [10] F. Karsch and K. Redlich, Phys. Lett. B 695, 136 (2011).
  * [11] C. Sasaki, B. Friman and K. Redlich, Phys. Rev. Lett. 99, 232301 (2007).
  * [12] C. Sasaki, B. Friman and K. Redlich, Phys. Rev. D 75, 074013 (2007).
  * [13] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov and N. Xu, Phys. Rept. 853, 1 (2020).
  * [14] V. A. Kuznietsov, O. Savchuk, M. I. Gorenstein, V. Koch and V. Vovchenko, bPhys. Rev. C 105 no.4, 044903 (2022).
  * [15] M. Barej and A. Bzdak, Phys. Rev. C 102, 064908 (2020).
  * [16] A. Bialas, Phys. Rev. C 75, 024904 (2007).
  * [17] A. Bzdak and V. Koch, Phys. Rev. C 100, no.5, 051902 (2019).
  * [18] M. Barej and A. Bzdak, [arXiv:2205.05497 [hep-ph]].
  * [19] L. Adamczyk, et al. [STAR Collaboration], Phys. Rev. Lett. 112, 032302 (2014).
  * [20] X. Luo [STAR], PoS CPOD2014, 019 (2015).
  * [21] X. Luo, Nucl. Phys. A 956, 75-82 (2016).
  * [22] M. Abdallah et al. [STAR], Phys. Rev. C 104, no.2, 024902 (2021).
  * [23] A. Rustamov for the ALICE Collaboration, Nucl. Phys. A 967, 453 (2017).
  * [24] M. Arslandok for the ALICE Collaboration, Nucl. Phys. A 1005, 121979 (2021).
  * [25] S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 807, 135564, (2020).
  * [26] X. Luo, J. Xu, B. Mohanty and N. Xu, J. Phys. G 40, 105104 (2013).
  * [27] V. Skokov, B. Friman and K. Redlich, Phys. Rev. C 88, 034911 (2013).
  * [28] P. Braun-Munzinger, A. Rustamov and J. Stachel, Nucl. Phys. A 960, 114 (2017).
  * [29] T. Sugiura, T. Nonaka and S. Esumi, Phys. Rev. C 100, 044904 (2019).
  * [30] P. Braun-Munzinger, A. Rustamov and J. Stachel, Nucl. Phys. A 982, 307 (2019).
  * [31] A. Bzdak, V. Koch and V. Skokov, Phys. Rev. C 87, 014901 (2013).
  * [32] P. Braun-Munzinger, A. Rustamov and J. Stachel, arXiv:1907.03032 [nucl-th].
  * [33] V. Begun, M. Gazdzicki, M. I. Gorenstein and O. Zozulya, Phys. Rev. C 70, 034901 (2004).
  * [34] P. Braun-Munzinger, B. Friman, K. Redlich, A. Rustamov and J. Stachel, Nucl. Phys. A 1008, 122141 (2021).
  * [35] V. Vovchenko, O. Savchuk, R. V. Poberezhnyuk, M. I. Gorenstein and V. Koch, Phys. Lett. B 811, 135868 (2020).
  * [36] V. Vovchenko, R. V. Poberezhnyuk and V. Koch, JHEP 10, 089 (2020).
  * [37] C. A. Pruneau, Phys. Rev. C 100, no.3, 034905 (2019).
  * [38] R. Dashen, S. K. Ma and H. J. Bernstein, Phys. Rev. 187, 345-370 (1969).
  * [39] R. Venugopalan and M. Prakash, Nucl. Phys. A 546, 718-760 (1992).
  * [40] W. Weinhold, B. Friman and W. Norenberg, Phys. Lett. B 433, 236-242 (1998).
  * [41] F. Giacosa, V. Begun and W. Broniowski, Acta Phys. Polon. Supp. 9, 213 (2016).
  * [42] P. M. Lo, Eur. Phys. J. C 77, no.8, 533 (2017).
  * [43] A. Dash, S. Samanta and B. Mohanty, Phys. Rev. C 99, no.4, 044919 (2019).
  * [44] B. Friman, P. M. Lo, M. Marczenko, K. Redlich and C. Sasaki, Phys. Rev. D 92, no.7, 074003 (2015).
  * [45] P. M. Lo, Phys. Rev. D 102, no.3, 034038 (2020).
  * [46] P. Braun-Munzinger, K. Redlich and J. Stachel, doi:10.1142/9789812795533_0008 [arXiv:nucl-th/0304013 [nucl-th]].
  * [47] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Nature 561, no.7723, 321-330 (2018), and references therein.
  * [48] J. Cleymans, P. M. Lo, K. Redlich and N. Sharma, Phys. Rev. C 103, no.1, 014904 (2021).
  * [49] A. Andronic, P. Braun-Munzinger, B. Friman, P. M. Lo, K. Redlich and J. Stachel, Phys. Lett. B 792, 304-309 (2019).
  * [50] A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno and P. Petreczky, et al. Phys. Rev. D 95, no.5, 054504 (2017).
  * [51] J. Noronha-Hostler, P. Parotto, C. Ratti and J. M. Stafford, Phys. Rev. C 100, no.6, 064910 (2019).
  * [52] P. M. Lo, B. Friman, K. Redlich and C. Sasaki, Phys. Lett. B 778, 454 (2018).
  * [53] G. Almasi, A. Andronic, P. Braun-Munzinger, B. Friman, P. M. Lo, K. Redlich, C. Sasaki and J. Stachel, Nucl. Phys. A 982, 295 (2019).
  * [54] J. Goswami, F. Karsch, C. Schmidt, S. Mukherjee and P. Petreczky, Acta Phys. Polon. Supp. 14, 251 (2021).
  * [55] M. Kitazawa and X. Luo, Phys. Rev. C 96, 024910 (2017).
  * [56] J. Riordan, Derivatives of composite functions, Bull. Amer. Math. Soc. 52, 664 (1946).
  * [57] L. Comtet, Advanced Combinatorics (D. Riedel, Dordrecht, (1974).
  * [58] A. Schumann, Multivariate Bell Polynomials and Derivatives of Composed Functions, arXiv:1903.03899 [math.CA]
  * [59] E. T. Bell, Annals of Mathematics, Second Series 29, 38 (1927).
  * [60] X. Luo, Phys. Rev. C 91, 034907 (2015); Erratum: Phys. Rev. C 94, 059901 (2016),
  * [61] X. Pan, M. Xu and Y. Wu, Int. J. Mod. Phys. E 30 (2021) no.05, 2150036 doi:10.1142/S0218301321500361 [arXiv:2101.02822 [nucl-th]].


# Diversity of Radial Spin Textures in Chiral Materials

Daniel Gosálbez-Martínez<EMAIL_ADDRESS>Departamento de Física
Aplicada, Universidad de Alicante, 03690 Alicante, Spain  Institute of
Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland  National Centre for Computational Design and Discovery of Novel
Materials MARVEL, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Alberto Crepaldi  Dipartimento di Fisica, Politecnico
di Milano, Piazza Leonardo da Vinci 32, Milan 20133, Italy  Oleg V. Yazyev
<EMAIL_ADDRESS>Institute of Physics, Ecole Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  National Centre for
Computational Design and Discovery of Novel Materials MARVEL, Ecole
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

###### Abstract

We introduce a classification of the radial spin textures in momentum space
that emerge at high-symmetry points in crystals characterized by non-polar
chiral point groups ($D_{2}$, $D_{3}$, $D_{4}$, $D_{6}$, $T$, $O$). Based on
the symmetry constraints imposed by these point groups in a vector field, we
study the general expression for the radial spin textures up to third order in
momentum. Furthermore, we determine the high-symmetry points of the 45 non-
polar chiral space groups supporting a radial spin texture. These two
principles are used to screen materials databases for archetypes that go
beyond the basic hedgehog radial spin texture. Among the selected materials we
highlight the axion insulator candidate
$\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$, the material proposed for dark
matter detection $\mathrm{Ag}_{3}\mathrm{Au}\mathrm{Te}_{2}$ and heazlewoodite
$\mathrm{Ni}_{3}\mathrm{S}_{2}$, a conventional metal predicted to exhibit
current-induced spin polarization. We point out that the symmetry analysis
proposed in this Letter is more general and extends to studying other vector
properties in momentum space.

Many spin-related phenomena in condensed matter physics are understood in
terms of the specific distribution in momentum space of
$\textbf{S}_{n}(\textbf{k})=\langle\psi_{n}(\textbf{k})|\hat{\textbf{S}}|\psi_{n}(\textbf{k})\rangle$,
the spin expectation value of the Bloch wave function $\psi_{n}(\textbf{k})$
of band $n$ with momentum k. This quantity provides a spin texture (ST) in
momentum space which is determined by the interplay between spin-orbit
coupling (SOC) and crystal symmetries [1, 2]. The archetypal examples of such
relationship are the effective $k\cdot p$ models of the Dresselhaus [3] and
Rashba [4, 5] spin-orbit couplings for non-centrosymmetric zinc blende and
wurtzite crystal structures, respectively. However, since a $k\cdot p$ model
depends on the crystal symmetries as well as the specific manifold of the
Bloch states that are considered, it is a challenging task to obtain effective
models to describe STs in a generic crystal. Such complex dependence on the
details of the material is well illustrated by the existence of exotic Rashba-
type spin textures with cubic momentum-dependent terms in the Bi/Ag(111)
surface alloy [6] and rare-earth ternary materials [7]. Recently, a new method
based on topological quantum chemistry [8] has been proposed to address this
problem. This method correlates the nature of the Bloch functions in a local
basis around a multifold fermion to predict its ST [9].

An alternative approach for obtaining information on the ST, or the momentum
dependence of any vector field, consists in examining its transformation under
the crystal symmetries. New STs have been identified on the basis of such
symmetry analysis. For instance, the persistent ST, a uniform distribution of
spins in momentum space, has been shown to be enforced by symmetry [10, 11].
This type of ST also presents an exotic version with cubic momentum dependence
around high-symmetry points with $C_{3h}$ and $D_{3h}$ point group (PG)
symmetries [12]. As highlighted in Ref. 13, the relevant symmetry operations
that the vector field must satisfy change across the Brillouin zone (BZ) and
are dictated by the PG symmetries of the wavevector around which the ST is
evaluated.

Radial spin distributions also emerge as new types of texture imposed by
crystal symmetry [14]. Such STs appear around high-symmetry points in the BZ
where spins are forced to point parallel to the momentum along several
rotation symmetry axes. The lack of inversion and mirror symmetries, and the
presence of three or more rotational axes are the necessary conditions for the
radial STs. The only PGs that respect these conditions are the non-polar
chiral $D_{2}$, $D_{3}$, $D_{4}$, $D_{6}$, $T$ and $O$.

Radial STs are commonly pictured either as a hedgehog configuration with the
spins pointing parallel to the momentum along all the rotation symmetry axes
or as a vector field where the spins point parallel and antiparallel to the
momentum at different rotation symmetry axes [15, 16, 17]. The first reported
examples of materials hosting such STs are CoSi [18] and elemental Te [14],
two well-known chiral materials. Later, these STs were associated with the
special Kramers-Weyl (KW) points located at time-reversal invariant momentum
(TRIM) points in chiral materials. Experimental confirmation of radial STs via
spin-resolved angle-resolved photoemission spectroscopy (ARPES) has been
provided in Te at both TRIM [19] and non-TRIM high-symmetry points [20], as
well as in PtGa [21]. Theory predicted that materials with such spin
configurations present a new magnetoelectric effect where a magnetization
parallel to the current is induced [22, 23]. This effect has been
experimentally confirmed in Te [24, 25, 26],
${\mathrm{CrNb}}_{3}{\mathrm{S}}_{6}$ [27] and chiral disilicides [28]. The
advance opens a new route for the generation of long-range spin accumulation
that can be used for spin manipulation by electrical means only [29, 26, 30].

In this Letter, we propose an alternative route for classifying spin textures
in momentum space. We systematically analyze the constraints imposed by the
symmetries with chiral non-polar little groups on a general vector field at
high-symmetry points in momentum space. We provide a list of all the high-
symmetry points (both TRIM and non-TRIM) from the 45 non-polar chiral space
groups where it is possible to find a radial ST. Using examples of several
known materials, we illustrate different types of radial STs according to the
PG symmetry, the presence of time-reversal symmetry around the high-symmetry
point, and the need to include higher-order momentum terms.

The ST
$\textbf{S}(\textbf{k})=(S^{x}(\textbf{k}),S^{y}(\textbf{k}),S^{z}(\textbf{k}))$
around $\mathbf{k}_{0}$ (we omit the band index for simplicity), can be
expressed as the Taylor expansion of its components

$S^{\alpha}(\textbf{q})=\sum_{i,j,k}S^{\alpha}_{ijk}q_{x}^{i}q_{y}^{j}q_{z}^{k},\,\,\,\,\,\alpha=x,y,z,$
(1)

where $\textbf{q}=\textbf{k}-\textbf{k}_{0}$ is the momentum measured from the
expansion point $\textbf{k}_{0}$, typically a high-symmetry point. The real
coefficients $S^{\alpha}_{ijk}$ are not arbitrary and must satisfy the
symmetry requirements of the vector field around $\textbf{k}_{0}$, _i.e.,_ the
symmetries of the group of $\mathcal{G}_{\textbf{k}_{0}}$. Since translations
do not affect momentum space properties, only PG symmetries determine how a
vector or pseudo-vector transforms. Let $g$ be the matrix representation of a
symmetry of $\mathcal{G}_{\textbf{k}_{0}}$. A pseudo vector must transform as

$\textbf{S}(\textbf{q})=\mathrm{det}(g)\,g\textbf{S}(g^{-1}\textbf{q}),\,\,\,\,\,g\in\mathcal{G}_{\mathbf{k}_{0}}.$
(2)

In the case of non-polar chiral PGs only proper rotations with
$\mathrm{det}(g)=1$ are present, then pseudo-vectors transform as normal
vectors. By imposing the symmetry constraints of each non-polar chiral PG, we
determine which $S^{\alpha}_{ijk}$ coefficients vanish up to the third order.
Table 1 summarizes the most general expressions of the vector field compatible
with each non-polar chiral PG. These expressions describe the ST up to a
global function $f(q)$ that depends on the modulus of q. Let
$\textbf{S}(\textbf{q})$ be a vector from Table 1 that satisfies Eq. (2), then
$\textbf{S'}(\textbf{q})=f(q)\textbf{S}(\textbf{q})$ satisfies Eq. (2) as
well. Real coefficients $A_{i}$, $B_{i}$ and $C_{i}$ correspond to non-zero
parameters of the linear, quadratic and cubic terms, respectively. .

Table 1:  Expansion terms of the spin texture up to the third order for the
non-polar chiral point groups.

[.5]0 Point group | $\alpha$ | $1^{st}$ | $2^{nd}$ | $3^{rd}$
---|---|---|---|---
[.5]0 $D_{2}$ | $x:$ | $A_{1}q_{x}$ | $B_{1}q_{y}q_{z}$ | $C_{1}q_{x}^{3}+C_{2}q_{x}q_{y}^{2}+C_{3}q_{x}q_{z}^{2}$
[.5]0 | $y:$ | $A_{2}q_{y}$ | $B_{2}q_{x}q_{z}$ | $C_{4}q_{y}^{3}+C_{5}q_{y}q_{z}^{2}+C_{6}q_{y}q_{z}^{2}$
[.5]0 | $z:$ | $A_{3}q_{z}$ | $B_{3}q_{x}q_{y}$ | $C_{7}q_{z}^{3}+C_{8}q_{z}q_{x}^{2}+C_{9}q_{z}q_{y}^{2}$
[.5]0 $D_{4}$ |  | $A_{1}q_{x}$ | $B_{1}q_{y}q_{z}$ | $C_{1}q_{x}^{3}+C_{2}q_{x}q_{y}^{2}+C_{3}q_{x}q_{z}^{2}$
[.5]0 |  | $A_{1}q_{y}$ | -$B_{1}q_{x}q_{z}$ | $C_{1}q_{y}^{3}+C_{2}q_{y}q_{x}^{2}+C_{3}q_{y}q_{z}^{2}$
[.5]0 |  | $A_{2}q_{z}$ |  | $C_{4}q_{z}^{3}+C_{5}q_{z}(q_{x}^{2}+q_{y}^{2})$
[.5]0 $D_{3}$ |  | $A_{1}q_{x}$ | -$2B_{1}q_{x}q_{y}+B_{2}q_{y}q_{z}$ | $C_{1}(q_{x}^{3}+q_{x}q_{y}^{2})+C_{2}q_{x}q_{z}^{2}$
[.5]0 |  | $A_{1}q_{y}$ | $B_{1}(q_{y}^{2}-q_{x}^{2})-B_{2}q_{x}q_{z}$ | $C_{1}(q_{y}^{3}+q_{y}q_{x}^{2})+C_{2}q_{y}q_{z}^{2}$
[.5]0 |  | $A_{2}q_{z}$ |  | $C_{3}(q_{x}^{3}-3q_{x}q_{y}^{2})+C_{4}q_{z}(q_{x}^{2}+q_{z}^{2})+C_{5}q_{z}^{3}$
[.5]0 $D_{6}$ |  | $A_{1}q_{x}$ | $B_{1}q_{y}q_{z}$ | $C_{1}(q_{x}^{3}+q_{x}q_{y}^{2})+C_{2}q_{x}q_{z}^{2}$
[.5]0 |  | $A_{1}q_{y}$ | -$B_{1}q_{x}q_{z}$ | $C_{1}(q_{y}^{3}+q_{y}q_{x}^{2})+C_{2}q_{y}q_{z}^{2}$
[.5]0 |  | $A_{2}q_{z}$ |  | $C_{3}q_{z}(q_{x}^{2}+q_{z}^{2})+C_{4}q_{z}^{3}$
[.5]0 $T$ |  | $A_{1}q_{x}$ | $B_{1}q_{y}q_{z}$ | $C_{1}q_{x}^{3}+C_{2}q_{x}q_{y}^{2}+C_{3}q_{x}q_{z}^{2}$
[.5]0 |  | $A_{1}q_{y}$ | $B_{1}q_{x}q_{z}$ | $C_{1}q_{y}^{3}+C_{2}q_{y}q_{z}^{2}+C_{3}q_{y}q_{x}^{2}$
[.5]0 |  | $A_{1}q_{z}$ | $B_{1}q_{x}q_{y}$ | $C_{1}q_{z}^{3}+C_{2}q_{z}q_{x}^{2}+C_{3}q_{z}q_{y}^{2}$
[.5]0 $O$ |  | $A_{1}q_{x}$ |  | $C_{1}q_{x}^{3}+C_{2}q_{x}(q_{y}^{2}+q_{z}^{2})$
[.5]0 |  | $A_{1}q_{y}$ |  | $C_{1}q_{y}^{3}+C_{2}q_{y}(q_{x}^{2}+q_{z}^{2})$
[.5]0 |  | $A_{1}q_{z}$ |  | $C_{1}q_{z}^{3}+C_{2}q_{z}(q_{x}^{2}+q_{y}^{2})$

Figure 1: Examples of radial spin textures with linear momentum dependence.
Respective materials and band indices are labeled. In all the cases the sphere
surrounds the $\Gamma$ point.

In order to identify all possible points in momentum space that can support
radial STs in generic crystals, we determine the high-symmetry points, for the
45 eligible space groups, with non-polar chiral PGs. Following the double-
value representation tables in Ref. 31, we identify 149 high-symmetry points
with a non-polar chiral PG. These points are listed in Table I of the
Supplemental Material [32]. Only the points not intersected by high-symmetry
lines or planes with symmetry-enforced degeneracies are selected, otherwise
the continuity criteria required for performing Taylor expansion is not
satisfied.

For a given space group with PG $\mathcal{G}$, the high-symmetry points of the
corresponding BZ have a point group
$\mathcal{G}_{\textbf{k}_{0}}\subseteq\mathcal{G}$. When
$\mathcal{G}_{\textbf{k}_{0}}$ is a subset of $\mathcal{G}$, the star of the
high-symmetry point $\textbf{k}_{0}$ contains more than one element. At each
of these non-equivalent points the spin vector field bears a different Taylor
expansion, but the coefficients $S^{\alpha}_{ijk}$ of the vector field
expansion at non-equivalent points of the star are related by the
complementary symmetry operations that form subgroup
$\mathcal{Q}_{\textbf{k}}=\mathcal{G}-\mathcal{G}_{\textbf{k}_{0}}$. These
additional relations among $S^{\alpha}_{ijk}$ at different high-symmetry
points of the star can introduce new constraints on the vector field and even
reduce the number of constants required to describe such field. The details of
constraints imposed by the PG of the space group are given in the Supplemental
Material [32].

Figure 2: Examples of radial spin textures with quadratic momentum dependence.
The corresponding high-symmetry points in the Brillouin zone are indicated in
parentheses. Left panels are the computed spin textures and right panels are
the normalized spin texture for the same band.

Next, we illustrate the diversity of spin textures around high-symmetry points
for each non-polar chiral PG 111We compute, using density functional theory,
spin textures around small spheres of radius 0.005 $2\pi/a$, where $a$ is the
lattice constant of the material, surrounding the high-symmetry points. See
the Supplemental Material [32] for details of calculations. using a selection
of chiral materials listed in Table 2. These materials are chosen according to
the following criteria. Both $\mathrm{CoSi}$ and $\mathrm{TaSi}_{2}$ highlight
the complexity of the radial ST in well-known chiral crystals that were
assumed to host a simple hedgehog spin configuration. We have also included
the axion insulator candidate $\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$ [34]
and the material proposed for dark matter detection
$\mathrm{Ag}_{3}\mathrm{AuTe}_{2}$ [35] to emphasize the possibility of
simultaneous manifestation of several exciting electronic phenomena in these
novel materials. Finally, in order to complete the list of STs generated by
the different non-polar chiral PGs, we screened the TopoMat database [36] with
the help of Table I in the Supplemental Material [32] to find new materials
that can host radial STs. We selected heazlewoodite
$\mathrm{Ni}_{3}\mathrm{S}_{2}$ [37] and the non-magnetic phase of
$\mathrm{CoTeMoO}_{6}$ as examples of the $D_{3}$ and $D_{2}$ PG symmetries.
Special interest must be paid to heazlewoodite, a conventional metal that is
susceptible to exhibit the chirality-induced spin polarization. In all cases,
we choose STs produced by the bands that are the most representatives of all
scenarios, bearing in mind that these are not isolated cases.

Table 2:  Selected materials and their space groups illustrating the spin texture imposed by each non-polar point group. [.5]0 Point group | TRIM | Non-TRIM
---|---|---
[.5]0 $D_{2}$ |  $\mathrm{CoTeMoO}_{6}$ (SG 18) |  $\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$ (SG 97)
[.5]0 $D_{3}$ |  $\mathrm{Ni}_{3}\mathrm{S}_{2}$ (SG 155) |  $\mathrm{TaSi}_{2}$ (SG 180)
[.5]0 $D_{4}$ |  $\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$ (SG 97) | -
[.5]0 $D_{6}$ |  $\mathrm{TaSi}_{2}$ (SG 180) | -
[.5]0 $T$ | CoSi (SG 198) |  $\mathrm{Ag_{3}AuTe}_{2}$ (SG 214)
[.5]0 $O$ |  $\mathrm{Ag}_{3}\mathrm{AuTe}_{2}$ (SG 214) | -

So far, the spin textures of non-polar chiral PGs have been studied only up to
linear order. This simplification gave rise to the notion of a hedgehog vector
field. However, hedgehog or purely radial or hedgehog STs are only present
around high-symmetry points with $T$ and $O$ PGs, while for the rest of PGs
the textures cannot be purely radial. Focusing on the linear terms in Table 1,
purely radial STs of the $T$ and $O$ PGs are described by a vector field with
all three components having the same coefficient $A_{1}$. Examples of such
purely radial STs are illustrated in Figs. 1(a) and 1(b) for CoSi ($T$) and
$\mathrm{Ag_{3}AuTe}_{2}$ ($O$) at the $\Gamma$ point.

As the number of rotation axes of the PGs is reduced, the number of parameters
describing the ST increases, and non-radial components along the directions
other than the rotation axes emerge. For the $D_{3}$, $D_{4}$ and $D_{6}$ PGs,
two different parameters $A_{1}\neq A_{2}$ are required. For the $D_{2}$ PG
each spin component has a different coefficient $A_{1}$, $A_{2}$ and $A_{3}$.
Examples of possible STs that correspond to specific bands at the $\Gamma$
point in $\mathrm{Ni}_{3}\mathrm{S}_{2}$,
$\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$, $\mathrm{TaSi}_{2}$ and
$\mathrm{CoTeMoO}_{6}$ are shown in Figs. 1(c)–1(f), respectively. While Fig.
1(c) shows an almost purely radial ST due to $A_{2}$ being slightly larger
than $A_{1}$, the rest of examples [Figs. 1(d)–1(f)] deviate strongly from the
ideal hedgehog configurations. This is due to $|A_{1}|\approx|A_{2}|$ but
$\mathrm{sign}(A_{1})\neq\mathrm{sign}(A_{2})$, $A_{2}\gg A_{1}$ and
$A_{1}\neq A_{2}\neq A_{3}$, respectively. We would like to emphasize the case
shown in Fig. 1(e), in which $A_{2}$ is significantly larger than $A_{1}$ and
a near-persistent ST is created between the north and south hemisphere.

Figure 3: Examples of radial spin textures with cubic momentum dependence. In
all the cases the sphere surrounds the $\Gamma$ point. Top panels are the
computed spin texture and bottom panels the normalized spin texture for the
same band.

Although radial STs have been associated with KW points [15] and multifold
band degeneracies at TRIM points, such configurations also appear at non-TRIM
points. The seminal work on elemental Te was devoted to the study of the band
structure and the ST around the $H$ point, a non-TRIM point at which the two
highest valence bands are non-degenerate. There is, however, a fundamental
difference since at TRIM points the quadratic terms of Eq. (1) are forced to
be zero due to time-reversal symmetry, while this condition is not present at
non-TRIM points. Non-TRIM high-symmetry points are less common, and we
identified a total of 24 distributed in the $D_{2}$, $D_{3}$ and $T$ point
groups. In order to illustrate the appearance of radial STs due to leading
quadratic momentum terms at non-TRIM points, we selected three materials from
the TRIM column in Table 2 that also have non-TRIM points in their BZ. These
are $\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$, $\mathrm{Ta}\mathrm{Si}_{2}$
and $\mathrm{Ag}_{3}\mathrm{Au}\mathrm{Te}_{2}$ for the $D_{2}$, $D_{3}$ and
$T$ point groups, respectively, and the corresponding STs are shown in Figs.
2(a)–2(c). In order to highlight the details of the STs, we divided each panel
into two parts: left part shows the computed ST, while the right displays the
normalized ST ($\textbf{S}_{n}(\textbf{k})/|\textbf{S}_{n}(\textbf{k})|$) for
the same band. In the following, we describe the contribution of the different
quadratic terms to the three cases. Figure 2(a) shows the ST of
$\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$ for a band around the $P$ point.
The quadratic behaviour of the spin is evident at the equator $q_{z}=0$, where
the $B_{3}q_{x}q_{y}$ term of the $S_{z}$ component is manifested. Similarly,
the selected ST of $\mathrm{Ta}\mathrm{Si}_{2}$ in Fig. 2(b) is governed by
the $B_{1}q_{x}q_{y}$ and $B_{1}(q_{y}^{2}-q_{x}^{2})$ terms of the $S_{x}$
and $S_{y}$ components, respectively. Finally, the case of
$\mathrm{Ag}_{3}\mathrm{Au}\mathrm{Te}_{2}$ at the $P$ point [Fig. 2(c)] shows
a competition between the linear and quadratic terms with
$\mathrm{sign}(A_{1})\neq\mathrm{sign}(B_{1})$.

Furthermore, Fig. 2(a) illustrates an example of ST at a high-symmetry point
$\textbf{k}_{0}$ with a PG $\mathcal{G}_{\textbf{k}_{0}}$, a subgroup of space
group $\mathcal{G}$, being conditioned by the subgroup
$\mathcal{Q}_{\textbf{k}}=\mathcal{G}-\mathcal{G}_{\textbf{k}}$. According to
Table 1, four parameters determine the ST up to quadratic order at the equator
encircling a point with the $D_{2}$ point group. As a consequence, low-
symmetry STs are possible, but in $\mathrm{Ta}_{2}\mathrm{Se}_{8}\mathrm{I}$
the $C_{4}$ rotation symmetry of the crystal point group $\mathcal{G}$ relates
the vector field at non-equivalent points $P$. As a result, at these points
$A_{1}=A_{2}$ and $B_{1}=B_{2}$ resulting in a more symmetric ST.

In general, the complexity of ST around high-symmetry points with non-polar
chiral PGs goes beyond the linear and quadratic terms. For each material in
Table 2 it is possible to find several bands, for which STs require cubic
terms to describe momentum dependence. Figure 3 shows some examples that we
found around the $\Gamma$ point of each material. The top panels show the
computed ST, while the bottom panels illustrate its normalized vectors. In the
example of the $D_{2}$ point group symmetry [Fig. 3(a)], the $S_{z}$ component
is dominated by the $q_{z}q_{x}^{2}$ term since the spin only points almost
parallel to $k_{z}$ as we move along the $k_{x}$-$k_{z}$ plane. Figure 3(b)
illustrates an example for the $D_{3}$ point group with $S_{z}$ being governed
by the $C_{3}(q_{x}^{3}-3q_{x}q_{y}^{2})$ term close to the equator of the
sphere. The ST shown in Fig. 3(c) corresponds to the $D_{4}$ point group. In
this case, the spin texture near the equator is the result of a competition
between the cubic terms $C_{1}$ and $C_{2}$ with
$\mathrm{sign}(C_{1})=\mathrm{sign}(C_{2})$ but $C_{2}>C_{1}$. In Figs. 3(d)
and 3(e), the nature of the cubic spin texture is revealed by a fast change in
the orientation of ST on the sphere. In the case of the $D_{6}$ point group
[Fig. 3(d)], the $S_{z}$ component is governed by the $q_{z}^{3}$ term because
it rapidly decays to zero with $k_{z}$. For the $T$ point group symmetry [Fig.
3(e)], the $q_{x}=0$ and $q_{y}=0$ planes are asymmetric, therefore
$\mathrm{sign}(C_{2})\neq\mathrm{sign}(C_{3})$. Finally, the selected example
for the $O$ point group [Fig. 3(f)] presents a ST with a high winding for each
high-symmetry plane. Such ST is generated by the competition of the $A_{1}$
and $C_{1}$ term with $\mathrm{sign}(A_{1})\neq\mathrm{sign}(C_{1})$.

To summarize, we have studied the radial STs generated by the non-polar chiral
PGs. Based on the symmetry transformation of a pseudo-vector under the
operations of these PGs, we provide the most general expression up to the
third order in momentum expansion. We show that beyond the simplicity that the
name radial spin texture alludes to, many exotic STs with a non-radial
behaviour appears. We also provide a comprehensive list of all the high-
symmetry points where such STs can be found among the 45 non-polar chiral
space groups. These results give a simple description of the radial STs based
on the PG symmetries of the high-symmetry points that can be used for predict
the radial ST and correlate it with the current-induced spin polarization. Our
results allow for straightforward screening of materials databases in order to
find novel compounds with radial STs. Finally, the symmetry analysis provided
in this Letter is not only valid for the spin texture, but for all the other
vector properties in momentum space.

###### Acknowledgements.

D.G.-M. and O.V.Y. acknowledge support by the NCCR MARVEL, a National Centre
of Competence in Research, funded by the Swiss National Science Foundation
(grant number 182892). First-principles calculations were performed at the
Swiss National Supercomputing Centre (CSCS) under project No. s1146. D.G.-M.
thanks the María Zambrano Programme at the University of Alicante founded by
the European union-Next Generation EU as well as J. J. Palacios for insightful
discussions. This study forms part of the Advanced Materials program and was
supported by the Spanish MCIN with funding from European Union
NextGenerationEU and by Generalitat Valenciana through MFA/2022/045.

## References

  * Samokhin [2009] K. Samokhin, Spin–orbit coupling and semiclassical electron dynamics in noncentrosymmetric metals, Ann. Phys. 324, 2385 (2009).
  * Samokhin [2019] K. Samokhin, On the pseudospin description of the electron Bloch bands, Ann. Phys. 407, 179 (2019).
  * Dresselhaus [1955] G. Dresselhaus, Spin-orbit coupling effects in zinc blende structures, Phys. Rev. 100, 580 (1955).
  * Rashba and Sheka [1959] E. I. Rashba and V. I. Sheka, Symmetry of energy bands in crystals of wurtzite type II. Symmetry of bands with spin-orbit interaction, Fiz. Tverd. Tela—Collected Papers II 162, 62–76 (1959).
  * Bihlmayer  _et al._ [2015] G. Bihlmayer, O. Rader, and R. Winkler, Focus on the Rashba effect, New J. Phys. 17, 050202 (2015).
  * Vajna  _et al._ [2012] S. Vajna, E. Simon, A. Szilva, K. Palotas, B. Ujfalussy, and L. Szunyogh, Higher-order contributions to the Rashba-Bychkov effect with application to the Bi/Ag(111) surface alloy, Phys. Rev. B 85, 075404 (2012).
  * Usachov  _et al._ [2020] D. Y. Usachov, I. A. Nechaev, G. Poelchen, M. Güttler, E. E. Krasovskii, S. Schulz, A. Generalov, K. Kliemt, A. Kraiker, C. Krellner, K. Kummer, S. Danzenbächer, C. Laubschat, A. P. Weber, J. Sánchez-Barriga, E. V. Chulkov, A. F. Santander-Syro, T. Imai, K. Miyamoto, T. Okuda, and D. V. Vyalikh, Cubic Rashba effect in the surface spin structure of rare-earth ternary materials, Phys. Rev. Lett. 124, 237202 (2020).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Lin  _et al._ [2022] M. Lin, I. Robredo, N. Schröter, C. Felser, M. G. Vergniory, and B. Bradlyn, Spin-momentum locking from topological quantum chemistry: Applications to multifold fermions, Phys. Rev. B 106, 245101 (2022).
  * Tao and Tsymbal [2018] L. L. Tao and E. Y. Tsymbal, Persistent spin texture enforced by symmetry, Nat. Commun. 9, 2763 (2018).
  * Shi and Li [2020] W. Shi and G. Li, Nodal plane and persistent spin texture in a Weyl semimetal without mirror symmetry, Phys. Rev. B 101, 125207 (2020).
  * Zhao  _et al._ [2020] H. J. Zhao, H. Nakamura, R. Arras, C. Paillard, P. Chen, J. Gosteau, X. Li, Y. Yang, and L. Bellaiche, Purely cubic spin splittings with persistent spin textures, Phys. Rev. Lett. 125, 216405 (2020).
  * Mera Acosta  _et al._ [2021] C. Mera Acosta, L. Yuan, G. M. Dalpian, and A. Zunger, Different shapes of spin textures as a journey through the Brillouin zone, Phys. Rev. B 104, 104408 (2021).
  * Hirayama  _et al._ [2015] M. Hirayama, R. Okugawa, S. Ishibashi, S. Murakami, and T. Miyake, Weyl node and spin texture in trigonal tellurium and selenium, Phys. Rev. Lett. 114, 206401 (2015).
  * Chang  _et al._ [2018] G. Chang, B. J. Wieder, F. Schindler, D. S. Sanchez, I. Belopolski, S.-M. Huang, B. Singh, D. Wu, T.-R. Chang, T. Neupert, _et al._ , Topological quantum properties of chiral crystals, Nat. Mater. 17, 978 (2018).
  * Dil [2020] J. H. Dil, Finding spin hedgehogs in chiral crystals, Physics 13, 45 (2020).
  * Tan  _et al._ [2022] W. Tan, X. Jiang, Y. Li, X. Wu, J. Wang, and B. Huang, A unified understanding of diverse spin textures of kramers–weyl fermions in nonmagnetic chiral crystals, Adv. Funct. Mater. 32, 2208023 (2022).
  * [18] F. Ishii, H. Kotaka, and T. Onishi, Spin–orbit interaction effects in the electronic structure of B20-type CoSi: First-principles density functional study, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2013)_.
  * Gatti  _et al._ [2020] G. Gatti, D. Gosálbez-Martínez, S. S. Tsirkin, M. Fanciulli, M. Puppin, S. Polishchuk, S. Moser, L. Testa, E. Martino, S. Roth, P. Bugnon, L. Moreschini, A. Bostwick, C. Jozwiak, E. Rotenberg, G. Di Santo, L. Petaccia, I. Vobornik, J. Fujii, J. Wong, D. Jariwala, H. A. Atwater, H. M. Rønnow, M. Chergui, O. V. Yazyev, M. Grioni, and A. Crepaldi, Radial spin texture of the Weyl fermions in chiral tellurium, Phys. Rev. Lett. 125, 216402 (2020).
  * Sakano  _et al._ [2020] M. Sakano, M. Hirayama, T. Takahashi, S. Akebi, M. Nakayama, K. Kuroda, K. Taguchi, T. Yoshikawa, K. Miyamoto, T. Okuda, K. Ono, H. Kumigashira, T. Ideue, Y. Iwasa, N. Mitsuishi, K. Ishizaka, S. Shin, T. Miyake, S. Murakami, T. Sasagawa, and T. Kondo, Radial spin texture in elemental tellurium with chiral crystal structure, Phys. Rev. Lett. 124, 136404 (2020).
  * Krieger  _et al._ [2022] J. A. Krieger, S. Stolz, I. Robredo, K. Manna, E. C. McFarlane, M. Date, E. B. Guedes, J. H. Dil, C. Shekhar, H. Borrmann, Q. Yang, M. Lin, V. N. Strocov, M. Caputo, B. Pal, M. D. Watson, T. K. Kim, C. Cacho, F. Mazzola, J. Fujii, I. Vobornik, S. S. P. Parkin, B. Bradlyn, C. Felser, M. G. Vergniory, and N. B. M. Schröter, Parallel spin-momentum locking in a chiral topological semimetal (2022), arXiv:2210.08221 [cond-mat.mes-hall] .
  * Yoda  _et al._ [2015] T. Yoda, T. Yokoyama, and S. Murakami, Current-induced orbital and spin magnetizations in crystals with helical structure, Sci. Rep. 5, 12024 (2015).
  * Yoda  _et al._ [2018] T. Yoda, T. Yokoyama, and S. Murakami, Orbital Edelstein effect as a condensed-matter analog of solenoids, Nano Lett. 18, 916 (2018).
  * Furukawa  _et al._ [2017] T. Furukawa, Y. Shimokawa, K. Kobayashi, and T. Itou, Observation of current-induced bulk magnetization in elemental tellurium, Nat. Commun. 8, 954 (2017).
  * Furukawa  _et al._ [2021] T. Furukawa, Y. Watanabe, N. Ogasawara, K. Kobayashi, and T. Itou, Current-induced magnetization caused by crystal chirality in nonmagnetic elemental tellurium, Phys. Rev. Research 3, 023111 (2021).
  * Calavalle  _et al._ [2022] F. Calavalle, M. Suárez-Rodríguez, B. Martín-García, A. Johansson, D. C. Vaz, H. Yang, I. V. Maznichenko, S. Ostanin, A. Mateo-Alonso, A. Chuvilin, _et al._ , Gate-tuneable and chirality-dependent charge-to-spin conversion in tellurium nanowires, Nat. Mater. 21, 526 (2022).
  * Inui  _et al._ [2020] A. Inui, R. Aoki, Y. Nishiue, K. Shiota, Y. Kousaka, H. Shishido, D. Hirobe, M. Suda, J.-i. Ohe, J.-i. Kishine, H. M. Yamamoto, and Y. Togawa, Chirality-induced spin-polarized state of a chiral crystal ${{\mathrm{CrNb}}_{3}{\mathrm{S}}_{6}}$, Phys. Rev. Lett. 124, 166602 (2020).
  * Shiota  _et al._ [2021] K. Shiota, A. Inui, Y. Hosaka, R. Amano, Y. Ōnuki, M. Hedo, T. Nakama, D. Hirobe, J.-i. Ohe, J.-i. Kishine, H. M. Yamamoto, H. Shishido, and Y. Togawa, Chirality-induced spin polarization over macroscopic distances in chiral disilicide crystals, Phys. Rev. Lett. 127, 126602 (2021).
  * He  _et al._ [2021] W.-Y. He, X. Y. Xu, and K. T. Law, Kramers Weyl semimetals as quantum solenoids and their applications in spin-orbit torque devices, Commun. Phys. 4, 66 (2021).
  * Roy  _et al._ [2022] A. Roy, F. T. Cerasoli, A. Jayaraj, K. Tenzin, M. B. Nardelli, and J. Sławińska, Long-range current-induced spin accumulation in chiral crystals, npj Comput. Mater. 8, 243 (2022).
  * Bradley and Cracknell [2010] C. Bradley and A. Cracknell,  _The mathematical theory of symmetry in solids: representation theory for point groups and space groups_ (Oxford University Press, 2010).
  * [32] See Supplemental Material at http://… for complete tables listing high-symmetry points, details of first-principles calculations, band structure of selected materials and the relation between the vector fields at different non-equivalent points of a given star.
  * Note [1] We compute, using density functional theory, spin textures around small spheres of radius 0.005 $2\pi/a$, where $a$ is the lattice constant of the material, surrounding the high-symmetry points. See the Supplemental Material [32] for details of calculations.
  * Gooth  _et al._ [2019] J. Gooth, B. Bradlyn, S. Honnali, C. Schindler, N. Kumar, J. Noky, Y. Qi, C. Shekhar, Y. Sun, Z. Wang, B. A. Bernevig, and C. Felser, Axionic charge-density wave in the Weyl semimetal (TaSe4)2I, Nature 575, 315 (2019).
  * Sánchez-Martínez  _et al._ [2019] M.-Á. Sánchez-Martínez, I. Robredo, A. Bidaurrazaga, A. Bergara, F. de Juan, A. G. Grushin, and M. G. Vergniory, Spectral and optical properties of $\mathrm{Ag}_{3}\mathrm{Au}(\mathrm{Se}_{2},\mathrm{Te}_{2}$) and dark matter detection, J. Phys. Mater. 3, 014001 (2019).
  * Autès  _et al._ [2019] G. Autès, Q. S. Wu, N. Mounet, and O. V. Yazyev, Topomat: a database of high-throughput first-principles calculations of topological materials (2019), https://www.materialscloud.org/discover/topomat.
  * Lu  _et al._ [1996] Z. W. Lu, B. M. Klein, and D. J. Singh, Electronic structure of heazlewoodite Ni3S2, Phys. Rev. B 54, 13542 (1996).