# Tests of the parametrizations of Fragmentation Functions using data on
inclusive pion and kaon production in unpolarized $pp$ collisions from the
STAR collaboration and at the NICA project.

E. Christova<EMAIL_ADDRESS>Institute for Nuclear Research and Nuclear
Energy, Bulgarian Academy of Sciences, Tzarigradsko chaussée 72, 1784 Sofia,
Bulgaria  D. Kotlorz<EMAIL_ADDRESS>Department of Physics, Opole
University of Technology, 45-758 Opole, Proszkowska 76, Poland  Bogoliubov
Laboratory of Theoretical Physics, JINR, 141980 Dubna, Russia  E. Leader
<EMAIL_ADDRESS>Imperial College London, London SW7 2AZ, United
Kingdom

###### Abstract

The goal of this study is to check which, if any, of the published versions of
the pion and kaon fragmentation functions is compatible with the STAR data on
semi-inclusive pion and kaon production in proton-proton collisions, and on
the basis of this analysis to make reliable predictions for the $p_{T}$
spectra of the pions and kaons in inclusive pion and kaon production at the
future NICA proton-proton collider. The calculations are carried out in next-
to-leading order of quantum chromodynamics, using the well tested CTEQ6 parton
distributions. We consider the following pion and kaon fragmentation functions
– DSEHS-14 [1], DSEHS-17 [2], LSS-15 [3], HKNS-07 [4] and AKK-08 [5]. Our
analysis shows that within the experimental errors all tested sets of
fragmentation functions provide a good fit to STAR data, and the best ones are
LSS-15 for pions and DSEHS-17 for kaons. From comparison of the LO and NLO
results it is clear that the latter fit data much better, specially in the
region of small $p_{T}$. The NLO cross sections are also less scale-$Q^{2}$
dependent, where $p_{T}/2\leqslant Q\leqslant 2p_{T}$, than the LO ones. It
seems that the NICA measurements will allow to discriminate among AKK and
other fragmentation functions for charged pions and among HKNS and DSEHS for
charged kaons. Unfortunately, due to the relatively large uncertainties,
measurements of the difference cross sections, $\sigma^{h^{+}-h^{-}}$, will
not be able to distinguish between different FFs.

##  I Introduction

Study of the spectra of identified hadrons at high $p_{T}$ in $pp$ collisions
is based on model calculations in perturbative quantum chromodynamics (pQCD)
[6]. Both in leading order (LO) [7, 8, 9, 10, 11] and in next-to-leading order
(NLO) [12, 13, 14] of pQCD calculations, inclusive production of single
hadrons is described in terms of parton distribution functions (PDFs), parton-
parton interaction cross sections calculated in the Standard Model (SM), and
fragmentation functions (FFs). The PDFs are fairly well known, while on the
contrary, the flavor-separated quark and gluon FFs, being relatively new
objects, required for a quantitative description of hard scattering processes
involving identified light hadrons in the final-state, are not so well
constrained. At present there are several sets of FFs that fairly well
describe data, but nevertheless differ quite a lot in different kinematic
regions. Most directly the FFs have been extracted from one-hadron production
in electron-positron collisions. However, this process in principle, cannot
distinguish the quark and anti-quark FFs and information only about
$D^{h}_{q+\bar{q}}$ is obtained. In order to obtain separate quark and anti-
quark FFs, one-hadron semi-inclusive $lN$ and $pp$ processes play an essential
role. However, as the hadron structure enters in comparing theoretical
calculations with experimental data, different model assumptions have to be
made. One of the goals in the experimental measurements by the STAR
collaboration of pions and kaons at high $p_{T}$ in $pp$ collisions is to
provide more information on the FFs.

In this paper, based fully on pQCD and taking into account all partonic cross
sections in NLO, we calculate the hadronic $p_{T}$-spectra of the final pions
and kaons in single hadron production in $pp$ collisions in the high $p_{T}$
region, and compare it to the STAR data from 2006 [15], 2007 [16] and 2012
[17]. We use various published parametrizations for the FFs: HKNS-07 [4],
AKK-08 [5], DSEHS-14 [1], DSEHS-17 [2] and LSS-15 [3], aiming to obtain
constraints on the available sets of FFs from the sensitivity of STAR data to
the used FFs. We use CTEQ6 parametrization for PDFs [18]. Using these data to
constrain the FFs, should lead to more reliable predictions for future
measurements. In particular we consider the effect on data expected from the
future NICA accelerator.

Our analysis is partially based on the model independent approach of the
difference cross sections of $h^{+}$ and $h^{-}$-production, developed in
[19]. It suggests that, both in $e^{+}e^{-}$, SIDIS and proton-proton
collisions, if instead of the X-sections $d\sigma_{e^{+}e^{-}}^{h}$,
$d\sigma_{N}^{h}$ and $d\sigma_{pp}^{h}$ for inclusive production of hadrons
and their antiparticles, one deals with their differences, i.e. with
$d\sigma_{e^{+}e^{-}}^{h-\bar{h}}$, $d\sigma_{N}^{h-\bar{h}}$ and
$d\sigma_{pp}^{h-\bar{h}}\equiv d\sigma_{pp}^{h}-d\sigma_{pp}^{\bar{h}}$, one
determines directly the non-singlet (NS) combinations of FFs
$D_{q}^{h-\bar{h}}=D_{q_{V}}^{h}$. This result follows solely from charge-
conjugation invariance of the strong interactions without any model
assumptions about the sea-quarks or favoured and unfavoured fragmentation
functions.

The STAR Collaboration has presented data for pions and kaons. In comparing
with the data we make two simplifications in our formula for the spectra of
the pions and the kaons. 1) We use $s=\bar{s}$ for the strange quark PDFs.
Note that this is not an assumption but follows from the strong limit
$|s-\bar{s}|<0.025$ obtained in neutrino experiments [20] and implies that the
contribution from terms proportional to $s-\bar{s}$ will be within the
experimental errors and thus negligible. 2) We assume SU(2) isospin symmetry
of the strong interactions. The theoretical framework for inclusive production
of single hadrons in $pp$ collisions is presented in Section II, and the
expressions for the difference cross sections and the conditions under which
they have been derived are given in Section III. The comparison of our results
to the STAR data is presented in Section IV. Here, we also compare the results
obtained within NLO and LO approaches and study the scale dependence of the
cross sections for the charged hadrons at LO and NLO where the scale $Q$ may
vary from $p_{T}/2$ to $2p_{T}$.

Unfortunately we cannot take advantage of the very simple expressions derived
in Section III for the cross-section differences, because the experimental
errors on these are too large.

Based on the obtained results for the FFs, we calculate the expected
$p_{T}$-spectra of the pions and kaons produced in $pp$ collisions at the
kinematics of the planned NICA accelerator. Specific to NICA are the
relatively low energies $\sqrt{S}\simeq 27\,GeV$ and the very large solid
angle $0.1\pi<\theta<0.9\pi$, almost the whole forward hemisphere in which the
final hadrons will be detected. This is presented in Section V. We summarize
the obtained results in a Conclusion.

##  II The theoretical framework for $pp\to h+X$

We consider the process:

$\displaystyle p(P_{A})+p(P_{B})\to h(P^{h})+X$ (1)

The basic concept [6] underlying the theoretical analysis of most high energy
interactions is the factorization theorem, which states that the cross-section
for large momentum-transfer reactions may be factorized into long-distance
pieces that contain the desired information on the structure of the nucleon in
terms of its parton densities such as $q(x)$ and fragmentation functions
$D^{h}_{q}(z)$, and short-distance parts which describe the hard interactions
of the partons. The strength of this interaction is controlled by the running
coupling constant evaluated at a large scale associated with the hard
interaction which we denote by $Q^{2}$. The two crucial points here are, on
one hand the long-distance contributions are universal, i.e., they are the
same in any inelastic reaction under consideration, and on the other hand the
short-distance pieces depend only on the large scales related to the large
momentum transfer in the overall reaction and, therefore, can be evaluated
using QCD perturbation theory.

Calculations in LO, Eq. (1), have been considered and compared to data in
several papers in [7, 8, 9, 10, 11]. The first calculations beyond LO, at
$(\alpha_{s}^{3})$ order in perturbative QCD, that include the virtual
corrections in the $2\to 2$ and the $2\to 3$ partonic subprocesses, have been
calculated in [12, 13, 14], and the full $\mathcal{O}(\alpha_{s}^{3})$
radiative NLO corrections are presented in [14].

The results presented for the unpolarized cross section for $pp\to\pi+X$ at
$\sqrt{S}=200$ GeV, that are well described by the NLO QCD calculation [21],
provide confidence that the theoretical framework based on perturbative-QCD
hard scattering is adequate.

Here we consider cross sections $pp\to h^{\pm}+X$ and also difference cross
sections $\sigma(pp\to h^{+}+X)-\sigma(pp\to h^{-}+X)$. The calculations are
done in NLO in the QCD improved parton model following the approach of [14]
which include all parton subprocesses involving quarks and gluons. In our
numerical analysis, we have adapted the INCNLO inclusive hadron production
code [22] and also a code for calculating the fragmentation functions [4].

In the simple parton model process (1) proceeds via the $2\to 2$ partonic
subprocesses:

$\displaystyle q_{a}(p_{a})+q_{b}(p_{b})\to c(p_{c})+d(p_{d})$ (2)

where $q_{a},q_{b}$ and $c,d$ can be either quarks or gluons. In what follows
all kinematic variables refer to the CM of the $pp$ collision. The expression
for the cross section for $pp\to hX$ in the c.m.s. of $pp$ has the factorized
form:

$\displaystyle E^{h}\frac{d\sigma_{pp}^{h}}{d^{3}P^{h}}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\,\sum_{ab\to
cd}\int^{1}_{x_{a,min}}dx_{a}\int^{1}_{x_{b,min}}dx_{b}\,\frac{1}{z}\,\times$
(3)
$\displaystyle\times\left\\{q_{a}(x_{a})q_{b}(x_{b})\,\left[\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{c}^{h}(z)+\frac{d\hat{\sigma}_{ab}^{cd}}{du}\,D_{d}^{h}(z)\right]+q_{a}(x_{b})q_{b}(x_{a})\,\left[\frac{d\hat{\sigma}_{ab}^{cd}}{du}\,D_{c}^{h}(z)+\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{d}^{h}(z)\right]\right\\}.$

where $d\hat{\sigma}_{ab}^{cd}$ are the Born cross sections to order
$\alpha^{2}_{s}$ and the parton densities and FFs are scale independent,
$z=E^{h}/E_{c}$, and where

$\displaystyle\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\equiv\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,(s,t,u),\qquad\frac{d\hat{\sigma}_{ab}^{cd}}{du}\equiv\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,(s,u,t).$
(4)

As usual $x_{a}$ ($x_{b}$) is the fraction of the proton momentum $P_{A}$
($P_{B}$) carried by the parton $q_{a}$ ($q_{b}$), collinear to the momentum
of the initial hadron $A$ ($B$).

In the QCD improved parton model, when QCD corrections are included, the
factorized form of (3) is preserved and the PDFs and FFs are replaced by
$Q^{2}$-dependent distribution functions:

$\displaystyle q(x)\to q(x,Q^{2}),\qquad D_{q}^{h}(z)\to D_{q}^{h}(z,Q^{2}),$
(5)

where $Q^{2}$ is some relevant large momentum scale. Furthermore the effective
quark-gluon-quark coupling constant $\alpha_{s}$ becomes the running coupling
given by

$\displaystyle\alpha_{s}(Q^{2})=\frac{4\pi}{\beta_{0}\ln(Q^{2}/\Lambda^{2})}\left(1-\frac{\beta_{1}\ln\ln(Q^{2}/\Lambda^{2})}{\beta_{0}^{2}\ln(Q^{2}/\Lambda^{2})}\right),$
(6)

where $\beta_{0}=11-2/3\,N_{f}$ for $N_{f}$ flavours,
$\beta_{1}=102-38/3\,N_{f}$ and $\Lambda_{\rm NLO}=0.32\,{\rm GeV}$. The
momentum scale $Q^{2}$ is set to the “natural” scale $Q^{2}=p_{T}$ [13]. In
order to use the PDFs and FFs, known from some other process at some other
“input” scale $Q_{0}^{2}$, we evolve the PDFs and FFs from $Q_{0}^{2}$ to
$Q^{2}$ via the NLO DGLAP-equations. In our approach when using the difference
cross sections, only NS combinations of FFs enter – see for example Eqs. (20)
and (29). In this case the DGLAP evolution equations for a given FF is
especially simple with no mixing with other FFs.

The Mandelstam variables of the partonic process (2) are:

$\displaystyle s$ $\displaystyle=$
$\displaystyle(p_{a}+p_{b})^{2}=(x_{a}P_{A}+x_{b}P_{B})^{2}=x_{a}\,x_{b}\,S$
$\displaystyle t$ $\displaystyle=$
$\displaystyle(p_{a}-p_{c})^{2}=(x_{a}P_{A}-\frac{P^{h}}{z})^{2}=\frac{x_{a}}{z}\,T,$
$\displaystyle u$ $\displaystyle=$
$\displaystyle(p_{b}-p_{c})^{2}=(x_{b}P_{B}-\frac{P^{h}}{z})^{2}=\frac{x_{b}}{z}\,U$
(7)

Here $c$ or $d$ are the fragmenting partons that are assumed collinear to the
final hadron $h$, and everywhere masses are neglected, so that $s+t+u=0$
holds, which determines the value of $z$:

$\displaystyle z=-\,\frac{x_{a}\,T+x_{b}\,U}{x_{a}\,x_{b}\,S}$ (8)

The letters $S,\,T$ and $U$ stand for the Mandelstam variables of the
inclusive hadronic process (1):

$\displaystyle S$ $\displaystyle=$
$\displaystyle(P_{A}+P_{B})^{2}=2\,(P_{A}\cdot P_{B})=4E^{2}$ (9)
$\displaystyle T$ $\displaystyle=$
$\displaystyle(P_{A}-P^{h})^{2}=-\,2\,(P_{A}\cdot
P^{h})=-2EE^{h}(1-\cos\theta)$ (10) $\displaystyle U$ $\displaystyle=$
$\displaystyle(P_{B}-P^{h})^{2}=-\,2\,(P_{B}\cdot
P^{h})=-2EE^{h}(1+\cos\theta)$ (11) $\displaystyle p_{T}^{2}$ $\displaystyle=$
$\displaystyle\frac{UT}{S}=(E^{h})^{2}\sin^{2}\theta\,,$ (12)

where, as mentioned, the kinematic variables refer to the CM of the $pp$
collision; $E$ is the energy of the colliding proton beams. The maximal value
of the transverse momentum of a hadron $h$ produced with energy $E^{h}$ in
$pp$ collisions is

$\displaystyle p_{T}^{max}$ $\displaystyle=$ $\displaystyle
E^{h},\,\,\textrm{at}\,\,\theta=\pi/2\,.$ (13)

Often the $p_{T}$-spectra are presented in terms of its relative value
$x_{T}=2p_{T}/\sqrt{S}$, where $\sqrt{S}/2=E$.

The lower limits of integration are determined by the conditions $z<1$ and
$s+t+u=0$ [9]:

$\displaystyle x_{a,min}$ $\displaystyle=$ $\displaystyle\frac{-U}{T+S}$ (14)
$\displaystyle x_{b,min}$ $\displaystyle=$
$\displaystyle\frac{-x_{a}T}{x_{a}S+U}$ (15)

There are 8 different $2\to 2$ partonic processes that contribute to
$d\sigma^{h}_{pp}$:

$\displaystyle\hat{\sigma}_{1}:\,$ $\displaystyle q_{i}q_{j}\to
q_{i}q_{j},\quad\bar{q}_{i}\bar{q}_{j}\to\bar{q}_{i}\bar{q}_{j},\quad
q_{i}\bar{q}_{j}\to q_{i}\bar{q}_{j},i\neq j$
$\displaystyle\hat{\sigma}_{2}:\,$ $\displaystyle q_{i}q_{i}\to
q_{i}q_{i},\quad\bar{q}_{i}\bar{q}_{i}\to\bar{q}_{i}\bar{q}_{i},$
$\displaystyle\hat{\sigma}_{3}:\,$ $\displaystyle q_{i}\bar{q}_{i}\to
q_{j}\bar{q}_{j},\quad i\neq j$ $\displaystyle\hat{\sigma}_{4}:\,$
$\displaystyle q_{i}\bar{q}_{i}\to q_{i}\bar{q}_{i}$
$\displaystyle\hat{\sigma}_{5}:\,$ $\displaystyle q_{i}\bar{q}_{i}\to gg$
$\displaystyle\hat{\sigma}_{6}:\,$ $\displaystyle gg\to q_{i}\bar{q}_{i}$
$\displaystyle\hat{\sigma}_{7}:\,$ $\displaystyle q_{i}g\to q_{i}g$
$\displaystyle\hat{\sigma}_{8}:\,$ $\displaystyle gg\to gg$ (16)

expressions for which can be found in many places. In NLO
$\mathcal{O}(\alpha_{s}^{3})$, additional $2\to 3$ scattering processes also
contribute [12, 13, 14].

The contributions of the partonic processes which occur in Eq. (1) are given
in Appendix A.

##  III The difference cross sections : $h^{+}-h^{-}$

Here we would like to introduce the idea of the difference cross sections and
discuss its possible advantages. We define the difference cross section:

$\displaystyle\sigma_{pp}^{h^{+}-h^{-}}\equiv\sigma_{pp}^{h^{+}}-\sigma_{pp}^{h^{-}}.$
(17)

C-invariance of strong interactions implies:

$\displaystyle D_{g}^{h^{+}-h^{-}}=0,$ (18) $\displaystyle
D_{q}^{h^{+}-h^{-}}=-D_{\bar{q}}^{h^{+}-h^{-}}=D_{q_{V}}^{h^{+}}.$ (19)

The cross section $\sigma_{pp}^{h^{+}-h^{-}}$ has especially simply form in
LO, where, without any assumptions about the sea quark PDFs or the favoured
and unfavoured FFs, we obtain:

$\displaystyle
E^{h}\frac{d\,\sigma_{pp}^{h^{+}-h^{-}}}{d^{3}P^{h}}=\frac{1}{\pi}\int
dx_{a}\,dx_{b}\,\frac{1}{z}\sum_{q=u,d,s}D_{q_{V}}^{h^{+}}(z)\big{(}q_{V}(x_{a})\,L_{q}(x_{b},t,u)+q_{V}(x_{b})L_{q}(x_{a},u,t)\big{)},$
(20)

where

$\displaystyle L_{u}(x_{b},t,u)$ $\displaystyle=$
$\displaystyle(\widetilde{d}+\widetilde{s})(x_{b})\,\frac{d\hat{\sigma}_{1}}{dt}+\frac{1}{2}\,\widetilde{u}(x_{b})\left[\frac{d\hat{\sigma}_{2}}{dt}+\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right]+g(x_{b})\,\frac{d\hat{\sigma}_{7}}{dt}$
(21) $\displaystyle L_{d}(x_{b},t,u)$ $\displaystyle=$
$\displaystyle(\widetilde{u}+\widetilde{s})(x_{b})\,\frac{d\hat{\sigma}_{1}}{dt}+\frac{1}{2}\,\widetilde{d}(x_{b})\left[\frac{d\hat{\sigma}_{2}}{dt}+\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right]+g(x_{b})\,\frac{d\hat{\sigma}_{7}}{dt}$
(22) $\displaystyle L_{s}(x_{b},t,u)$ $\displaystyle=$
$\displaystyle(\widetilde{u}+\widetilde{d})(x_{b})\,\frac{d\hat{\sigma}_{1}}{dt}+\frac{1}{2}\,\widetilde{s}(x_{b})\left[\frac{d\hat{\sigma}_{2}}{dt}+\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right]+g(x_{b})\,\frac{d\hat{\sigma}_{7}}{dt}$
(23)

and for example

$\displaystyle
L_{u}(x_{a},u,t)=(\widetilde{d}+\widetilde{s})(x_{a})\,\frac{d\hat{\sigma}_{1}}{du}+\frac{1}{2}\,\widetilde{u}(x_{a})\left[\frac{d\hat{\sigma}_{2}}{dt}+\frac{d\hat{\sigma}_{4}}{du}-\frac{d\hat{\sigma}_{4}}{dt}\right]+g(x_{a})\,\frac{d\hat{\sigma}_{7}}{du}.$
(24)

In the last formula we have used the fact that $\hat{\sigma}_{2}$ is symmetric
under $t\leftrightarrow u$. Here $u_{V}$ and $d_{V}$ are the usual valence
quark PDFs:

$\displaystyle u_{V}=u-\bar{u},\quad d_{V}=d-\bar{d}$ (25)

and we have used the notation:

$\displaystyle s_{V}=s-\bar{s},\qquad\widetilde{q}=q+\bar{q}\,.$ (26)

Using the strong bound on $(s-\bar{s})$ obtained from neutrino experiments
$s_{V}\leq 0.025$ [20] one may safely neglect the contribution of
$s_{V}D_{s_{V}}^{h}$ and we have:

$\displaystyle
E^{h}\frac{d\,\sigma_{pp}^{h^{+}-h^{-}}}{d^{3}P^{h}}=\frac{1}{\pi}\int
dx_{a}\,dx_{b}\,\frac{1}{z}\sum_{q=u,d}D_{q_{V}}^{h^{+}}(z)\big{(}q_{V}(x_{a})\,L_{q}(x_{b},t,u)+q_{V}(x_{b})L_{q}(x_{a},u,t)\big{)}$
(27)

This implies that the contribution of $(s-\bar{s})D_{s}^{h^{+}-h^{-}}$ is
expected to be within the experimental error and negligible, also the large
uncertainties in $D_{s}^{h}$ should not affect the results for
$\sigma_{pp}^{h^{+}-h^{-}}$.

Note that in LO only 4 partonic cross sections contribute to the difference
$d\sigma_{pp}^{h^{+}-h^{-}}$:

– scattering of flavour-unlike quarks $\hat{\sigma}_{1}$,

– scattering of flavor-like quarks $\hat{\sigma}_{2}$ and $\hat{\sigma}_{4}$
and

– quark-gluon scattering $\hat{\sigma}_{7}$.

In NLO, the expressions for the difference cross sections are much more
complicated as all radiative corrections $\mathcal{O}(\alpha_{s}^{3})$ must be
accounted for. Nevertheless, again, one obtains only nonsinglet contributions
proportional to $D_{q_{V}}^{h^{+}}(z,Q^{2})$.

###  III.1 The difference cross section for $\pi^{+}-\pi^{-}$

Eq. (20) considerably simplifies for pions if one assumed SU(2) isospin
symmetry, as used in all present analyses. It implies:

$\displaystyle D_{u}^{\pi^{+}-\pi^{-}}=-D_{d}^{\pi^{+}-\pi^{-}},\qquad
D_{s}^{\pi^{+}-\pi^{-}}=0\cdot$ (28)

Then $\sigma_{pp}^{\pi^{+}-\pi^{-}}$ is expressed solely in terms of
$D_{u_{V}}^{\pi^{+}}$, enhanced by the best known valence-quark $u_{V}$:

$\displaystyle E^{\pi}\frac{d\,\sigma_{pp}^{\pi^{+}-\pi^{-}}}{d^{3}P^{\pi}}$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\int
dx_{a}\,dx_{b}\,\frac{1}{z}\,D_{u_{V}}^{\pi^{+}}(z)\times$ (29)
$\displaystyle\times\left[u_{V}(x_{a})\,L_{u}(x_{b},t,u)-d_{V}(x_{a})\,L_{d}(x_{b},t,u)+u_{V}(x_{b})L_{u}(x_{a},u,t)-d_{V}(x_{b})L_{d}(x_{a},u,t)\right]$

###  III.2 The difference cross section for $K^{+}-K^{-}$

The formula for $\sigma_{pp}^{K^{+}-K^{-}}$ strongly simplifies if we assume:

$\displaystyle D_{d}^{K^{+}-K^{-}}=0$ (30)

This seems a reasonable physical assumption that follows from the quark
content of $K^{\pm}=(\bar{s},u)$, $(s,\bar{u})$ and is used in all current
analyses. Then $\sigma_{pp}^{K^{+}-K^{-}}$ depends solely on one non-singlet
combination of FFs $D_{u_{V}}^{K^{+}}$, multiplied by the large valence
$u_{V}$-quark distributions:

$\displaystyle E^{K}\frac{d\,\sigma_{pp}^{K^{+}-K^{-}}}{d^{3}P^{K}}$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\int
dx_{a}\,dx_{b}\,\frac{1}{z}\,D_{u_{V}}^{K^{+}}(z)\big{[}u_{V}(x_{a})\,L_{u}(x_{b},t,u)+u_{V}(x_{b})L_{u}(x_{a},u,t)\big{]}$
(31)

This expression may be used as a test for the assumption (30).

###  III.3 The difference cross section for $K^{\pm}$ and $K_{s}^{0}$

If in addition to the charged $K^{\pm}$ also neutral kaons
$K_{s}^{0}=(K^{0}+\bar{K}^{0})/\sqrt{2}$ are measured, no new FFs are
introduced into the cross-sections. This is a consequence of SU(2) isospin
symmetry, according to which $(K^{+},K^{0})$ and $(K^{-},\bar{K}^{0})$ form
isospin doublets, and we have:

$\displaystyle
D_{u}^{K^{+}+K^{-}}(z,Q^{2})=D_{d}^{K^{0}+\bar{K}^{0}}(z,Q^{2})=2D_{d}^{K_{s}^{0}}(z,Q^{2})$
$\displaystyle
D_{d}^{K^{+}+K^{-}}(z,Q^{2})=D_{u}^{K^{0}+\bar{K}^{0}}(z,Q^{2})=2D_{u}^{K_{s}^{0}}(z,Q^{2})$
$\displaystyle
D_{s}^{K^{+}+K^{-}}(z,Q^{2})=D_{s}^{K^{0}+\bar{K}^{0}}(z,Q^{2})=2D_{s}^{K_{s}^{0}}(z,Q^{2}).$
(32)

Note that the first relation involves favoured and unfavoured FFs, while the
second – only unfavoured FFs. As explained in [5], due to the nature of the
DGLAP evolution, these constraints are independent of $Q^{2}$, as constraints
that follow from symmetry should be.

We write Eqs. (32) in the form:

$\displaystyle
D_{u}^{K^{+}+K^{-}-2K_{s}^{0}}=-D_{d}^{K^{+}+K^{-}-2K_{s}^{0}}={(D_{u}-D_{d})}^{K^{+}+K^{-}}$
(33) $\displaystyle
D_{s}^{K^{+}+K^{-}-2K_{s}^{0}}=D_{c}^{K^{+}+K^{-}-2K_{s}^{0}}=D_{b}^{K^{+}+K^{-}-2K_{s}^{0}}=D_{g}^{K^{+}+K^{-}-2K_{s}^{0}}=0.$
(34)

Then, as shown in [19], the combination

$\displaystyle\sigma^{K^{+}}+\sigma^{K^{-}}-2\sigma^{K_{s}^{0}}$ (35)

in the three types of inclusive processes, $K=K^{\pm},K_{s}^{0}$:

$\displaystyle e^{+}+e^{-}\to K+X,$ (36) $\displaystyle e+N\to e+K+X,\qquad
N=p,d,$ (37) $\displaystyle p+p\to K+X,$ (38)

measure the same NS combination of FFs, namely $(D_{u}-D_{d})^{K^{+}+K^{-}}$.
This result relies only on SU(2) invariance for the kaons and does not involve
any assumptions about PDFs or FFs; it holds in any order in QCD.

For inclusive hadron production, from Eqs. (3) and (33) we obtain:

$\displaystyle E^{K}\frac{d\sigma_{pp}^{K^{+}+K^{-}-2K_{S}^{0}}}{d^{3}P^{K}}$
$\displaystyle=\frac{1}{\pi}\int\,dx_{a}\int\,dx_{b}\,\frac{1}{z}\,\times$
(39)
$\displaystyle\times\Big{\\{}\left[\tilde{u}(x_{a})[\tilde{d}(x_{b})+\tilde{s}(x_{b})]-\tilde{d}(x_{a})[\tilde{u}(x_{b})+\tilde{s}(x_{b})]\right]\frac{d\hat{\sigma}_{1}}{dt}+$
$\displaystyle+\frac{1}{2}\,\left[u(x_{a})u(x_{b})+\bar{u}(x_{a})\bar{u}(x_{b})-[d(x_{a})d(x_{b})+\bar{d}(x_{a})\bar{d}(x_{b})]\right]\frac{d\hat{\sigma}_{2}}{dt}+$
$\displaystyle+\left[d(x_{a})\bar{d}(x_{b})-u(x_{a})\bar{u}(x_{b})\right]\left[\frac{d\hat{\sigma}_{3}}{dt}-\frac{d\hat{\sigma}_{4}}{dt}\right]+[\tilde{u}(x_{a})-\tilde{d}(x_{a})]g(x_{b})\,\frac{d\hat{\sigma}_{7}}{dt}$
$\displaystyle+\left[(x_{a}\leftrightarrow x_{b}),(t\leftrightarrow
u)\right]\Big{\\}}\,D_{u-d}^{K^{+}+K^{-}}(z).$

Note that $D_{q}^{h+\bar{h}}$ can be determined directly from $e^{+}e^{-}\to
h+X$ and thus without assumptions, which implies that $D_{q}^{h+\bar{h}}$
should be roughly the same for all sets of FFs. If this appears the case, then
use of Eq. (39) to compare the STAR data on $K^{+}+K^{-}-2K_{s}^{0}$ can be
used as a test for the consistency of the data with SU(2) invariance.

##  IV Comparison to STAR data

In this section we study numerically the impact of the different sets of FFs
on the $p_{T}$-spectra of the final hadrons in inclusive $pp$-collisions and
compare it to STAR data.

Based on the theoretical NLO pQCD-formulas, we examine the sensitivity of the
$p_{T}$-spectra to the FFs. We use STAR data from 2006, 2007 and 2012 to try
to distinguish among the different FFs. First, we compare the cross sections
for $d\sigma_{pp}^{h^{\pm}}$ and $d\sigma_{pp}^{h^{+}-h^{-}}$, calculated with
different sets for $D_{q}^{h}$, to the measured hadron yield. Second, we try
to obtain a fit of the FFs and compare it to the available parametrizations.

Note that there is a principle difference between the commonly used formula
for $\pi^{\pm}$-production and the one for $\pi^{+}-\pi^{-}$. In $\pi^{\pm}$
production FFs with all flavours
$D_{i}^{h},\,\,i=g,u,d,s,\bar{u},\bar{d},\bar{s}$ enter the cross section,
while in the formula for $\pi^{+}-\pi^{-}$ it is only the valence-quark FFs
that enter – the non-singlet combinations $D_{u_{V}}^{h}$ and $D_{d_{V}}^{h}$.
Note also that, being non-singlets, in their $Q^{2}$-evolution they do not mix
with other FFs and this property holds in all orders in QCD.

Finally we use the preferred sets of FFs obtained in the above studies to make
predictions for the $p_{T}$ spectra of the charged pions and kaons produced in
inclusive $pp$ collisions at the planned NICA accelerator.

The results are presented separately for final pions and kaons.

We consider the following sets of NLO in pQCD FFs for pions: HKNS-07 [4],
AKK-08 [5], DSEHS-14 [1], LSS-15 [3] and for kaons: HKNS-07 [4], AKK-08 [5],
DSEHS-17 [2]. These sets of FFs differ by the different processes, whose data
they use, and by the used assumptions. Only AKK and DSEHS use the inclusive
$pp$-process as a source of information about the FFs.

We apply NLO approach [14] to STAR data from 2006 [15], 2007 [16] and 2012
[17]. The kinematics of the used measurements of STAR Collaborations in terms
of $\sqrt{S}$, $p_{T}$ and the rapidity $y$ are:
STAR-2006, $\pi^{\pm}$:

$\displaystyle\qquad\qquad\qquad\sqrt{S}=200\,{\rm GeV},\qquad
0.3<p_{T}<10\,{\rm GeV}/c,\qquad\qquad|y|<0.5$ (40)

STAR-2007

$\displaystyle K^{\pm}:\qquad\qquad\sqrt{S}=200\,{\rm GeV},\qquad
0.25<p_{T}<2.2\,{\rm GeV}/c,\qquad\qquad|y|<0.5$ (41) $\displaystyle
K_{s}^{0}:\qquad\qquad\sqrt{S}=200\,{\rm GeV},\qquad 0.26<p_{T}<4.7\,{\rm
GeV}/c,\qquad\qquad|y|<0.5$ (42)

STAR-2012, $\pi^{\pm},K^{\pm},K_{s}^{0}$:

$\displaystyle\qquad\qquad\qquad\sqrt{S}=200\,{\rm GeV},\qquad
3.0<p_{T}<15\,{\rm GeV}/c,\qquad\qquad|y|<0.5$ (43)

The STAR Collaborations present their results for the hadronic yield
$d^{2}N/dp_{T}\,dy$, extracted from the measured non-single diffractive (NSD)
$p+p$ cross section, while our formulae are for the differential cross
sections $d\sigma^{h}/d^{3}P_{h}$. The relation between them reads:

$\displaystyle\sigma_{NSD}\,\frac{d^{2}N}{(2\pi)\,p_{T}\,dp_{T}\,dy}=\,E_{h}\,\frac{d\sigma^{h}_{pp}}{d^{3}P_{h}}(p_{T},\theta(y))$
(44)

where $\sigma_{NSD}$ is the measured total NSD $p+p$ cross section
($\sigma_{NSD}=30.0\pm 3.5$ mb) and $E_{h}\,d\sigma^{h}_{pp}/(d^{3}P_{h})$ is
our formula for the corresponding cross section – Eqs. (3), (29), (31) or
(39).

In high energy hadron collider physics for the particle with the negligible
mass the rapidity $y$ coincides with the pseudo-rapidity $\eta$:

$\displaystyle
y\approx\eta\equiv\frac{1}{2}\ln\frac{1+\cos\theta}{1-\cos\theta}=-\ln\tan\frac{\theta}{2}\,,$
(45)

and hence

$\displaystyle\theta=2\arctan(e^{-\eta})\approx 2\arctan(e^{-y})\,.$ (46)

As data from the STAR collaboration is presented averaged over rapidity:
$|y|<0.5$, which corresponds to $\theta\approx[\pi/3;\pi/2]$, we compare data
to the following averaged over $y$ expressions:

$\displaystyle\left.\sigma_{NSD}\,\frac{d^{2}N}{(2\pi)\,p_{T}\,dp_{T}\,dy}\right|_{|y|<0.5}=\,\int_{-0.5}^{0.5}dy\,E_{h}\,\frac{d\sigma^{h}_{pp}}{d^{3}P_{h}}(p_{T},\theta(y))\,.$
(47)

###  IV.1 Results for charged pions

We calculate the $p_{T}$ spectra of $\pi^{\pm}$ using the parametrizations for
$D_{i}^{\pi^{\pm}}$ obtained in HKNS-07, AKK-08, DSEHS-14 and LSS-15, and
compare it to data on $\pi^{\pm}$-spectra measured by the STAR Collaborations
STAR-2006 and STAR-2012 at $\sqrt{S}=200$ GeV for different ranges of $p_{T}$
and rapidity $|y|<0.5$ – see Eqs. (40) and (43). All parametrizations are
obtained within the same theoretical framework – NLO in QCD, assuming that
isospin SU(2) symmetry of the $u$ an $d$ quarks for the favored and unfavored
fragmentation functions holds at the initial scale $Q^{2}_{0}$:

$\displaystyle D_{u}^{\pi^{+}}=D_{\bar{d}}^{\pi^{+}},\qquad
D_{\bar{u}}^{\pi^{+}}=D_{d}^{\pi^{+}}=D_{s}^{\pi^{+}}=D_{\bar{s}}^{\pi^{+}}\cdot$
(48)

Using charge conjugation invariance of strong interactions we have, in total,
3 independent FFs:

$\displaystyle D_{u}^{\pi^{+}},\quad D_{d}^{\pi^{+}},\quad
D_{g}^{\pi^{+}}\cdot$ (49)

The obtained FFs differ by the used data and the used parametrizations at
$Q^{2}_{0}$. Bellow we give the main sources of data of the used here FFs.

$\bullet$ In HKNS-07 the parametrizations and their uncertainties are obtained
both in LO and NLO from analysis of the data on $e^{+}e^{-}\to\pi^{\pm}+X$. We
use the NLO parameters and show the error band based on the Hessian method
[23] and the error matrix provided by HKNS [4]. The analysis for kaons
proceeds in the same manner.

$\bullet$ In AKK-08 the parametrizations are obtained from analysis of the
data on $e^{+}e^{-}\to\pi^{\pm}+X$ and $pp\to\pi^{\pm}+X$.

$\bullet$ In DSEHS-14 the parametrizations are obtained in a global analysis
of the data on $e^{+}e^{-}\to\pi^{\pm}+X$, $l+N\to\pi^{\pm}+X$ and
$pp\to\pi^{\pm}+X$ including the STAR-2012 data. The uncertainties estimated
within the Hessian approach are presented as shaded areas in their figures and
the error matrix can be provided upon request. In turn, in our analysis, both
for pions and kaons, we estimate the scale uncertainty when $Q$ varies in the
range $p_{T}/2\leqslant Q\leqslant 2p_{T}$.

$\bullet$ In LSS-15 the parametrizations are obtained using the published
HERMES data from 2013 on pion multiplicities from SIDIS $l+N\to\pi^{\pm}+X$ on
proton and deutron targets.

The results of our comparison to pion production are shown in Figs. 1–5 and in
Table 1.

Figure 1: The cross sections $\sigma_{pp}^{\pi^{+}}$ (left) and
$\sigma_{pp}^{\pi^{-}}$ (right) calculated for different sets of FFs, compared
to STAR-2006 and STAR-2012 data. The one-$\sigma$ uncertainty band for HKNS
FFs estimated by the Hessian method is shown.

Figure 2: The best fit LSS together with $\sigma_{pp}^{\pi^{+}}$ (left) and
$\sigma_{pp}^{\pi^{-}}$ (right) STAR data. The theoretical uncertainties bands
indicate the scale $Q$ variation in the range $p_{T}/2\leqslant Q\leqslant
2p_{T}$.

Figure 3: The ratio (data-theory)/theory for the cross sections:
$\sigma_{pp}^{\pi^{+}}$ (left) and $\sigma_{pp}^{\pi^{-}}$ (right) for
different FFs. The experimental error bars and theoretical uncertainties for
HKNS (dashed) are shown.

In Figs. 1–3, we show the agreement of the $p_{T}$-spectra for the charged
pions $\sigma_{pp}^{\pi^{\pm}}$ calculated for the different sets of FFs with
the STAR-2006 and STAR-2012 data.

Table 1: Values of $\chi^{2}$, Eq. (50), calculated for the single and difference cross sections in the pion production obtained for different FFs parametrizations. FFs | $~{}~{}~{}\pi^{+}~{}~{}~{}$ | $~{}~{}~{}\pi^{-}~{}~{}~{}$ | $~{}~{}\pi^{+}-\pi^{-}~{}~{}$
---|---|---|---
HKNS | 10.2 | 11.7 | 0.136
AKK | 77.2 | 83.3 | -
DSEHS | 4.55 | 4.39 | 0.148
LSS | 2.83 | 2.46 | 0.143

Using the values of $\chi^{2}$ defined as

$\displaystyle\chi^{2}=\frac{1}{N}\sum_{i}^{N}\frac{(\sigma_{i}^{exp}-\sigma_{i}^{th})^{2}}{\Delta_{i}^{2}},$
(50)

where $\Delta_{i}^{2}$ incorporates both statistical and systematic
experimental errors, we are able to determine the best fit to the data. In our
analysis for the pions, we took into account the combined STAR data for
$p_{T}>1\,{\rm GeV}^{2}$ with $N=30$ experimental points (STAR 2006 and STAR
2012). We find that the NLO LSS results agree the best with the data and the
DSEHS fit which is partially based on the STAR-2012 data is comparably good as
LSS. In Table 1, we collect the values of $\chi^{2}$, Eq. (50), for the single
charged pions $\pi^{\pm}$ and also for the difference $\pi^{+}-\pi^{-}$
obtained with the use of HKNS, AKK, DSEHS and LSS sets of FFs.

Figure 4: Left: The difference cross section $\sigma_{pp}^{\pi^{+}-\pi^{-}}$
calculated for different sets of FFs, compared to the STAR data. The
one-$\sigma$ uncertainty band for HKNS FFs is shown. Right: The ratio (data-
theory)/theory for the difference cross section. The experimental error bars
are shown. For better visibility points for DSEHS (LSS) are displaced by -0.1
(+0.1) in $p_{T}$.

Fig. 4 shows a comparison of the difference cross section
$\sigma_{pp}^{\pi^{+}-\pi^{-}}$ calculated for different sets of FFs with the
STAR-2006 and STAR-2012 data. Since the experimental errors are relatively
large, all FFs parametrizations are comparably good for this nonsinglet case.
Finally, in Fig. 5, using LSS FFs, we compare NLO and LO results for
$\sigma_{pp}^{\pi+}$. The right panel shows the $K$ factor, defined as

$\displaystyle K=\frac{\left[d\sigma_{pp}^{h}\right]^{\rm
NLO}}{\left[d\sigma_{pp}^{h}\right]^{\rm LO}},$ (51)

where NLO calculations include the full $\alpha_{s}^{3}$ radiative NLO
corrections presented by Aversa et. al. in 1988, [14], and LO stands for the
LO Born terms. It is seen that NLO predictions provide a better match to the
data than LO ones, and the NLO corrections are more essential in the small
$p_{t}$ range. The middle panel of Fig. 5 shows the scale dependence of
$\sigma_{pp}^{\pi^{+}}$ at LO and NLO where the scale $Q$ may vary from
$p_{T}/2$ to $2p_{T}$. It is seen that the scale dependence is clearly smaller
at NLO. The same occurs for the polarized cross section for $\pi^{0}$
production, [6].

Figure 5: Left: a comparison of the NLO and LO results for the cross section
$\sigma_{pp}^{\pi^{+}}$ with the use of LSS FFs together with the STAR data.
Middle: the scale dependence of $\sigma_{pp}^{\pi^{+}}$ at LO and NLO for LSS
FFs.. The shaded bands indicate the uncertainties from varying the $Q$ scale
in the range $p_{T}/2\leqslant Q\leqslant 2p_{T}$. The solid (dashed) lines
are for the choice $Q=p_{T}$. For better visibility, the LO results (and the
STAR data, respectively) are rescaled by 0.1. Right: the $K$ factor, Eq. (51),
vs $p_{T}$ for the HKNS and LSS results at the choice $Q=p_{T}$. The error
band for HKNS is shown (dashed).

###  IV.2 Results for charged kaons

We calculate the $p_{T}$ spectra of $K^{+}$ using the parametrizations for
$D_{i}^{K^{+}}$ obtained in HKNS-07, AKK-08 and DSEHS-17, and compare it to
data on $K^{\pm}$-spectra measured by the STAR Collaborations STAR-2007 and
STAR-2012 at $\sqrt{S}=200$ GeV for different ranges of $p_{T}$ and rapidity
$-0.5<y<0.5$ – see Eqs. (41) and (43). All considered parametrizations of
$D_{q}^{K^{+}}$ are obtained within the same theoretical framework – NLO in
QCD111In an interesting study [24], Demirci and Ahmadov have considered
various models of FFs and also the effect of HT in charged kaon production and
suggested that HT effects may be large in certain kinematic regions., using
the simplifying assumption that at the initial scale $Q^{2}_{0}$ all unfavored
fragmentation functions are equal:

$\displaystyle
D_{\bar{u}}^{K^{+}}(z,Q_{0}^{2})=D_{\bar{d}}^{K^{+}}(z,Q_{0}^{2})=D_{d}^{K^{+}}(z,Q_{0}^{2})=D_{s}^{K^{+}}(z,Q_{0}^{2}).$
(52)

If neutral kaons are considered then SU(2) invariance relates the FFs of
$K^{0}$ to those for charged kaons and thus no new FFs are introduced, and, as
shown, this leads to relation (39). If only charged kaons are involved there
are no SU(2) constraints. SU(2) invariance enters the game only when $K^{0}$
are measured as well. Taking into account also charge conjugation invariance
of strong interactions we end up, in total, with 4 independent kaon FFs:

$\displaystyle D_{u}^{K^{+}},\quad D_{\bar{u}}^{K^{+}},\quad
D_{\bar{s}}^{K^{+}},\quad D_{g}^{K^{+}}\cdot$ (53)

Below we give the main characteristics of the FFs used here.

$\bullet$ In HKNS-07 the parametrizations are obtained both in LO and NLO in
QCD from an analysis of the data on $e^{+}e^{-}\to K^{\pm}+X$. We use the NLO
parameters.

$\bullet$ In AKK-08 the parametrizations are obtained in NLO from an analysis
of the data on $e^{+}e^{-}\to K^{\pm}+X$ and $pp\to K^{\pm}+X$. No
uncertainties of the obtained parametrizations are presented.

$\bullet$ In DSEHS-17 the parametrizations are obtained in a global analysis
of the data on $e^{+}e^{-}\to K^{\pm}+X$, $l+N\to K^{\pm}+X$ and $pp\to
K^{\pm}+X$ including the STAR-2012 data. The parameters of the FFs are
presented in NLO, using the same assumption Eq. (52) that all unfavoured FFs
at the initial scale $Q^{2}_{0}$ are the same, but here at NLO.

Our results for kaons are presented in Figs. 6–10 and in Table 2 in a similar
way as those for pions.

Figure 6: The cross sections $\sigma_{pp}^{K^{+}}$ (left) and
$\sigma_{pp}^{K^{-}}$ (right) calculated for different sets of FFs, compared
to STAR-2007 and STAR-2012 data. The one-$\sigma$ uncertainty band for HKNS
FFs estimated by the Hessian method is shown.

Figure 7: The best fit DSEHS together with $\sigma_{pp}^{K^{+}}$ (left) and
$\sigma_{pp}^{K^{-}}$ (right) STAR data. The theoretical uncertainties bands
indicate the scale $Q^{2}$ variation in the range $p_{T}/2\leqslant Q\leqslant
2p_{T}$.

Figure 8: The ratio (data-theory)/theory for the cross sections: $\sigma_{pp}^{K^{+}}$ (left) and $\sigma_{pp}^{K^{-}}$ (right) for different FFs. The experimental error bars and theoretical uncertainties for HKNS (dashed) are shown. Table 2: Values of $\chi^{2}$, Eq. (50), calculated for the single and difference cross sections in the kaon production obtained for different FFs parametrizations. FFs | $~{}~{}~{}K^{+}~{}~{}~{}$ | $~{}~{}~{}K^{-}~{}~{}~{}$ | $~{}~{}K^{+}-K^{-}~{}~{}$
---|---|---|---
HKNS | 124 | 222 | 0.335
AKK | 16.0 | 34.1 | -
DSEHS | 1.18 | 1.31 | 0.329

From Figs. 6–8 and also from the $\chi^{2}$ analysis, Eq. (50), it is seen
that the DSEHS-17 FFs which are based also on the STAR-2012 data provide the
best fit. We took into account the combined STAR data for $p_{T}>1\,{\rm
GeV}^{2}$ with $N=17$ experimental points (STAR 2007 and STAR 2012). In Table
2, we collect the values of $\chi^{2}$, Eq. (50), for the single charged kaons
$K^{\pm}$ and also for the difference $K^{+}-K^{-}$ obtained with the use of
HKNS, AKK and DSEHS sets of FFs.

Figure 9: The difference cross section $\sigma_{pp}^{K^{+}-K^{-}}$ calculated
for different sets of FFs, compared to the STAR data. Left: for the small
$p_{T}$ range. The one-$\sigma$ uncertainty band for HKNS FFs is shown
(dashed). Middle: for $p_{T}>3\,{\rm GeV}/c$. Since the theoretical error for
HKNS is larger than the value of $\sigma_{pp}^{K^{+}-K^{-}}$ itself, we show
only the upper line of the uncertainty (dashed). Right: The ratio (data-
theory)/theory for the difference cross section. The experimental error bars
are shown. For better visibility points for DSEHS are displaced by -0.1 in
$p_{T}$.

A comparison of the difference cross section $\sigma_{pp}^{K^{+}-K^{-}}$
obtained for different sets of FFs with the STAR-2007 and STAR-2012 data is
presented in Fig. 9. For some experimental points of STAR-2007 and also for
the theoretical predictions at $p_{T}<3\,{\rm GeV}/c$ for HKNS and
$p_{T}<4.5\,{\rm GeV}/c$ for DSEHS $\sigma_{pp}^{K^{+}-K^{-}}<0$. Again,
similarly as for pions, the large experimental errors make useless the idea to
find the best fit.

A comparison of the NLO and LO predictions to the data for $\sigma_{pp}^{K+}$
presented in Fig. 10 shows, similarly as in the pion case, better agreement of
the NLO results with the data. From the middle panel one can see the smaller
scale dependence at NLO. The $K$ factor, Eq. (51), obtained for $Q=p_{T}$,
presented in the right panel of Fig. 10, varies between about 2 for
$p_{T}=10\,{\rm GeV}/c$ to 3 for $p_{T}=1\,{\rm GeV}/c$.

Figure 10: Left: a comparison of the NLO and LO results for the cross section
$\sigma_{pp}^{K^{+}}$ with the use of DSEHS FFs together with the STAR data.
Middle: the scale dependence of $\sigma_{pp}^{K^{+}}$ at LO and NLO for DSEHS
FFs. The shaded bands indicate the uncertainties from varying the $Q$ scale in
the range $p_{T}/2\leqslant Q\leqslant 2p_{T}$. The solid (dashed) lines are
for the choice $Q=p_{T}$. For better visibility, the LO results (and the STAR
data, respectively) are rescaled by 0.1. Right: the $K$ factor, Eq. (51), vs
$p_{T}$ for the HKNS and DSEHS results at the choice $Q=p_{T}$. The error band
for HKNS is shown (dashed).

###  IV.3 Results for neutral and charged kaons $K^{+}+K^{-}-2K_{s}^{0}$

It is important to recall that if not only charged, but neutral kaons are
measured, as well, and if SU(2) isospin symmetry holds no new kaon FFs are
introduced. In this section, we shall test SU(2) invariance of the $u$ ad $d$
quarks for kaons.

If the FFs of both charged and neutral kaons are extracted independently from
the data, then they should satisfy Eq. (32). Only in AKK [5] the neutral and
charged kaons FFs are obtained separately, without imposing SU(2)-invariance.
This allows to test Eq. (32) and study the SU(2)- symmetry directly.

Figure 11: Left: Test of the SU(2) invariance via Eq. (32) for AKK FFs.
Middle: A comparison of the AKK predictions with the STAR data for
$K^{+}+K^{-}$ and $2K_{s}^{0}$. Right: [Data - Theory]/Theory for
$K^{+}+K^{-}-2K_{s}^{0}$. The experimental error bars are shown.

In the left panel of Fig. 11, we compare $D_{u}^{K^{+}+K^{-}}$ with
$2D_{d}^{K_{s}^{0}}$ and $D_{d}^{K^{+}+K^{-}}$ with $2D_{u}^{K_{s}^{0}}$. This
clearly shows that isospin symmetry is badly broken for the AKK-08 kaon FFs.
It is also visible in the middle and right panels of the same figure. The
experimental STAR data for kaons where
$\sigma^{K^{+}+K^{-}}=\sigma^{2K_{s}^{0}}$ confirm the SU(2)-invariance. The
theoretical predictions for HKNS and DSEHS based on FFs for which the
SU(2)-invariance is assumed also provides, contrary to AKK,
$\sigma^{K^{+}+K^{-}}-\sigma^{2K_{s}^{0}}=0$.

In summary, the above studies clearly show that one can distinguish two types
of parametrizations of the FFs that describe the $p_{T}$-spectra of the
produced charged kaons in $pp$ collisions: HKNS and DSEHS, for which
SU(2)-invariance holds - that agrees with the data - and AKK for which it does
not. For this reason, making the predictions for NICA, we shall not longer
consider AKK a reasonable fit for kaons.

##  V Results for NICA Project – pions and kaons

In anticipation of the construction of the NICA accelerator we have produced
estimates of various semi-inclusive cross-sections based on the discussed
above FFs. These should be useful for planning experiments.

The NICA kinematics are:

$\displaystyle 5$ $\displaystyle<$ $\displaystyle\sqrt{S}<27\,GeV$ (54)
$\displaystyle 0.1$ $\displaystyle<$ $\displaystyle p_{T}<2\,GeV$ (55)
$\displaystyle 0.1\,\pi$ $\displaystyle<$
$\displaystyle\theta<0.9\,\pi\,\,\to\,\,|\eta|<2\,,$ (56)

where the pseudo-rapidity $\eta$ is defined in Eq. (45).

The predictions for NICA based on the viable sets of FFs elaborated in the
previous section are presented in Figs. 12–16.

Figure 12: Predictions for NICA: $p_{T}$ spectra for $\sigma_{pp}^{\pi^{+}}$
for different parametrizations at different $\sqrt{S}=5$, 10 and 25 GeV and
rapidity $|y|<2$.

Figure 13: Predictions for NICA: the ratio of $p_{T}$ spectra for
$\sigma_{pp}^{\pi^{+}}$ for different parametrizations with respect to HKNS at
different $\sqrt{S}=5$, 10 and 25 GeV and rapidity $|y|<2$.

Figure 14: The $p_{T}$ spectra for $\sigma_{pp}^{\pi^{+}}$ in a linear scale
with visible differences between AKK and other parametrizations for different
$\sqrt{S}$. The error bands for HKNS are shown.

Figure 12 shows the $p_{T}$ spectra of the pion cross section for different
$\sqrt{S}=5$, 10 and 25 GeV. We consider HKNS, AKK, LSS and DSEHS FFs. In
Fig.13, we plot the ratio of $\sigma_{pp}^{\pi^{+}}$ obtained for different
parametrizations with respect to HKNS at the same energies as in Fig.12. In
Fig.14, we present the $p_{T}$ spectra for $\sigma_{pp}^{\pi^{+}}$ in a linear
scale for small $p_{T}$ and at different $\sqrt{S}$. It is seen that the
difference between AKK and other sets of FFs is best pronounced at the smaller
$p_{T}$ for all considered values of $\sqrt{S}$ and increases with the energy.

Figure 15: Predictions for NICA: $p_{T}$ spectra for $\sigma_{pp}^{K^{+}}$ for
HKNS and DSEHS FFs at different $\sqrt{S}=5$, 10 and 25 GeV and rapidity
$|y|<2$. The error bands are shown.

Figure 16: The difference cross section $\sigma_{pp}^{h^{+}-h^{-}}$ cannot
distinguish between different parametrizations of FFs. The error bands for
HKNS are shown.

In Fig. 15, we show the $p_{T}$ spectra of the kaon cross section for
different $\sqrt{S}=5$, 10 and 25 GeV. We consider HKNS and DSEHS FFs. For the
charged kaons, the difference between HKNS and DSEHS FFs are clearly visible
in the whole range of $p_{T}$.

Finally, the $p_{T}$ spectra of the difference cross sections
$\sigma_{pp}^{h^{+}-h^{-}}$ for the pions and kaons are presented in Fig. 16.
One can see that all predictions are indistinguishable due to the theoretical
uncertainties which are relatively large despite of some cancellations. This
effect may be additionally magnified by the large experimental errors.

In summary, measurements at NICA would make it possible to distinguish between
different groups of FFs. Unfortunately, for the difference cross sections the
gap between predictions for different FFs looks undetectably small within the
theoretical and expected experimental error bands.

## Conclusions

There exist in the literature several different relatively up to date
expressions for the pion and kaon fragmentation functions: HKNS-07 [4], AKK-08
[5], DSEHS-14 [1], DSEHS-17 [2] and LSS-15 [3] and we have attempted, in this
paper, to examine which, if any, of these is compatible with the STAR data on
semi-inclusive pion and kaon production in proton-proton collisions. We have
thus compared the predictions for the semi-inclusive cross-sections
$pp\rightarrow\textrm{pions}+X$ and $pp\rightarrow\textrm{kaons}+X$ based on
these FFs with the data from the STAR collaboration. Of the published FFs all
except AKK respect isospin invariance. As seen in Fig. 11, AKK fails badly to
describe the data on $pp\rightarrow K^{+}+K^{-}-2K_{s}^{0}$, so is no longer
considered a viable FF for kaons. All the others provide a good fit to the
STAR data, with LSS for pions and DSEHS for kaons significantly better than
the others.

Furthermore, we have compared the results for the cross sections obtained
within LO and NLO approaches and studied their scale-$Q^{2}$ dependence, where
$Q$ may vary from $p_{T}/2$ to $2p_{T}$. We found that the NLO predictions fit
clearly better the data, specially in the small $p_{T}$ range, and are less
scale dependent than the LO ones.

In our paper, we have also produced estimates of various semi-inclusive cross-
sections based on the above viable FFs for the future run of the NICA
accelerator. It seems that measurements at NICA will have an ability to
determine which of the various FFs versions is the best. For the charged
pions, the differences between AKK and the other FFs are best pronounced at
the smaller $p_{T}$ (see Figs. 12-14). For the charged kaons, the difference
between HKNS and DSEHS FFs are clearly visible in the whole range of $p_{T}$
(see Fig. 15).

We have checked that due to the relatively large uncertainties, measurements
of the difference cross sections, $\sigma^{h^{+}-h^{-}}$, will not be able to
distinguish between different FFs.

## Acknowledgements

E.C. and D.K. acknowledge the support of the INRNE-BAS (Bulgaria) - JINR
(Russia) collaborative Grant, E.C. – of Grant KP-06-N58/5 of the Bulgarian
Science Foundation and D.K. – of the Bogoliubov-Infeld Program. D.K. is
grateful to S.V. Mikhailov and A. Kotlorz for useful discussions on numerical
analysis.

##  Appendix A The partonic subprocesses contributing to $pp\to hX$

Following the paper [9], we obtain the expression for the inclusive cross
section $pp\to hX$:

$\displaystyle
E^{h}\frac{d\sigma_{pp}^{h}}{d^{3}P^{h}}=\frac{1}{\pi}\,\sum_{ab\to
cd}\int^{1}_{x_{a,min}}dx_{a}\int^{1}_{x_{b,min}}dx_{b}\,\frac{1}{z}\,q_{a}(x_{a})q_{b}(x_{b})\,\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{c}^{h}(z)\,.$
(57)

where

$\displaystyle s=x_{a}x_{b}\,S,\quad t=\frac{x_{a}}{z}\,T,\quad
u=\frac{x_{b}}{z}\,U,\quad z=-\,\frac{x_{a}\,T+x_{b}\,U}{x_{a}\,x_{b}\,S}$
(58)

The lower limits of integration are determined by the conditions $z<1$ and
$s+t+u=0$ [9]:

$\displaystyle x_{a,min}$ $\displaystyle=$ $\displaystyle\frac{-U}{T+S}$ (59)
$\displaystyle x_{b,min}$ $\displaystyle=$
$\displaystyle\frac{-x_{a}T}{x_{a}S+U}$ (60)

Hence, we can calculate the contributions from the different partonic
subprocesses of Eq. (57),

$\displaystyle
q_{a}(x_{a})q_{b}(x_{b})\,\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{c}^{h}(z)\,.$
(61)

Below, we give formulas for all possible $ab\to cd$ contributions together
with detailed examples.

The general formula is:

$\displaystyle
q_{a}(x_{a})q_{b}(x_{b})\,\left[\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{c}^{h}(z)+\frac{d\hat{\sigma}_{ab}^{cd}}{du}\,D_{d}^{h}(z)\right]+q_{a}(x_{b})q_{b}(x_{a})\,\left[\frac{d\hat{\sigma}_{ab}^{cd}}{du}\,D_{c}^{h}(z)+\frac{d\hat{\sigma}_{ab}^{cd}}{dt}\,D_{d}^{h}(z)\right],$
(62)

where

  1. 1.

$\hat{\sigma}_{1}$: $q_{i}q_{j}\to
q_{i}q_{j},\quad\bar{q}_{i}\bar{q}_{j}\to\bar{q}_{i}\bar{q}_{j},\quad
q_{i}\bar{q}_{j}\to
q_{i}\bar{q}_{j},\quad\bar{q}_{i}q_{j}\to\bar{q}_{i}q_{j},\quad i\neq j$

$\displaystyle\hat{\sigma}_{1}:~{}~{}~{}~{}~{}~{}\sum_{i\neq j}$
$\displaystyle
q_{i}(x_{a})q_{j}(x_{b})\left[D_{q_{i}}^{h}\,\frac{d\hat{\sigma}_{1}}{dt}+D_{q_{j}}^{h}\,\frac{d\hat{\sigma}_{1}}{du}\right]+\bar{q}_{i}(x_{a})\bar{q}_{j}(x_{b})\left[D_{\bar{q}_{i}}^{h}\,\frac{d\hat{\sigma}_{1}}{dt}+D_{\bar{q}_{j}}^{h}\,\frac{d\hat{\sigma}_{1}}{du}\right]$
$\displaystyle+\,q_{i}(x_{a})\bar{q}_{j}(x_{b})\left[D_{q_{i}}^{h}\,\frac{d\hat{\sigma}_{1}}{dt}+D_{\bar{q}_{j}}^{h}\,\frac{d\hat{\sigma}_{1}}{du}\right]+\bar{q}_{i}(x_{a})q_{j}(x_{b})\left[D_{\bar{q}_{i}}^{h}\,\frac{d\hat{\sigma}_{1}}{dt}+D_{q_{j}}^{h}\,\frac{d\hat{\sigma}_{1}}{du}\right]$
$\displaystyle=\sum_{i\neq j}$
$\displaystyle\left[q_{i}(x_{a})D_{q_{i}}^{h}+\bar{q}_{i}(x_{a})D_{\bar{q}_{i}}^{h}\right]\,\widetilde{q}_{j}(x_{b})\,\frac{d\hat{\sigma}_{1}}{dt}+\left[q_{j}(x_{b})D_{q_{j}}^{h}+\bar{q}_{j}(x_{b})D_{\bar{q}_{j}}^{h}\right]\,\tilde{q}_{i}(x_{a})\,\frac{d\hat{\sigma}_{1}}{du}$
(63) $\displaystyle i\neq j:\,i=u,j=d,s;\,\,i=d,j=u,s,\quad i=s,j=u,d\,.$
$\displaystyle\sigma_{1}^{h-\bar{h}}:~{}~{}~{}~{}~{}\sum_{i\neq j}$
$\displaystyle\left\\{q_{iV}(x_{a})\,\widetilde{q}_{j}(x_{b})\,D_{q_{iV}}^{h}\,\frac{d\hat{\sigma}_{1}}{dt}+q_{jV}(x_{b})\,\widetilde{q}_{i}(x_{a})\,D_{q_{jV}}^{h}\,\frac{d\hat{\sigma}_{1}}{du}\right\\},$
(64) $\displaystyle i\neq j:\,i=u,j=d,s;\,\,i=d,j=u,s,\quad i=s,j=u,d\,.$
$\displaystyle\sigma_{1}^{h+\bar{h}}:~{}~{}~{}~{}~{}=\sum_{i\neq j}$
$\displaystyle\left[\widetilde{q}_{i}(x_{a})\widetilde{q}_{j}(x_{b})\,D_{q_{i}}^{h+\bar{h}}\,\,\frac{d\hat{\sigma}_{1}}{dt}+\widetilde{q}_{j}(x_{b})\,\widetilde{q}_{i}(x_{a})\,D_{q_{j}}^{h+\bar{h}}\,\frac{d\hat{\sigma}_{1}}{du}\right]$
(65) $\displaystyle i\neq j:\,i=u,j=d,s;\,\,i=d,j=u,s,\quad i=s,j=u,d\,.$

  2. 2.

$\hat{\sigma}_{2}$: $q_{i}q_{i}\to
q_{i}q_{i},\quad\bar{q}_{i}\bar{q}_{i}\to\bar{q}_{i}\bar{q}_{i}$

$\displaystyle\hat{\sigma}_{2}:~{}~{}~{}~{}~{}~{}\sum_{i=u,d,s}\left[q_{i}(x_{a})q_{i}(x_{b})\,\,D_{q_{i}}^{h}(z)+\bar{q}_{i}(x_{a})\bar{q}_{i}(x_{b})\,D_{\bar{q}_{i}}^{h}(z)\right]\frac{d\hat{\sigma}_{2}}{dt}\,,$
(66)

symmetric in $t-u$.

$\displaystyle\hat{\sigma}_{2}^{h-\bar{h}}:~{}~{}~{}~{}~{}~{}$
$\displaystyle\sum_{q=u,d,s}\left[q(x_{a})q(x_{b})\,-\bar{q}(x_{a})\bar{q}(x_{b})\right]\,D_{q_{V}}^{h}(z)\frac{d\hat{\sigma}_{2}}{dt}$
(67) $\displaystyle=$
$\displaystyle\sum_{q=u,d,s}\frac{1}{2}\left[q_{V}(x_{a})\widetilde{q}(x_{b})\,+\widetilde{q}(x_{a})\,q_{V}(x_{b})\right]\,D_{q_{V}}^{h}(z)\frac{d\hat{\sigma}_{2}}{dt}\,.$
$\displaystyle\hat{\sigma}_{2}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}$
$\displaystyle\sum_{q=u,d,s}\,\left[\,q(x_{a})q(x_{b})\,+\bar{q}(x_{a})\bar{q}(x_{b})\,\right]\,\frac{d\hat{\sigma}_{2}}{dt}\,D_{q}^{h+\bar{h}}(z)\,.$
(68)

  3. 3.

$\hat{\sigma}_{3}$: $q_{i}\bar{q}_{i}\to q_{j}\bar{q}_{j},\quad i\neq j$

$\displaystyle\hat{\sigma}_{3}:~{}~{}~{}~{}~{}~{}\sum_{i\neq
j}\left[\,q_{i}(x_{a})\bar{q}_{i}(x_{b})+\bar{q}_{i}(x_{a})q_{i}(x_{b})\,\right]\,\frac{d\hat{\sigma}_{3}}{dt}\,D_{j+\bar{j}}^{h}(z)\,.$
(69)

symmetric in $t-u$.

$\displaystyle\hat{\sigma}_{3}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}\sum_{i\neq
j}2\,\left[\,q_{i}(x_{a})\bar{q}_{i}(x_{b})+\bar{q}_{i}(x_{a})q_{i}(x_{b})\,\right]\,\frac{d\hat{\sigma}_{3}}{dt}\,D_{j}^{h+\bar{h}}(z)\,.$
(70)

  4. 4.

$\hat{\sigma}_{4}$: $q_{i}\bar{q}_{i}\to q_{i}\bar{q}_{i}$:

$\displaystyle\hat{\sigma}_{4}:~{}~{}~{}~{}~{}~{}\sum_{i}\left\\{q_{i}(x_{a})\bar{q}_{i}(x_{b})\,\left(\frac{d\hat{\sigma}_{4}}{dt}\,D_{q_{i}}^{h}(z)+\frac{d\hat{\sigma}_{4}}{du}\,D_{\bar{q_{i}}}^{h}(z)\right)+\,\bar{q}_{i}(x_{a})q_{i}(x_{b})\,\left(\frac{d\hat{\sigma}_{4}}{dt}\,D_{\bar{q}_{i}}^{h}(z)+\frac{d\hat{\sigma}_{4}}{du}\,D_{q_{i}}^{h}(z)\right)\right\\}$
$\displaystyle\,.$ (71)
$\displaystyle\hat{\sigma}_{4}^{h-\bar{h}}:~{}~{}~{}~{}~{}~{}$
$\displaystyle\sum_{q=u,d,s}\left[q(x_{a})\bar{q}(x_{b})\,-\,\bar{q}(x_{a})q(x_{b})\right]\,\left(\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right)\,D_{q_{V}}^{h}(z)$
(72)
$\displaystyle=\sum_{q=u,d,s}\left[q_{V}(x_{a}){q}(x_{b})\,-\,{q}(x_{a})q_{V}(x_{b})\right]\,\left(\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right)\,D_{q_{V}}^{h}(z)$
$\displaystyle=\sum_{q=u,d,s}\left[q_{V}(x_{a})\bar{q}(x_{b})\,-\,\bar{q}(x_{a})q_{V}(x_{b})\right]\,\left(\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right)\,D_{q_{V}}^{h}(z)$
$\displaystyle=\sum_{q=u,d,s}\frac{1}{2}\left[q_{V}(x_{a})\tilde{q}(x_{b})\,-\,\tilde{q}(x_{a})q_{V}(x_{b})\right]\,\left(\frac{d\hat{\sigma}_{4}}{dt}-\frac{d\hat{\sigma}_{4}}{du}\right)\,D_{q_{V}}^{h}(z)\,.$
$\displaystyle\hat{\sigma}_{4}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}$
$\displaystyle\sum_{q=u,d,s}\left[q(x_{a})\bar{q}(x_{b})\,+\,\bar{q}(x_{a})q(x_{b})\right]\,\left(\frac{d\hat{\sigma}_{4}}{dt}+\frac{d\hat{\sigma}_{4}}{du}\right)\,D_{q}^{h+\bar{h}}(z)\
$ (73)

  5. 5.

$\hat{\sigma}_{5}$: $q_{i}\bar{q}_{i}\to gg$

$\displaystyle\hat{\sigma}_{5}:~{}~{}~{}~{}~{}~{}\sum_{i=u,d,s}\left[\,q_{i}(x_{a})\bar{q}_{i}(x_{b})+q_{i}(x_{b})\bar{q}_{i}(x_{a})\,\right]\frac{d\hat{\sigma}_{5}}{dt}\,D_{g}^{h}(z)\,,$
(74)

symmetric in $t-u$.

$\displaystyle\hat{\sigma}_{5}^{h-\bar{h}}=0$ (75)
$\displaystyle\hat{\sigma}_{5}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}\sum_{i=u,d,s}2\left[\,q_{i}(x_{a})\bar{q}_{i}(x_{b})+q_{i}(x_{b})\bar{q}_{i}(x_{a})\,\right]\,\frac{d\hat{\sigma}_{5}}{dt}\,D_{g}^{h}(z)\,,$
(76)

  6. 6.

$\hat{\sigma}_{6}$: $gg\to q_{i}\bar{q}_{i}$

$\displaystyle\hat{\sigma}_{6}:~{}~{}~{}~{}~{}~{}g(x_{a})g(x_{b})\,\frac{d\hat{\sigma}_{6}}{dt}\sum_{i=u,d,s}\left[D_{q_{i}}^{h}(z)+D_{\bar{q}_{i}}^{h}(z)\right]\,,$
(77)

symmetric in $t-u$.

$\displaystyle\hat{\sigma}_{6}^{h-\bar{h}}=0$ (78)
$\displaystyle\hat{\sigma}_{6}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}2\,g(x_{a})g(x_{b})\,\frac{d\hat{\sigma}_{6}}{dt}\,\sum_{i=u,d,s}D_{q_{i}}^{h+\bar{h}}(z)\,,$
(79)

where we have used that C-inv. implies $D_{g}^{h}=D_{g}^{\bar{h}}$

  7. 7.

$\hat{\sigma}_{7}$: $q_{i}g\to q_{i}g$

$\displaystyle\hat{\sigma}_{7}:~{}~{}$
$\displaystyle\sum_{i}\left\\{\left[\left(q_{i}(x_{a})D_{q_{i}}^{h}+\bar{q}_{i}(x_{a})D_{\bar{q}_{i}}^{h}\right)\,g(x_{b})\frac{d\hat{\sigma}_{7}}{dt}+\left(q_{i}(x_{a})+\bar{q}_{i}(x_{a})\right)g(x_{b})\,D_{g}^{h}\,\frac{d\hat{\sigma}_{7}}{du}\right]\right.$
(81) $\displaystyle\hskip
14.22636pt+\left.\left[\left(q_{i}(x_{b})D_{q_{i}}^{h}+\bar{q}_{i}(x_{b})D_{\bar{q}_{i}}^{h}\right)\,g(x_{a})\frac{d\hat{\sigma}_{7}}{du}+\left(q_{i}(x_{b})+\bar{q}_{i}(x_{b})\right)g(x_{a})\,D_{g}^{h}\,\frac{d\hat{\sigma}_{7}}{dt}\right]\right\\}\,.$
$\displaystyle\hat{\sigma}_{7}^{h-\bar{h}}:~{}~{}$
$\displaystyle\sum_{q=u,d,s}\left\\{q_{V}(x_{a})\,g(x_{b})\frac{d\hat{\sigma}_{7}}{dt}+q_{V}(x_{b})\,g(x_{a})\frac{d\hat{\sigma}_{7}}{du}\right\\}D_{q_{V}}^{h}\,.$
(82) $\displaystyle\hat{\sigma}_{7}^{h+\bar{h}}:~{}~{}$
$\displaystyle\sum_{q=u,d,s}\left\\{\left[\widetilde{q}(x_{a})\,g(x_{b})\frac{d\hat{\sigma}_{7}}{dt}+\widetilde{q}(x_{b})\,g(x_{a})\frac{d\hat{\sigma}_{7}}{du}\right]D_{q}^{h+\bar{h}}\right.$
(83) $\displaystyle\hskip
34.14322pt+\left.2\left[\widetilde{q}(x_{a})\,g(x_{b})\frac{d\hat{\sigma}_{7}}{du}+\widetilde{q}(x_{b})\,g(x_{a})\frac{d\hat{\sigma}_{7}}{dt}\right]D_{g}^{h}\right\\}\,.$

  8. 8.

$\hat{\sigma}_{8}$: $gg\to gg$

$\displaystyle\hat{\sigma}_{8}:~{}~{}~{}~{}~{}~{}g(x_{a})g(x_{b})\,\frac{d\hat{\sigma}_{8}}{dt}\,D_{g}^{h}(z)\,,$
(84)

symmetric in $t-u$.

$\displaystyle\hat{\sigma}_{8}^{h-\bar{h}}=0$ (85)
$\displaystyle\hat{\sigma}_{8}^{h+\bar{h}}:~{}~{}~{}~{}~{}~{}2\,g(x_{a})g(x_{b})\,\frac{d\hat{\sigma}_{8}}{dt}\,D_{g}^{h}(z)\,,$
(86)

The partonic cross sections that enter the difference cross sections, averaged
over initial and summed over final spin and colour, are [9]:

$\displaystyle\frac{d\hat{\sigma}_{i}(ab\to cd)}{dt}$ $\displaystyle=$
$\displaystyle\frac{\pi\alpha_{s}^{2}(Q^{2})}{s^{2}}\,|M_{i}(s,t,u)|^{2}$ (87)

where the matrix elements $|M_{i}(s,t,u)|^{2}$ are:

$\displaystyle|M_{1}(s,t,u)|^{2}$ $\displaystyle=$
$\displaystyle\,\frac{4}{9}\,\frac{s^{2}+u^{2}}{t^{2}}$ (88)
$\displaystyle|M_{2}(s,t,u)|^{2}$ $\displaystyle=$
$\displaystyle\,\frac{4}{9}\,\left(\frac{s^{2}+u^{2}}{t^{2}}+\frac{s^{2}+t^{2}}{u^{2}}\right)-\frac{8}{27}\,\frac{s^{2}}{tu}$
(89) $\displaystyle|M_{4}(s,t,u)|^{2}$ $\displaystyle=$
$\displaystyle\,\frac{4}{9}\,\left(\frac{s^{2}+u^{2}}{t^{2}}+\frac{u^{2}+t^{2}}{s^{2}}\right)-\frac{8}{27}\,\frac{u^{2}}{st}$
(90) $\displaystyle|M_{7}(s,t,u)|^{2}$ $\displaystyle=$
$\displaystyle\,-\frac{4}{9}\,\frac{s^{2}+u^{2}}{us}+\,\frac{s^{2}+u^{2}}{t^{2}}\cdot$
(91)

## References

  * de Florian et al. [2015] D. de Florian, R. Sassot, M. Epele, R. J. Hernández-Pinto, and M. Stratmann, Phys. Rev. D 91, 014035 (2015), eprint 1410.6027.
  * de Florian et al. [2017] D. de Florian, M. Epele, R. J. Hernandez-Pinto, R. Sassot, and M. Stratmann, Phys. Rev. D 95, 094019 (2017), eprint 1702.06353.
  * Leader et al. [2016] E. Leader, A. V. Sidorov, and D. B. Stamenov, Phys. Rev. D 93, 074026 (2016), eprint 1506.06381.
  * Hirai et al. [2007] M. Hirai, S. Kumano, T. H. Nagai, and K. Sudoh, Phys. Rev. D 75, 094009 (2007), eprint hep-ph/0702250.
  * Albino et al. [2008] S. Albino, B. A. Kniehl, and G. Kramer, Nucl. Phys. B 803, 42 (2008), eprint 0803.2768.
  * Jager et al. [2003] B. Jager, A. Schafer, M. Stratmann, and W. Vogelsang, Phys. Rev. D 67, 054005 (2003), eprint hep-ph/0211007.
  * Combridge et al. [1977] B. L. Combridge, J. Kripfganz, and J. Ranft, Phys. Lett. B 70, 234 (1977).
  * Cutler and Sivers [1978] R. Cutler and D. W. Sivers, Phys. Rev. D 17, 196 (1978).
  * Owens et al. [1978] J. F. Owens, E. Reya, and M. Gluck, Phys. Rev. D 18, 1501 (1978).
  * Feynman et al. [1978] R. P. Feynman, R. D. Field, and G. C. Fox, Phys. Rev. D 18, 3320 (1978).
  * Baier et al. [1979] R. Baier, J. Engels, and B. Petersson, Z. Phys. C 2, 265 (1979).
  * Ellis et al. [1980] R. K. Ellis, M. A. Furman, H. E. Haber, and I. Hinchliffe, Nucl. Phys. B 173, 397 (1980).
  * Ellis and Sexton [1986] R. K. Ellis and J. C. Sexton, Nucl. Phys. B 269, 445 (1986).
  * Aversa et al. [1989] F. Aversa, P. Chiappetta, M. Greco, and J. P. Guillet, Nucl. Phys. B 327, 105 (1989).
  * Adams et al. [2006] J. Adams et al. (STAR), Phys. Lett. B 637, 161 (2006), eprint nucl-ex/0601033.
  * Abelev et al. [2007] B. I. Abelev et al. (STAR), Phys. Rev. C 75, 064901 (2007), eprint nucl-ex/0607033.
  * Agakishiev et al. [2012] G. Agakishiev et al. (STAR), Phys. Rev. Lett. 108, 072302 (2012), eprint 1110.0579.
  * Pumplin et al. [2002] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky, and W. K. Tung, JHEP 07, 012 (2002), eprint hep-ph/0201195.
  * Christova and Leader [2009] E. Christova and E. Leader, Phys. Rev. D 79, 014019 (2009), eprint 0809.0191.
  * Bourrely et al. [2007] C. Bourrely, J. Soffer, and F. Buccella, Phys. Lett. B 648, 39 (2007), eprint hep-ph/0702221.
  * Torii [2003] H. Torii (PHENIX), Nucl. Phys. A 715, 753 (2003), eprint nucl-ex/0210005.
  * [22] P. Aurenche, T. Binoth, M. Fontannaz, J.-P. Guillet, G. Heinrich, E. Pilon, and M. Werlen, https://lapth.cnrs.fr/PHOX_FAMILY/readme_inc.html.
  * Pumplin et al. [2001] J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Phys. Rev. D 65, 014013 (2001), eprint hep-ph/0101032.
  * Demirci and Ahmadov [2018] M. Demirci and A. I. Ahmadov, Int. J. Mod. Phys. A 33, 1850166 (2018), eprint 1810.01491.


Sensor Applications

# Detecting Anomaly in Chemical Sensors via L1-Kernel based Principal
Component Analysis

Hongyi Pan1, Diaa Badawi1, Ishaan Bassi2, Sule Ozev2, Ahmet Enis Cetin1
11Department of Electrical and Computer Engineering, University of Illinois
Chicago, Chicago, IL, 60607, USA
2School of Electrical, Computer and Energy Engineering, Arizona State
University, Tempe, AZ, 85281, USA
1Fellow, IEEE This work is supported by NSF grants 1739396, 1934915 (UIC) and
1739451 (ASU). Corresponding author is H. Pan (e-mail: hpan21@uic.edu).

###### Abstract

We propose a kernel-PCA based method to detect anomaly in chemical sensors. We
use temporal signals produced by chemical sensors to form vectors to perform
the Principal Component Analysis (PCA). We estimate the kernel-covariance
matrix of the sensor data and compute the eigenvector corresponding to the
largest eigenvalue of the covariance matrix. The anomaly can be detected by
comparing the difference between the actual sensor data and the reconstructed
data from the dominant eigenvector. In this paper, we introduce a new
multiplication-free kernel, which is related to the $\ell_{1}$-norm for the
anomaly detection task. The $\ell_{1}$-kernel PCA is not only computationally
efficient but also energy-efficient because it does not require any actual
multiplications during the kernel covariance matrix computation. Our
experimental results show that our kernel-PCA method achieves a higher area
under curvature (AUC) score (0.7483) than the baseline regular PCA method
(0.7366).

Anomalous sensor detection, Principal component analysis, $\ell_{1}$-kernel,
multiplication-free method.

##  1 Introduction

Chemical sensors are widely used to detect ammonia, methane, and other
Volatile Organic Compounds (VOCs) [1, 2, 3]. The life and performance of
chemical gas detection sensors can be affected by various factors, including
temperature, humidity, other interfering chemical gases, physical factors etc.
Anomalous sensors can produce drifting waveforms and it is a fatal problem for
reliable gas identification and concentration estimation [4, 5]. In this work,
we determine anomalous sensors and sensor measurements in an array of
uncalibrated sensors by using robust $\ell_{1}$-Principal Component Analysis
(PCA) without using a reference time-series data.

PCA is used in anomalous sensor and sensor signal detection [6, 7, 8]. In this
approach, the covariance matrix is constructed from a set of data vectors and
the anomalous items (outliers) or vectors are found by using the
reconstruction difference [9, 10]. The principal components of the data
covariance matrix is computed and the original data vectors are reconstructed
using only the first few principal components. In general, the reconstructed
data vector is similar to the original data vectors, and the reconstructed
data items that are different from the corresponding original items are
considered to be anomalous.

In this paper, we propose to use $\ell_{1}$-kernel PCA based on a
multiplication-free kernel to detect anomaly in chemical sensors. Although the
conventional PCA, which is based on the $\ell_{2}$-norm has successfully
solved many problems it is sensitive to outliers in the data because the
effects of the outliers are over-amplified by the $\ell_{2}$-norm. Recently,
it has been shown that $\ell_{1}$-norm based methods produce better results in
practical problems compared to the $\ell_{2}$-norm-based methods in several
real-world signal, image, and video processing problems. In particular,
$\ell_{1}$-kernel PCA usually is more robust against outliers in data compared
to the $\ell_{2}$-PCA [11]. In the ordinary $\ell_{2}$-PCA, the principal
vector is calculated as the dominant eigenvector of the data covariance
matrix, which itself is calculated using the standard outer products. In this
paper, we propose the eigen-decomposition of the $\ell_{1}$-kernel covariance
matrix obtained using a new vector product, which induces the $\ell_{1}$-norm
without performing any multiplications. Because of low computational
complexity, the $\ell_{1}$-kernel PCA can be implemented in edge devices
directly connected to the chemical sensors.

Related work includes the regular $\ell_{2}$-PCA, kernel-PCA methods [12] and
the recursive $\ell_{1}$-PCA [13]. The recursive $\ell_{1}$-PCA method
requires some parameters to be properly adjusted. On the contrary, the
proposed $\ell_{1}$-kernel PCA approach does not need any hyper-parameters
adjustments. Its implementation is as straightforward as the regular PCA
because we construct a sample kernel-covariance matrix using the proposed
$\ell_{1}$-kernel and obtain the eigenvalues and eigenvectors of the kernel
covariance matrix to define the linear transformation instead of solving an
optimization problem.

The rest of the paper is organized as follows: In Section 2, we formally
introduce the $\ell_{1}$-kernel PCA and describe its use in an anomalous
chemical sensor detection task. In Section 3, we compare our methods with the
regular PCA and the recursive $\ell_{1}$-PCA [13]. Finally, in Section 4, we
draw our main conclusions.

##  2 L1-Kernel PCA

In our recent work [14, 11], we proposed a set of kernel-based PCA methods
related with the $\ell_{1}$-norm. These Mercer-type kernels are obtained from
multiplication-free (MF) dot products.

Let $\mathbf{w}=[w_{1}\cdots w_{n}]^{T}\in\mathbb{R}^{n\times 1}$ and
$\mathbf{x}=[x_{1}\cdots x_{n}]^{T}\in\mathbb{R}^{n\times 1}$ be two
$n$-dimensional column vectors. Similar to the regular dot product of vectors
we defined the multiplication-free (MF) vector product [11]. In vector data
correlation operations, we use the following vector product

$\displaystyle\mathbf{w}^{T}\odot\mathbf{x}$
$\displaystyle\triangleq\sum_{i=1}^{n}\text{sign}(w_{i}\times
x_{i})\min(|w_{i}|,|x_{i}|).$ (1)

which turns out to be a Mercer type kernel:
$K(\mathbf{w},\mathbf{x}))=\mathbf{w}^{T}\odot\mathbf{x}$. In Eq. (1),
$\text{sign}(w_{i}\times x_{i})$ can be computed without performing any actual
multiplications and $\min$ operation can be implemented by subtraction and
checking the sign of the result of the subtraction. For this reason, we call
(1) a Multiplication-Free (MF) dot product. The dot product defined in Eq. (1)
induces the $\ell_{1}$-norm as
$\mathbf{x}^{T}\odot\mathbf{x}=\sum_{i=1}^{n}\min(|x_{i}|,|x_{i}|)=\|\mathbf{x}\|_{1}$
and it induces a Mercer-type kernel [11].

Suppose that we collect vectors of sensor data and form a dataset
$\mathbf{X}=[\mathbf{x}_{1}\ \mathbf{x}_{2}\ ...\
\mathbf{x}_{D}]\in\mathbb{R}^{N\times D}$. The well-known $\ell_{2}$-PCA
(regular PCA) method relies on the eigendecomposition of the sample covariance
matrix $\mathbf{C}=\mathbf{X}^{T}\mathbf{X}$. Similarly, we estimate the
kernel covariance matrix as follows:

$\displaystyle\mathbf{A}=\mathbf{X}^{T}\odot\mathbf{X}$ (2)

where the matrix $\mathbf{A}$ is constructed using the dot products of the
form $\mathbf{x}_{i}^{T}\odot\mathbf{x}_{j}$. As a result, the construction of
the kernel-covariance matrix $\mathbf{A}$ is straightforward. We name the
kernel-PCA based on the vector product in Eq. (1) as the $\ell_{1}$-kernel
PCA.

###  2.1 Anomaly Detection Using $L_{1}$-Kernel PCA

We first assume that there are $D$ sensors and some of them are anomalous.
Sensors are assumed to be close to each other and they produce correlated
output waveforms as shown in Fig. 1. We have the measurement data
$\mathbf{X}=[\mathbf{x}_{1}\ \mathbf{x}_{2}\ ...\
\mathbf{x}_{D}]\in\mathbb{R}^{N\times D}$. The measurement data is generated
by normalizing the raw sensor measurement values to $[-1,1]$ and subtracting
the mean. We construct the covariance matrix
$\mathbf{A}=\mathbf{X}^{T}\odot\mathbf{X}\in\mathbb{R}^{D\times D}$ and
calculate its eigenvectors. Let $\mathbf{v}_{1}\in\mathbb{R}^{D\times 1}$ be
the dominant principal component vector. Finally, we reconstruct the data
segment using the vector $\mathbf{v}_{1}$:
$\hat{\mathbf{x}}_{i}=\mathbf{v}_{1}\mathbf{v}_{1}^{T}\mathbf{x}_{i},i=1,2,..,D$
and compute the error vector ${\mathbf{x}}_{i}-\hat{\mathbf{x}}_{i}$. The
sensor measurement in the $i$-th segment is assumed to be anomalous if the
Cumulative Squared Difference (CSD) between $\hat{\mathbf{x}}_{i}$ and
${\mathbf{x}}_{i}$ is larger than a threshold. The threshold can be set as
$T=\mu+\alpha\sigma$, where $\mu$ and $\sigma$ are the mean and standard
deviation of CSD values learned from a training data set. The parameter
$\alpha$ is usually selected as 3 with the Gaussianity assumption of the CSD
values. If a sensor produces successive anomalous measurement vectors it is
considered to be anomalous.

(a) Sample data 1

(b) Sample data 2

Figure 1: Data produced by three ammonia sensors in two different ammonia gas
recording experiments. The sensor 2 (in red) is the obstructed sensor. This
data constitute our test set.

In the second case, we assume that we only have a single sensor. For example,
the sensor 3 produces impulsive spikes between 0 and 1200$s$ as shown in
Fig.1a (in orange). Similar to the PCA-based denoising methods [12, 15] we
form data vectors from neighboring temporal data windows and form the
measurement data matrix $\mathbf{X}=[\mathbf{x}_{1}\ \mathbf{x}_{2}\ ...\
\mathbf{x}_{L}]\in\mathbb{R}^{N\times L}$ where $L$ is the number of data
windows and we have $N$ samples in each window. We form the kernel covariance
matrix $\mathbf{A}=\mathbf{X}^{T}\odot\mathbf{X}\in\mathbb{R}^{L\times L}$ and
compute its eigenvalues and eigenvectors $\mathbf{v}_{i}\in\mathbb{R}^{L\times
1}$. We reconstruct the data matrix
$\hat{\mathbf{X}}=\mathbf{V}\mathbf{V}^{T}{\mathbf{X}}$ using the first $l$
eigenvectors $\mathbf{V}=[{\bf v_{1},...,v}_{l}]$ where $\mathbf{v}_{1}$ is
the eigenvector corresponding to the largest eigenvalue. After this step we
compare the actual data vectors $\mathbf{x}_{i}$ with the reconstructed ones
$\hat{\mathbf{x}}_{i},i=1,2,...,N$. The vectors significantly different from
the reconstructed vectors are considered to be anomalous. We let $l=1$ and
observed that it is sufficient for anomaly detection.

Complexity Analysis: To calculate the covariance matrix $\mathbf{C}$, we
perform $D^{2}N$ multiplications and $D^{2}(N-1)$ additions because we perform
$N$ multiplications in each dot product. On the other hand, to calculate the
$\ell_{1}$-kernel covariance matrix $\mathbf{A}$, we perform $D^{2}N$ sign
operations, $D^{2}N$ min operations and $D^{2}(N-1)$ additions. A
multiplication operation consumes about 4 times more energy compared to the
MF-operations (Eq. (1)) in compute-in-memory (CIM) implementation at 1 GHz
operating frequency [16], [17]. Since the value of $D=3$ or $5$ is much
smaller than the vector length $N=125$ or $224$ the eigenvector computation is
negligible compared to the covariance matrix construction in our applications.
As a result, the $\ell_{1}$-kernel PCA is about 4 times more energy efficient
in CIM implementation. It is also more energy efficient in many other
processors because multiplications consume more energy than additions and
subtractions.

##  3 Experimental Results

We collected data using three ammonia MQ137 Tin oxide (SnO2) based sensors
[18]. Sensors are connected to an Arduino Uno board, and the sampling rate is
2 samples per second. Sensors and a cylindrical ammonia source are placed in
an airtight chamber. Sensors are pre-heated for 48 hours before collecting the
data. When SnO2 is heated and exposed to the air, it reacts with the oxygen
present in the air and form a layer of negative ion on the surface and reduce
the surface conductivity [18]. When ammonia vapor comes in contact with the
surface, it combines with the oxide ion layer on the top and releases
electrons for conduction. As a result, the conductivity of the surface
increases. This change in surface resistance can be measured in the form of
voltage.

Multiple Sensor Anomaly Detection: The three sensors are placed close to each
other and one of the sensors (sensor 2) is obstructed with a cylindrical cover
with multiple holes. The cover causes the sensor to react more slowly than the
other sensors to the ammonia build-up and release. The obstruction level of
the outlier sensor is adjusted in each trial to avoid over-fitting to one
condition. Moreover, to generate a more realistic environment with varying
levels of ammonia concentration, the chamber lid is opened at random intervals
and with random duration. Opening and closing the lid is repeated multiple
times to create different rise and fall responses. Sensor waveforms from two
experiments are shown in Fig. 1.

We apply the $\ell_{1}$-kernel PCA based anomaly detection method described in
Section 2.1 to the data obtained from the three sensors. We compared the
proposed method with the regular $\ell_{2}$-PCA and the recursive $L_{1}$-PCA
[13], which computes
$\mathbf{v}_{1}=\arg\max_{\mathbf{v}:\|\mathbf{v}\|=1}\textstyle\sum_{i=1}^{N}|\mathbf{v}^{T}\mathbf{x}_{i}|$.
Tolerance parameter of the recursive $\ell_{1}$-PCA method is set as $1\times
10^{-8}$ as suggested by the authors [13]. We plot the Receiver Operating
Characteristic (ROC) curve in Fig. 2 and compute the Area Under Curve (AUC)
scores for each method. As shown in Table 1 states, the $\ell_{1}$-kernel PCA
provides the highest AUC score, and the recursive $\ell_{1}$-PCA [13] provides
the lowest AUC score in this case. This is probably due to the non-convex
optimization method that they use to compute the principle vector, and the
method requires a suitable tolerance parameter. On the other hand, the
$\ell_{1}$-kernel PCA does not need any tolerance parameters and the
eigenvector computations are equivalent to the computational load of the
regular PCA.

Figure 2: ROC of different PCAs on the sensor data. Method | AUC score | AUC Increment
---|---|---
Regular PCA (baseline) | 0.7366 | -
Recursive $\ell_{1}$-PCA [13] | 0.7203 | -0.0163
$L_{1}$-kernel PCA | 0.7483 | 0.0117
Table 1: Area under curvature (AUC) for various methods. The $\ell_{1}$-kernel
PCA provides the best AUC score.

CSD values of the sensors’ response in Fig. 1a using the regular PCA and the
$\ell_{1}$-kernel PCA are listed in Tables 2 and 3, and CSD values of the
sensors response in Fig. 1b using the regular PCA and the $\ell_{1}$-kernel
PCA are listed in Tables 4 and 5, respectively. The threshold values 21.63 for
the regular PCA and 21.61 for the $\ell_{1}$-kernel PCA are computed from the
normal cases obtained from other experiments. We have different threshold
values for the regular PCA and the $\ell_{1}$-kernel PCA because they have
different principle eigenvectors. Even if we use the same threshold the
classification result in Tables 4 and 5 will not change.

Segment (Time/S) | Sensor 1 | Sensor 2 | Sensor 3
---|---|---|---
0–249 | 7.54 | 9.44 | 2.20
250–499 | 25.61 | 18.48 | 9.23
500–749 | 11.81 | 12.19 | 17.13
750–999 | 0.71 | 24.75 | 1.08
1000–1249 | 3.03 | 9.11 | 11.88
1250–1499 | 4.57 | 7.02 | 1.30
1500–1749 | 5.71 | 31.50 | 2.79
1750–1999 | 3.77 | 28.11 | 3.52
2000–2249 | 1.29 | 14.74 | 2.04
2250–2499 | 2.09 | 26.11 | 3.55
2500–2749 | 5.08 | 6.08 | 7.39
2750–2999 | 1.09 | 10.67 | 1.87
3000–3249 | 0.98 | 12.89 | 3.51
Table 2: CSD of the sample data 1 in Fig. 1a using regular PCA. Sensor 2 is obstructed. Values larger than threshold 21.63 are in bold. Segment (Time/S) | Sensor 1 | Sensor 2 | Sensor 3
---|---|---|---
0–249 | 8.01 | 9.08 | 2.21
250–499 | 19.54 | 27.64 | 7.24
500–749 | 14.06 | 10.96 | 16.59
750–999 | 0.75 | 25.26 | 0.90
1000–1249 | 3.13 | 8.70 | 12.35
1250–1499 | 4.64 | 6.92 | 1.38
1500–1749 | 7.84 | 30.97 | 1.73
1750–1999 | 3.32 | 29.27 | 3.01
2000–2249 | 1.46 | 14.61 | 2.13
2250–2499 | 1.98 | 25.94 | 3.88
2500–2749 | 4.95 | 6.32 | 7.37
2750–2999 | 1.17 | 10.66 | 1.86
3000–3249 | 0.96 | 12.94 | 3.51
Table 3: CSD of the sample data 1 in Fig. 1a using the $\ell_{1}$-kernel PCA.
Sensor 2 is obstructed. Values larger than threshold 21.61 are in bold.

The AUC scores are relatively low in Table 1 because Sensor 2 does not always
exhibit anomalous behavior. In general, its response increases due to ammonia
gas exposure but not decrease as fast as the other sensors when there is no
gas as shown in Fig. 1. The regular PCA detects the anomalous behavior of
Sensor 2 in 9 out of 28 data segments. On the other hand the $\ell_{1}$-kernel
PCA detects the anomalous behavior of Sensor 2 in 10 out of 28 data segments.
Moreover, the regular PCA produces a false alarm in the second data segment
(250 - 499 seconds) of Sensor 1 as shown in Table 4. In conclusion,
$\ell_{1}$-kernel PCA produces better results than the regular PCA and
Recursive $\ell_{1}$-PCA [13].

Anomaly Detection Using a Single Sensor: During the first three ammonia gas
exposures, the sensor 3 positively responds but it also produces spikes up to
1200s as shown in Fig. 1a. Multisensor PCA-based anomaly detection cannot
detect this behavior because only the sensor 1 works properly before 1200s.
However, we compare the current sensor window of Sensor 2 with its neighboring
data windows we can identify the anomalous behavior. We used $L=5$ data
segments to construct the $5\times 5$ covariance and $\ell_{1}$-kernel
covariance matrices. In each data segment we have $N=224$ measurement. We used
only the dominant eigenvector to estimate the data segments.

Segment (Time/S) | Sensor 1 | Sensor 2 | Sensor 3
---|---|---|---
0–249 | 1.82 | 11.52 | 1.32
250–499 | 2.55 | 1.30 | 1.06
500–749 | 13.63 | 11.07 | 6.61
750–999 | 4.62 | 6.93 | 1.36
1000–1249 | 3.29 | 1.76 | 1.85
1250–1499 | 4.76 | 50.56 | 1.51
1500–1749 | 1.68 | 28.23 | 0.92
1750–1999 | 0.988 | 21.09 | 3.91
2000–2249 | 2.65 | 18.01 | 18.88
2250–2499 | 5.67 | 4.72 | 1.41
2500–2749 | 3.62 | 1.58 | 0.61
2750–2999 | 1.97 | 1.55 | 1.81
3000–3249 | 5.54 | 50.49 | 1.79
3250–3499 | 3.58 | 37.25 | 1.29
3500–3749 | 2.53 | 35.25 | 4.92
Table 4: CSD of the sample data 2 in Fig. 1b using regular PCA. Sensor 2 is obstructed. Values larger than threshold 21.63 are in bold.  Segment (Time/S) | Sensor 1 | Sensor 2 | Sensor 3
---|---|---|---
0–249 | 2.09 | 11.29 | 1.45
250–499 | 2.52 | 1.35 | 1.06
500–749 | 13.55 | 10.91 | 6.93
750–999 | 4.37 | 7.42 | 1.25
1000–1249 | 3.25 | 1.81 | 1.86
1250–1499 | 5.12 | 50.12 | 1.92
1500–1749 | 2.10 | 27.71 | 1.20
1750–1999 | 1.22 | 20.74 | 4.20
2000–2249 | 2.41 | 20.60 | 17.05
2250–2499 | 5.52 | 4.95 | 1.39
2500–2749 | 3.57 | 1.65 | 0.62
2750–2999 | 1.98 | 1.61 | 1.78
3000–3249 | 5.72 | 50.23 | 1.98
3250–3499 | 3.13 | 38.00 | 1.06
3500–3749 | 2.42 | 35.53 | 4.78
Table 5: CSD of the sample data 2 in Fig. 1b using the $\ell_{1}$-kernel PCA.
Sensor 2 is obstructed. Values larger than threshold 21.61 are in bold.

The regular PCA and the $\ell_{1}$-kernel PCA reconstructed waveforms do not
have spikes and that is how we can identify the anomaly in sensor readings as
shown in Fig. 3. Table 6 shows the CSD values of both methods and they
correctly identified the anomolous segments.

Segment | Regular PCA |  $L_{1}$-Kernel PCA
---|---|---
1 | 9.62 | 9.79
2 | 13.59 | 13.55
3 | 7.21 | 7.26
4 | 4.35 | 4.38
5 | 3.48 | 3.65

Table 6: CSD values due to ammonia exposure of Sensor 3. Threshold values are
5.50 for the regular PCA and 5.64 for the $\ell_{1}$-kernel PCA.

.

(a) Segment 1

(b) Segment 2

(c) Segment 3

Figure 3: Response of the sensor 3 to ammonia gas exposures in orange. The
first three segments contains impulsive spikes. The regular PCA and the
$\ell_{1}$-kernel PCA reconstructed waveforms do not have spikes.

##  4 Conclusion

In this paper, we presented a framework for detecting anomalous sensors and
sensor measurements in a chemical sensory system using the $\ell_{1}$-kernel
PCA. We collected data from three commercial Tin Oxide (SnO2) sensors by
exposing them to ammonia in an environment-controlled experiment. The proposed
$\ell_{1}$-kernel PCA is more robust than the regular PCA in our experiments.
This is due to the fact that $\ell_{1}$-kernel PCA is related with the
$\ell_{1}$-norm and it gives less emphasis to anomalous spikes in sensor
measurements while computing the correlation matrix. The computational energy
cost of the $\ell_{1}$-kernel PCA is much lower than the regular PCA on many
processors. Because of low energy complexity, the $\ell_{1}$-kernel PCA can be
implemented in low-cost edge devices directly connected to the chemical
sensors.

## References

  * [1] Alexander Vergara, Shankar Vembu, Tuba Ayhan, Margaret A Ryan, Margie L Homer, and Ramón Huerta.  Chemical gas sensor drift compensation using classifier ensembles.  Sensors and Actuators B: Chemical, 166:320–329, 2012.
  * [2] Fatih Erden, E Birey Soyer, B Ugur Toreyin, and A Enis Cetin.  Voc gas leak detection using pyro-electric infrared sensors.  In 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 1682–1685. IEEE, 2010.
  * [3] Diaa Badawi, Hongyi Pan, Sinan Cem Cetin, and A Enis Çetin.  Computationally efficient spatio-temporal dynamic texture recognition for volatile organic compound (voc) leakage detection in industrial plants.  IEEE Journal of Selected Topics in Signal Processing, 14(4):676–687, 2020.
  * [4] U Mittal, Jitender Kumar, AT Nimal, MU Sharma, et al.  Single chip readout electronics for saw based gas sensor systems.  In 2017 IEEE SENSORS, pages 1–3. IEEE, 2017.
  * [5] Hilmi E Egilmez and Antonio Ortega.  Spectral anomaly detection using graph-based filtering for wireless sensor networks.  In 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1085–1089. IEEE, 2014.
  * [6] Ling Huang, XuanLong Nguyen, Minos Garofalakis, Michael I Jordan, Anthony Joseph, and Nina Taft.  In-network pca and anomaly detection.  In NIPS, volume 2006, pages 617–624, 2006.
  * [7] Vassilis Chatzigiannakis and Symeon Papavassiliou.  Diagnosing anomalies and identifying faulty nodes in sensor networks.  IEEE Sensors Journal, 7(5):637–645, 2007.
  * [8] Laura Erhan, M Ndubuaku, Mario Di Mauro, Wei Song, Min Chen, Giancarlo Fortino, Ovidiu Bagdasar, and Antonio Liotta.  Smart anomaly detection in sensor systems: A multi-perspective review.  Information Fusion, 67:64–79, 2021.
  * [9] Yuan Yao, Abhishek Sharma, Leana Golubchik, and Ramesh Govindan.  Online anomaly detection for sensor systems: A simple and efficient approach.  Performance Evaluation, 67(11):1059–1075, 2010.
  * [10] A Sree and K Venkata.  Anomaly detection using principal component analysis.  Journal of Computer Science and Technology, 5 (4), pages 124–126, 2014.
  * [11] Hongyi Pan, Diaa Badawi, Erdem Koyuncu, and A. Enis Cetin.  Robust principal component analysis using a novel kernel related with the $\ell_{1}$-norm.  In 2021 29th European Signal Processing Conference (EUSIPCO), pages 2189–2193, 2021.
  * [12] B. Schölkopf, A. Smola, and K. Müller.  Kernel principal component analysis.  In Artificial Neural Networks — ICANN’97, Lecture Notes in Computer Science, vol 1327, Heidelberg, Berlin, Germany, 1997.
  * [13] Panos P Markopoulos, George N Karystinos, and Dimitris A Pados.  Optimal algorithms for $l_{1}$-subspace signal processing.  IEEE Transactions on Signal Processing, 62(19):5046–5058, 2014\.
  * [14] Cem Emre Akbaş, Osman Günay, Kasım Taşdemir, and A Enis Çetin.  Energy efficient cosine similarity measures according to a convex cost function.  Signal, Image and Video Processing, 11(2):349–356, 2017.
  * [15] SB Mike, B Scholkopf, and AJ Smola.  Kernel pca and denoising in feature space.  Advances in Neural Information Processing System, pages 524–536, 1999.
  * [16] Shamma Nasrin, Diaa Badawi, Ahmet Enis Cetin, Wilfred Gomes, and Amit Ranjan Trivedi.  Mf-net: Compute-in-memory sram for multibit precision inference using memory-immersed data conversion and multiplication-free operators.  IEEE Trans. Circuits and Systems, 2021.
  * [17] Shamma Nasrin, Ahish Shylendra, Yuti Kadakia, Nick Iliev, Wilfred Gomes, Theja Tulabandhula, and Amit Ranjan Trivedi.  Enos: Energy-aware network operator search for hybrid digital and compute-in-memory dnn accelerators.  arXiv preprint arXiv:2104.05217, 2021.
  * [18] Joseph Watson, Kousuke Ihokura, and Gary SV Coles.  The tin dioxide gas sensor.  Measurement Science and Technology, 4(7):711, 1993.