# Spin Polarization Formula for Dirac Fermions at Local Equilibrium

Yu-Chen Liu  Physics Department and Center for Field Theory and Particle
Physics, Fudan University, Shanghai 200433, China  Xu-Guang Huang
<EMAIL_ADDRESS>Physics Department and Center for Field Theory and
Particle Physics, Fudan University, Shanghai 200433, China  Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai
200433, China

###### Abstract

We derive a Cooper-Frye type spin polarization formula for Dirac fermions at
local thermal equilibrium described by a grand canonical ensemble specified by
temperature, fluid velocity, chemical potential, and spin potential. We
discuss the physical meaning of different contributions to spin polarization
and compare them with previous works. The present formula provides machinery
to convert the spin potential computed in, e.g., relativistic spin
hydrodynamics to the spin polarization observable in, e.g., heavy-ion
collisions.

Keywords: Spin polarization, Local equilibrium, Spin kinetic theory

###### pacs:

12.38.Mh, 25.75.-q, 05.30.-d

##  I Introduction

The phenomena of spin polarization attract broad attention in different
branches of physics from condensed matter physics to particle physics to
nuclear physics. In low-energy nuclear physics, the spin polarization is an
important characteristic of the nuclear equation of state Fantoni et al.
(2001); Vidana and Bombaci (2002); Isayev and Yang (2004); Stein et al.
(2016). In heavy-ion collisions, the measurement of global spin polarization
of $\Lambda$ (and $\bar{\Lambda}$) hyperons with respect to the direction of
the reaction plane has been reported recently Adamczyk et al. (2017); Adam et
al. (2018); Abdallah et al. (2021) and that of $\Xi^{-}$ and $\Omega^{-}$
hyperons was reported in 2021 Adam et al. (2021).

In a fluid composed of massive spin-1/2 Dirac fermions at global equilibrium,
a Cooper-Frye type spin polarization formula has been derived, which shows
that the mean spin vector $\bar{S}^{\mu}$ at global equilibrium is determined
solely by thermal vorticity (assuming that other external torques, like the
magnetic field, do not exist) Becattini et al. (2013); Fang et al. (2016); Liu
et al. (2020a):

$\displaystyle\bar{S}_{\mu}({\bm{p}})=-\frac{1}{8m}\epsilon_{\mu\nu\alpha\beta}p^{\nu}\frac{\int
d\Xi\cdot p\,n_{F}\left(1-n_{F}\right)\varpi^{\alpha\beta}}{\int d\Xi\cdot
p\,n_{F}}\,,\;\;\;$ (1)

where $p^{\mu}$ is the on-shell 4-momentum, $\Xi$ is a spacelike hypersurface
on which the particles are frozen out, $n_{F}=n_{F}(p\cdot\beta)$ with
$n_{F}(x)=1/[\exp(x)+1]$ is the Fermi-Dirac function, and

$\displaystyle\varpi_{\mu\nu}=-\partial_{[\mu}\beta_{\nu]}\,\;\;\;$ (2)

is the thermal vorticity with $X_{[\mu\nu]}\equiv(X_{\mu\nu}-X_{\nu\mu})/2$
and $\beta^{\mu}=\beta u^{\mu}$ ($\beta=1/T$ with $T$ temperature) the thermal
flow velocity. The spin polarization vector is given by

$\displaystyle
P_{\mu}({\bm{p}})=\frac{1}{s}\bar{S}_{\mu}({\bm{p}})=2\bar{S}_{\mu}({\bm{p}})\,,\;\;\;$
(3)

where $s=1/2$ is the spin quantum number. Based on this formula, the
experimental data for global spin polarization of hyperons can be well
reproduced from theoretical computations Karpenko and Becattini (2017); Xie et
al. (2016, 2017); Li et al. (2017); Shi et al. (2019); Xia et al. (2018); Wei
et al. (2019); Vitiuk et al. (2020); Ivanov et al. (2019); Fu et al. (2021a);
Guo et al. (2021); Li et al. (2021a); Deng et al. (2021). However, when
applying the same formula to the calculation of local spin polarization (i.e.,
the differential spin polarization as a function of the azimuthal angle), the
results deviate substantially from the data Becattini and Karpenko (2018);
Florkowski et al. (2019a); Wu et al. (2019); Xie et al. (2020); Xia et al.
(2019); Becattini et al. (2019a); Liu et al. (2020b); Adam et al. (2019);
Acharya et al. (2021). This is a big puzzle that challenges the thermal-
vorticity interpretation of spin polarization and suggests that the global
equilibrium of spin degree of freedom may not be reached in real heavy-ion
collisions; see more discussions in reviews Liu and Huang (2020); Becattini
and Lisa (2020); Gao et al. (2020); Huang et al. . Possible resolutions have
been proposed including admitting different spin potentials (see the
definition below) to govern the amount of spin polarization Florkowski et al.
(2019a); Wu et al. (2019), assuming strong chiral anomaly effects Liu et al.
(2020b), and so on. New theoretical tools like spin kinetic theory Gao and
Liang (2019); Weickgenannt et al. (2019); Hattori et al. (2019a); Wang et al.
(2019); Liu et al. (2020a); Sheng et al. (2020); Guo (2020); Weickgenannt et
al. (2021a); Yang et al. (2020); Manuel and Torres-Rincon (2021); Sheng et al.
(2021); Weickgenannt et al. (2021b); Lin (2021) and spin hydrodynamics
Florkowski et al. (2018); Hattori et al. (2019b); Montenegro and Torrieri
(2019); Florkowski et al. (2019b); Bhadury et al. (2021a, b); Garbiso and
Kaminski (2020); Gallegos and Gürsoy (2020); Shi et al. (2021); Fukushima and
Pu (2021); Li et al. (2021b); She et al. (2021); Gallegos et al. (2021); Hongo
et al. (2021); Wang et al. (2021); Peng et al. (2021) that can be used to
study the spin polarization out of global equilibrium are also under
development though no numerical results have been reported so far.

Recently, a generalization of the spin polarization formula to local
equilibrium 111We note that the term “local equilibrium” has been given
different meanings in different context. In this article, we follow the
convention of Refs. Zubarev et al. (1979); van Weert (1982); Hongo (2017);
Becattini et al. (2019b, c); Speranza and Weickgenannt (2021) to use “local
equilibrium” to name the Gibbs state described by a density operator that
maximizes the von Neumann entropy under the constraint that the ensemble-
averaged densities of conserved charges (energy, momentum, angular momentum,
and baryon number, etc.) are fixed to the physical values. This is different
from, e.g., Ref. Hongo et al. (2021), where “local equilibrium” is used to
name the local state in which the strict hydrodynamic description can apply.
has been proposed in which not only the thermal vorticity but also the thermal
shear tensor contributes to the mean spin vector Becattini et al. (2021a); Liu
and Yin (2021). It was found that the thermal-shear contribution could be an
important component for local spin polarization under certain conditions: For
example, for hyperons that are isothermally frozen out Becattini et al.
(2021b), or for hyperons whose spin polarization is dominantly determined by
its strange quark constituents Fu et al. (2021b), the calculated spin
polarization can match the experimental data. But it was also argued that the
results could be sensitive to the equation of state of the matter Yi et al.
(2021) and to the procedure of how the thermal vorticity and thermal shear
contributions are combined Florkowski et al. (2021).

It is known that the description of local equilibrium state is not unique; it
depends on the choices of the energy-momentum tensor and spin current in a
compatible way Becattini et al. (2019c); Speranza and Weickgenannt (2021);
Tinti and Florkowski (2020). For example, when choosing the Belinfante
symmetric energy-momentum tensor for Dirac fermions, the corresponding spin
current vanishes, and the local equilibrium is specified by local thermal flow
vector $\beta^{\mu}$. On the other hand, if the canonical form for the energy-
momentum tensor and spin current are chosen, the local equilibrium state
admits a new, independent, rank-two tensor $\omega^{\mu\nu}$ (called spin
potential, which is conjugate to local canonical spin density) in addition to
$\beta^{\mu}$ 222Note that, generically, the $\beta^{\mu}$ in this case does
not coincide with the $\beta^{\mu}$ in the Belinfante case.. The previous
studies did not consider the contribution of $\omega^{\mu\nu}$ to the spin
polarization formula though the canonical spin current was adopted to describe
the local equilibrium state Becattini et al. (2021a); Liu and Yin (2021);
Becattini et al. (2021b); Fu et al. (2021b); Yi et al. (2021). In this
article, we will re-visit the spin polarization at local equilibrium and
derive a new spin polarization formula by using the canonical forms of energy-
momentum and spin current tensors. We will show that the spin potential makes
an additional contribution in the new spin polarization formula and we will
discuss the physical content of the new formula.

We will take metric $\eta_{\mu\nu}={\rm{diag}}(1,-1,-1,-1)$ and the natural
units $k_{B}=c=1$. The Levi-Civita tensor $\epsilon^{\mu\nu\rho\sigma}$ is
normalized as $\epsilon^{0123}=-\epsilon_{0123}=1$.

##  II Derivation using kinetic theory

We recall that for spin-1/2 Dirac fermions the canonical spin current operator
is given by (we restore $\hbar$ in this section)

$\displaystyle\begin{split}\hat{S}^{\mu,\rho\sigma}(x)&\equiv\frac{\hbar}{4}{\bar{\hat{\psi}}}(x)\left\\{\gamma^{\mu},\sigma^{\rho\sigma}\right\\}\hat{\psi}(x)\\\
&=\int\frac{d^{4}p}{(2\pi)^{4}}\hat{S}^{\mu,\rho\sigma}(x,p)\,,\end{split}$
(4)

where $\sigma^{\mu\nu}=i[\gamma^{\mu},\gamma^{\nu}]/2$ and
$\hat{S}^{\mu,\rho\sigma}(x,p)$ is the canonical spin current operator in
phase space which is related to Wigner operator as follows,

$\displaystyle\hat{S}^{\mu,\rho\sigma}(x,p)$ $\displaystyle=$
$\displaystyle\frac{\hbar}{4}{\rm Tr}_{\rm
D}\left[\left\\{\gamma^{\mu},\sigma^{\rho\sigma}\right\\}\hat{W}(x,p)\right]\,,$
(5) $\displaystyle\hat{W}(x,p)$ $\displaystyle=$ $\displaystyle\int
d^{4}se^{-ip\cdot
s/\hbar}\bar{\hat{\psi}}\left(x+\frac{s}{2}\right)\otimes\hat{\psi}\left(x-\frac{s}{2}\right),\;\;\;\;$
(6)

where ${\rm Tr}_{\rm D}$ is trace over Dirac space and
$[\bar{\hat{\psi}}\otimes\hat{\psi}]_{ab}\equiv\bar{\hat{\psi}}_{b}\hat{\psi}_{a}$
with $a,b$ spinor indices. The Wigner function is defined as the ensemble
average of the Wigner operator: $W(x,p)=\langle\hat{W}(x,p)\rangle={\rm
Tr}\left[\hat{\rho}\hat{W}(x,p)\right]$ with $\hat{\rho}$ the density operator
and ${\rm Tr}$ the trace over a complete set of microstates of the system. For
free Dirac fermions at local equilibrium, $\hat{\rho}=\hat{\rho}_{\rm LE}$,
the Wigner function can be directly calculated order by order in $\hbar$ using
second quantization, as we will see in next section. In this section, however,
we will use the results from chiral kinetic theory and spin kinetic theory
which provide an alternative way to calculate the Wigner function. First, we
decompose the Wigner function on the basis of Clifford algebra:

$W=\frac{1}{4}[\mathcal{F}+i\gamma^{5}\mathcal{P}+\gamma^{\mu}\mathcal{V}_{\mu}+\gamma^{5}\gamma^{\mu}\mathcal{A}_{\mu}+\frac{1}{2}\sigma^{\mu\nu}\mathcal{S}_{\mu\nu}]$
(7)

with all the Clifford coefficients being real. Further, one can find that

$S^{\mu,\rho\sigma}(x,p)=-\frac{\hbar}{2}\epsilon^{\mu\rho\sigma\lambda}\mathcal{A}_{\lambda}(x,p)\,,\\\
$ (8)

where $S^{\mu,\rho\sigma}(x,p)=\langle\hat{S}^{\mu,\rho\sigma}(x,p)\rangle$.
The above relation identifies the canonical spin vector $S^{\mu}$ with the
axial current $\mathcal{A}^{\mu}$ in phase space:

$S^{\mu}(x,p)=-\frac{1}{6}\epsilon^{\mu\nu\rho\sigma}S_{\nu,\rho\sigma}(x,p)=\frac{\hbar}{2}\mathcal{A}^{\mu}(x,p).$
(9)

The dynamic equation of Wigner function reads (see, e.g., Refs. Liu et al.
(2020a); Gao and Liang (2019); Weickgenannt et al. (2019); Hattori et al.
(2019a); Wang et al. (2019); Sheng et al. (2020))

$\displaystyle\left[\gamma^{\mu}\left(\frac{i\hbar}{2}\partial_{\mu}+p_{\mu}\right)-m\right]W=0\,,$
(10)

from which one can solve for $\mathcal{A}_{\mu}$ and other Clifford
coefficients order by order in $\hbar$.

###  II.1 Massless case

Let us first consider the massless limit, $m=0$. The solution for
$\mathcal{A}_{\mu}$ up to $O(\hbar)$ is (see, e.g., Refs. Hidaka et al.
(2017); Huang et al. (2018); Liu et al. (2019); Gao et al. (2018) for more
details)

$\begin{split}\mathcal{A}_{\mu}&=4\pi\delta(p^{2})\left(p_{\mu}f_{5}+\hbar\Sigma^{n}_{\mu\nu}\partial^{\nu}f\right),\end{split}$
(11)

where $f_{5}=f_{R}-f_{L}$ and $f=f_{R}+f_{L}$ are the axial and vector
distribution functions with $f_{R}$ and $f_{L}$ the distribution functions of
right-handed and left-handed fermions, respectively,
$\Sigma_{n}^{\mu\nu}\equiv\epsilon^{\mu\nu\rho\sigma}p_{\rho}n_{\sigma}/(2p\cdot
n)$ represents the single-particle spin tensor with $n^{\mu}$ a normalized
timelike frame-fixing vector, $n^{2}=1$. Similarly, the vector current reads

$\begin{split}\mathcal{V}_{\mu}=4\pi\delta(p^{2})\left(p_{\mu}f+\hbar\Sigma^{n}_{\mu\nu}\partial^{\nu}f_{5}\right).\end{split}$
(12)

The concrete forms of the distribution functions $f$ and $f_{5}$ cannot be
determined without knowing the collision terms. However, for collisional
processes that conserve energy, momentum, particle numbers, and angular
momentum, the local-equilibrium distributions should depend only on the linear
combination of these quantities. Therefore, we assume the distribution
functions of right-handed and left-handed fermions to be of Fermi-Dirac form
and read

$\displaystyle f_{R/L}$ $\displaystyle=$
$\displaystyle\theta(p_{0})n_{F}(g_{R/L})+\theta(-p_{0})[1-n_{F}(g_{R/L})]\,,$
(13)

where
$g_{R/L}=\alpha+p\cdot\beta\mp\frac{\hbar}{2}\Sigma_{n}^{\mu\nu}\omega_{\mu\nu}$,
$n_{F}(x)=(e^{x}+1)^{-1}$ is the Fermi-Dirac function, and the Lagrangian
multipliers $\beta_{\mu},\alpha$, $\omega_{\mu\nu}$ have the following
physical meanings: $\omega_{\mu\nu}=-\omega_{\nu\mu}$ is the spin potential
tensor 333Note that $\omega^{\mu\nu}$ is transverse to $n^{\mu}$ so that it
contains only three independent components in accordance with the three
components of spin vector., $\beta_{\mu}$ is the thermal flow velocity, and
$\alpha=-\beta\mu$ with $\mu$ the chemical potential. We will regard
$\beta^{\mu},\alpha$, and $\omega_{\mu\nu}$ as $O(1)$ in $\hbar$ power
counting. Substituting $f_{R/L}$ into $\mathcal{A}_{\mu}$ and keeping terms up
to $O(\hbar)$, we obtain

$\begin{split}\mathcal{A}_{\mu}&=-8\pi\delta(p^{2})\bigg{\\{}\frac{\hbar}{4}\epsilon_{\mu\nu\alpha\beta}p^{\nu}\omega^{\alpha\beta}+\hbar\Sigma^{n}_{\mu\nu}\left[p_{\lambda}(\xi^{\nu\lambda}+\Delta\omega^{\nu\lambda})+\partial^{\nu}\alpha\right]\bigg{\\}}\\\
&\quad\times[\theta(p_{0})-\theta(-p_{0})]n_{F}\left(1-n_{F}\right),\end{split}$
(14)

where $\xi^{\nu\lambda}=\partial^{(\nu}\beta^{\lambda)}$ is the thermal shear
tensor, $\Delta\omega^{\nu\lambda}=\omega^{\nu\lambda}-\varpi^{\nu\lambda}$ is
the deviation of spin potential from thermal vorticity, and we have used the
abbreviation $n_{F}=n_{F}(\alpha+p\cdot\beta)$. The first term is the spin
polarization induced by spin potential, which is not necessarily equal to
thermal vorticity at local equilibrium. The second term is the spin
polarization induced by the thermal shear tensor and the deviation of spin
potential from the thermal vorticity. The last term is due to the spatial
inhomogeneity of $\alpha$.

###  II.2 Massive case

Next, we consider massive fermions. Unlike the massless fermion whose spin is
slaved to the momentum, the spin of a massive fermion is an independent degree
of freedom. This makes the kinetic description of the spin dynamics of massive
fermion more complicated: Comparing to the massless case where two kinetic
equations are enough, we need four kinetic equations (one for the particle
number and three for the spin vector) in the massive case Gao and Liang
(2019); Weickgenannt et al. (2019); Hattori et al. (2019a); Wang et al.
(2019); Liu et al. (2020a); Sheng et al. (2020). Instead of presenting all
these kinetic equations, we here list a few relations for the Clifford
coefficients in Eq. (7) which can be obtained through the $\hbar$ expansion of
Eq. (10) and will be enough for our purpose:

$\displaystyle\begin{split}\mathcal{F}&=4\pi
m\delta(p^{2}-m^{2})f+O(\hbar^{2}),\\\ \mathcal{S}^{\mu\nu}&=4\pi
m\delta(p^{2}-m^{2}){\mathcal{M}}^{\mu\nu}+O(\hbar^{2}),\\\
\hbar\partial_{\mu}\mathcal{F}&=2p^{\nu}\mathcal{S}_{\mu\nu}+O(\hbar^{2}),\\\
\mathcal{A}_{\mu}&=-\frac{1}{2m}\epsilon_{\mu\nu\rho\sigma}p^{\nu}{\mathcal{S}}^{\rho\sigma}+O(\hbar^{2}),\\\
\mathcal{V}_{\mu}&=4\pi\delta(p^{2}-m^{2})p_{\mu}f+\frac{\hbar}{2m}\partial^{\nu}\mathcal{S}_{\mu\nu}+O(\hbar^{2}),\end{split}$
(15)

where $\mathcal{F}$ and $\mathcal{S}^{\mu\nu}$ are parameterized by $f$ and
$\mathcal{M}^{\mu\nu}$ which are identified as the on-shell particle number
and tensorial spin densities (distribution functions) in phase space. We count
$f$ as $O(1)$ and $\mathcal{M}^{\mu\nu}$ as $O(\hbar)$ quantities. Multiplying
the second equation in (15) by $p\cdot n$, then substituting
$\mathcal{S}^{\mu\nu}$ into the forth equation, and using the first and third
equations and the Schouten identity
$p_{\lambda}\epsilon_{\mu\nu\rho\sigma}+p_{\mu}\epsilon_{\nu\rho\sigma\lambda}+p_{\nu}\epsilon_{\rho\sigma\lambda\mu}+p_{\rho}\epsilon_{\sigma\lambda\mu\nu}+p_{\sigma}\epsilon_{\lambda\mu\nu\rho}=0$,
we express $\mathcal{A}^{\mu}$ in terms of $f$ and $\mathcal{M}^{\mu\nu}$, up
to $O(\hbar)$ 444This equation has been used in Ref. Sheng et al. (2020) to
show that in the massless limit, $m\rightarrow 0$, $\mathcal{A}^{\mu}$
smoothly goes into its massless counterpart in Eq. (11) by recognizing
$\Sigma^{n}_{\rho\sigma}{\mathcal{M}}^{\rho\sigma}$ as the axial distribution
function $f_{5}$.,

$\begin{split}\mathcal{A}_{\mu}&=4\pi\delta(p^{2}-m^{2})\bigg{(}p_{\mu}\Sigma^{n}_{\rho\sigma}{\mathcal{M}}^{\rho\sigma}+\hbar\Sigma^{n}_{\mu\nu}\partial^{\nu}f-\frac{m^{2}}{2p\cdot
n}\epsilon_{\mu\nu\rho\sigma}n^{\nu}{\mathcal{M}}^{\rho\sigma}\bigg{)}.\end{split}$
(16)

This expression shows that only the $n^{\mu}$-transverse components of
$\mathcal{M}^{\mu\nu}$ (corresponding to the three components of spin vector)
contribute to $\mathcal{A}^{\mu}$ at $O(\hbar)$ order. In fact, by decomposing
$\mathcal{M}^{\mu\nu}=\mathcal{E}^{\mu}n^{\nu}-\mathcal{E}^{\nu}n^{\mu}+\epsilon^{\mu\nu\alpha\beta}n_{\alpha}\mathcal{M}_{\beta}$,
it can be shown that $\mathcal{E}^{\mu}$ is determined by $f$ and
$\mathcal{M}^{\mu}$ Sheng et al. (2020). Using once again the Schouten
identity for the first term in Eq. (16), one can re-write $\mathcal{A}^{\mu}$
as

$\begin{split}\mathcal{A}_{\mu}&=-8\pi\delta(p^{2}-m^{2})\bigg{(}\frac{1}{4}\epsilon_{\mu\nu\alpha\beta}p^{\nu}\mathcal{M}^{\alpha\beta}_{\perp}-\frac{\hbar}{2}\Sigma^{n}_{\mu\nu}\partial^{\nu}f+\Sigma^{n}_{\mu\nu}\mathcal{M}_{\perp}^{\nu\lambda}p_{\lambda}\bigg{)},\end{split}$
(17)

where
$\mathcal{M}^{\mu\nu}_{\perp}=\epsilon^{\mu\nu\alpha\beta}n_{\alpha}\mathcal{M}_{\beta}$.
At local equilibrium, for an arbitrary space-like direction $r^{\mu}$ (with
$r^{2}=-1$) regarded as the spin quantization axis, both $f$ and
$\mathcal{M}_{\perp}^{\mu\nu}r_{\mu\nu}$ (the tensorial spin density projected
onto $r^{\mu}$) with
$r_{\mu\nu}=\epsilon_{\mu\nu\rho\sigma}n^{\rho}r^{\sigma}/2$ should be
determined by $\alpha,\beta$, and $\omega^{\mu\nu}$ via $f=f_{+}+f_{-}$ and
$\mathcal{M}_{\perp}^{\mu\nu}r_{\mu\nu}=f_{+}-f_{-}$ where

$\displaystyle f_{\pm}$ $\displaystyle=$
$\displaystyle\theta(p_{0})n_{F}(g_{\pm})+\theta(-p_{0})[1-n_{F}(g_{\pm})]\;,$
(18)

with
$g_{\pm}=\alpha+p\cdot\beta\mp\frac{\hbar}{2}\omega^{\alpha\beta}r_{\alpha\beta}$.
We have used the same assumption as that for Eq. (13) to obtain the above
equation. One thus finds that
$f=2[\theta(p_{0})n_{F}+\theta(-p_{0})(1-n_{F})]+O(\hbar^{2})$ with
$n_{F}=n_{F}(\alpha+p\cdot\beta)$ and
$\mathcal{M}_{\perp}^{\mu\nu}=-\hbar\omega^{\mu\nu}[\theta(p_{0})-\theta(-p_{0})]n_{F}^{\prime}(\alpha+p\cdot\beta)+O(\hbar^{2})$.
Substituting $f$ and $\mathcal{M}_{\perp}^{\mu\nu}$ into Eq. (17), we obtain

$\begin{split}\mathcal{A}_{\mu}&=-8\pi\hbar\delta(p^{2}-m^{2})\bigg{\\{}\frac{1}{4}\epsilon_{\mu\nu\alpha\beta}p^{\nu}\omega^{\alpha\beta}+\Sigma^{n}_{\mu\nu}\left[p_{\lambda}(\xi^{\nu\lambda}+\Delta\omega^{\nu\lambda})+\partial^{\nu}\alpha\right]\bigg{\\}}\\\
&\quad\times[\theta(p_{0})-\theta(-p_{0})]n_{F}\left(1-n_{F}\right).\end{split}$
(19)

Note that Eq. (19) is almost identical with its massless counterpart, Eq.
(14), except for the different on-shell condition.

##  III Derivation using local equilibrium density operator

In this section, we present another derivation of Eq. (19) by using the local
equilibrium density operator Zubarev et al. (1979); van Weert (1982). We will
focus on massive Dirac fermions and we will follow closely Refs. Becattini et
al. (2021a); Becattini . We set $\hbar=1$ in this section. The starting point
is the following density operator describing a local equilibrium state
specified by $\alpha(x),\beta_{\mu}(x)$, and $\omega_{\mu\nu}(x)$,

$\begin{split}\hat{\rho}_{\rm LE}&=\frac{1}{Z_{\rm LE}}\exp\bigg{\\{}-\int
d\Xi_{\mu}(y)\big{[}\hat{T}^{\mu\nu}(y)\beta_{\nu}(y)+\hat{J}^{\mu}(y)\alpha(y)-\frac{1}{2}\hat{S}^{\mu,\rho\sigma}(y)\omega_{\rho\sigma}(y)\big{]}\bigg{\\}},\end{split}$
(20)

where $\hat{T}^{\mu\nu}$, $\hat{J}^{\mu}$, and $\hat{S}^{\mu,\rho\sigma}$ are
the canonical operators for energy-momentum tensor, vector current, and spin
current, respectively, and $Z_{\rm LE}$ is the partition function. For a local
Heisenberg operator $\hat{O}(x)$, we have its ensemble average given by

$\displaystyle O(x)$ $\displaystyle\equiv$
$\displaystyle\langle\hat{O}(x)\rangle\equiv{\rm Tr}[\hat{\rho}_{\rm
LE}\hat{O}(x)]=\frac{1}{Z_{\rm LE}}{\rm
Tr}\left[e^{\hat{A}+\hat{B}}\hat{O}(x)\right],\qquad$ (21)

where we introduced the abbreviations

$\displaystyle\hat{A}$ $\displaystyle=$
$\displaystyle-\hat{P}^{\mu}\beta_{\mu}(x)-\hat{N}\alpha(x),$ (22)
$\displaystyle\hat{B}$ $\displaystyle=$ $\displaystyle-\int
d\Xi_{\mu}(y)\big{[}\hat{T}^{\mu\nu}(y)\Delta\beta_{\nu}(y)+\hat{J}^{\mu}(y)\Delta\alpha(y)-\frac{1}{2}\hat{S}^{\mu,\rho\sigma}(y)\omega_{\rho\sigma}(y)\big{]},$
(23)

with $\hat{P}^{\mu}=\int d\Xi_{\nu}(y)\hat{T}^{\nu\mu}(y)$, $\hat{N}=\int
d\Xi_{\nu}(y)\hat{J}^{\nu}(y)$,
$\Delta\beta_{\mu}(y)=\beta_{\mu}(y)-\beta_{\mu}(x)$, and
$\Delta\alpha(y)=\alpha(y)-\alpha(x)$.

In the following, we assume that $\omega_{\rho\sigma}(y)$ is small on
hypersurface $\Xi$ and $\Delta\beta_{\mu}(y)$ and $\Delta\alpha(y)$ are also
small on $\Xi$. This is justified when $\beta_{\mu}$ and $\alpha$ vary slowly
within the correlation lengths between $\hat{O}$ and the operators
$\hat{T}^{\mu\nu}$ and $\hat{J}^{\mu}$. Thus Eq. (21) can be calculated in a
perturbative way in derivative expansion using the power counting rule
$\Delta\beta_{\nu}(y)=(y-x)^{\lambda}\partial_{\lambda}\beta_{\nu}(x)+\cdots=O(\partial)$,
$\Delta\alpha(y)=(y-x)^{\lambda}\partial_{\lambda}\alpha(x)+\cdots=O(\partial)$,
and $\omega_{\rho\sigma}(y)=O(\partial)$. It should be emphasized that in the
current content the gradient expansion is equivalent to the $\hbar$ expansion
that we adopted in last section. Adopting
$e^{\hat{A}+\hat{B}}=e^{\hat{A}}+e^{\hat{A}}\int^{1}_{0}e^{-\lambda\hat{A}}\hat{B}e^{\lambda\hat{A}}+\cdot\cdot\cdot$,
we expand the ensemble average $O(x)$ into the following form

$\displaystyle O(x)$ $\displaystyle=$ $\displaystyle
O_{0}(x)+O_{1}(x)+\cdots\,,$ (24)

where

$\displaystyle O_{0}(x)=\langle\hat{O}(x)\rangle_{0}\equiv\frac{1}{Z_{0}}{\rm
Tr}\left(e^{\hat{A}}\hat{O}(x)\right),\qquad$ (25) $\displaystyle
O_{1}(x)\equiv\langle\hat{O}(x)\rangle_{(T)}+\langle\hat{O}(x)\rangle_{(J)}+\langle\hat{O}(x)\rangle_{(S)}\,\,,$
(26)

with

$\begin{split}\langle\hat{O}(x)\rangle_{(T)}&\equiv-\int_{0}^{1}d\lambda\int
d\Xi_{\mu}(y)\Delta\beta_{\nu}(y)\langle\hat{T}^{\mu\nu}(y-i\lambda\beta(x))\hat{O}(x)\rangle_{0,c}\,,\quad\\\
\langle\hat{O}(x)\rangle_{(S)}&\equiv\frac{1}{2}\int_{0}^{1}d\lambda\int
d\Xi_{\mu}(y)\omega_{\rho\sigma}(y)\langle\hat{S}^{\mu,\rho\sigma}(y-i\lambda\beta(x))\hat{O}(x)\rangle_{0,c}\,,\qquad\\\
\langle\hat{O}(x)\rangle_{(J)}&\equiv-\int_{0}^{1}d\lambda\int
d\Xi_{\mu}(y)\Delta\alpha(y)\langle\hat{J}^{\mu}(y-i\lambda\beta(x))\hat{O}(x)\rangle_{0,c}\,,\end{split}$
(27)

and $Z_{0}={\rm Tr}e^{\hat{A}}$. In the above expressions,
$\langle\cdots\rangle_{0,c}$ means the connected part of the correlator.

We now apply the above procedure to the Wigner function for massive Dirac
fermions:

$W(x,p)=W_{0}(x,p)+W_{1}(x,p)+\cdots.$ (28)

We first express the Wigner operator in Eq. (6) using the free field operator

$\begin{split}\hat{\psi}(x)&=\sum^{2}_{\sigma=1}\frac{1}{(2\pi)^{3/2}}\int\frac{d^{3}{\bm{k}}}{2E_{k}}\big{[}u_{\sigma}({\bm{k}})e^{-ik\cdot
x}\hat{a}_{\sigma}({\bm{k}})+v_{\sigma}({\bm{k}})e^{ik\cdot
x}\hat{b}^{\dagger}_{\sigma}({\bm{k}})\big{]}\,,\end{split}$ (29)

where $E_{k}=\sqrt{{\bm{k}}^{2}+m^{2}}$, the creation and annihilation
operators for particles and antiparticles satisfy the anti-commutation
relation
$\\{\hat{a}_{\sigma}({\bm{k}}),\hat{a}^{\dagger}_{\sigma^{\prime}}({\bm{q}})\\}=\\{\hat{b}_{\sigma}({\bm{k}}),\hat{b}^{\dagger}_{\sigma^{\prime}}({\bm{q}})\\}=2E_{k}\delta_{\sigma\sigma^{\prime}}\delta^{3}({\bm{k}}-\bm{q})$,
and the spinors satisfy the relation
$\bar{u}_{\sigma}({\bm{k}})u_{\sigma^{\prime}}({\bm{k}})=2m\delta_{\sigma\sigma^{\prime}}=-\bar{v}_{\sigma}({\bm{k}})v_{\sigma^{\prime}}({\bm{k}})$.
Substituting the above equation into Eq. (6), we arrive at the following form
for the Wigner operator (Here and below, we focus on particle branch only):

$\begin{split}\hat{W}(x,p)&=2\pi\sum^{2}_{\sigma,\sigma^{\prime}=1}\int\frac{d^{3}{\bm{q}}}{2E_{q}}\int\frac{d^{3}{\bm{k}}}{2E_{k}}e^{i(k-q)\cdot
x}\delta^{4}\left(p-\frac{k}{2}-\frac{q}{2}\right)\bar{u}_{\sigma}({\bm{k}})\otimes
u_{\sigma^{\prime}}({\bm{q}})\hat{a}^{\dagger}_{\sigma}({\bm{k}})\hat{a}_{\sigma^{\prime}}({\bm{q}}).\end{split}$
(30)

With the above preparations, we now calculate the Wigner function. Using
$\langle\hat{a}^{\dagger}_{\sigma}({\bm{k}})\hat{a}_{\sigma^{\prime}}({\bm{q}})\rangle_{0}=2E_{k}\delta_{\sigma\sigma^{\prime}}\delta^{3}({\bm{k}}-{\bm{q}})n_{F}(\alpha+p\cdot\beta)$,
one find that the zeroth order Wigner function takes the form
$W_{0}(x,p)=2\pi(p\\!\\!\\!/+m)\theta(p_{0})\delta(p^{2}-m^{2})n_{F}(\alpha+p\cdot\beta)$,
which is spin independent: ${\rm
Tr}_{D}\left[\gamma^{\mu}\gamma_{5}W_{0}(x,p)\right]=0$. This gives the
leading-order vector current as

$\displaystyle\mathcal{V_{\mu}}=8\pi\delta(p^{2}-m^{2})\theta(p_{0})p_{\mu}n_{F}(\alpha+p\cdot\beta),$
(31)

which coincides with the last expression in Eq. (15) at local equilibrium with
particle branch picked up. The first order Wigner function can be derived by
using Eq. (27), with the help of
$\hat{T}^{\mu\nu}(x)=\int\frac{d^{4}p}{(2\pi)^{4}}p^{\nu}{\rm
Tr}_{D}\left[\gamma^{\mu}\hat{W}(x,p)\right]$,
$\hat{J}^{\mu}(x)=\int\frac{d^{4}p}{(2\pi)^{4}}{\rm
Tr}_{D}\left[\gamma^{\mu}\hat{W}(x,p)\right]$, and Eq. (5). Substituting the
Wigner operator in Eq. (30) into the right-hand side of Eq. (27), we obtain

$\displaystyle\langle\hat{W}(x,p)\rangle_{(T/S/J)}$ $\displaystyle=$
$\displaystyle(2\pi)\int_{0}^{1}d\lambda\int\frac{d^{3}{\bm{k}}}{2E_{k}}\int\frac{d^{3}{\bm{q}}}{2E_{q}}\delta^{4}\left(p-\frac{q+k}{2}\right)(\gamma\cdot
k+m)\hat{t}_{\mu}I_{(T/S/J)}^{\mu}(\gamma\cdot q+m)$ (32) $\displaystyle\times
e^{\lambda(k-q)\cdot\beta(x)}n_{F}(k)\left[1-n_{F}(q)\right]\,,$

where $\hat{t}^{\mu}$ is the normal vector of the hypersurface $\Xi$,
$n_{F}(p)=n_{F}(\alpha+\beta(x)\cdot p)$, and

$\begin{split}I_{(T)}^{\mu}&\equiv-\gamma^{\mu}p^{\nu}\left[\partial_{\lambda}\beta_{\nu}(x)\right]\Delta^{\lambda}_{\beta}\left[i\partial_{q}^{\beta}\delta^{3}(\bm{q}-\bm{k})\right]\,,\\\
I_{(S)}^{\mu}&\equiv\frac{1}{4}\epsilon^{\mu\nu\rho\sigma}\gamma^{5}\gamma_{\nu}\omega_{\rho\sigma}\delta^{3}(\bm{q}-\bm{k})\,,\\\
I_{(J)}^{\mu}&\equiv\gamma^{\mu}\left[\partial_{\lambda}\alpha(x)\right]\Delta^{\lambda}_{\beta}\left[i\partial_{q}^{\beta}\delta^{3}(\bm{q}-\bm{k})\right]\,,\end{split}$
(33)

with $\Delta^{\mu\nu}=\eta^{\mu\nu}-\hat{t}^{\mu}\hat{t}^{\nu}$. In order to
derive Eq. (32), we have assumed that $\hat{t}_{\mu}$ does not deviate
strongly from the coordinate time direction so that approximately

$\begin{split}\int
d\Xi_{\mu}(y)(y-x)^{\alpha}e^{-i(p-q)\cdot(y-x)}\approx(2\pi)^{3}\hat{t}_{\mu}\Delta^{\alpha}_{\beta}\frac{i\partial}{\partial
p_{\beta}}\delta^{3}({\bm{p}}-{\bm{q}})+(2\pi)^{3}\hat{t}\cdot(y-x)\hat{t}^{\alpha}\hat{t}_{\mu}\delta^{3}({\bm{p}}-{\bm{q}}).\end{split}$
(34)

In heavy-ion collisions, the true freeze-out hypersurface $\Xi$ may not
satisfy the above condition and in this case one can choose $\hat{t}^{\mu}$ as
the time direction of an earlier-stage fluid occupying a flat three-
dimensional hypersurface $\Xi_{B}$, if the integration over the four-volume
encompassed by $\Xi$ and $\Xi_{B}$ does not contribute to the spin vector
Becattini et al. (2021a).

With the first order solution of Wigner function in Eq. (32), we directly
derive the following axial vector following similar procedures as Ref.
Becattini et al. (2021a):

$\begin{split}\mathcal{A}_{\mu}&=-8\pi\delta(p^{2}-m^{2})\theta(p_{0})n_{F}(1-n_{F})\bigg{\\{}\frac{1}{4}\epsilon_{\mu\nu\rho\sigma}p^{\nu}\omega^{\rho\sigma}+\Sigma^{\hat{t}}_{\mu\nu}\left[(\xi^{\nu\lambda}+\Delta\omega^{\nu\lambda})p_{\lambda}+\partial^{\nu}\alpha\right]\bigg{\\}},\end{split}$
(35)

which is exactly the particle branch of Eq. (19) with $n^{\mu}=\hat{t}^{\mu}$.
The anti-particle branch can be similarly obtained by replacing Eq. (30) with
the anti-particle Wigner operator.

##  IV A spin Cooper-Frye formula

With having $\mathcal{A}^{\mu}(x,p)$, the spin vector $S^{\mu}(x,p)$ is given
by Eq. (9). Let us focus on the spin polarization of particles; the anti-
particle spin polarization can be similarly discussed. Integrating Eq. (14) or
Eq. (19) [or Eq. (35)] over $p_{0}$ from $0$ to $\infty$ we obtain the on-
shell spin vector (we set $\hbar=1$):

$\displaystyle S^{\mu}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\infty}\frac{dp_{0}}{2\pi}\mathcal{A}^{\mu}(x,p)$
(36) $\displaystyle=$
$\displaystyle-\frac{1}{4E_{p}}\left\\{\epsilon^{\mu\nu\rho\sigma}p_{\nu}\omega_{\rho\sigma}+4\Sigma_{n}^{\mu\nu}\left[(\xi_{\nu\lambda}+\Delta\omega_{\nu\lambda})p^{\lambda}+\partial_{\nu}\alpha\right]\right\\}n_{F}(1-n_{F}),$

where, on the right-hand side, $p_{0}=E_{p}=\sqrt{{\bm{p}}^{2}+m^{2}}$ for
both $m>0$ and $m=0$ cases. In order to more clearly understand the different
terms in Eq. (36), it is instructive to choose $n^{\mu}$ as the coordinate
time direction, i.e., $n^{\mu}=(1,0,0,0)$, and consider the low-velocity
limit, i.e., $u^{\mu}\approx(1,{\bm{v}})$ with $|{\bm{v}}|\ll 1$. The zeroth
component of the spin vector is (half) the helicity induced by finite spin
potential:

$\displaystyle
S^{0}(x,{\bm{p}})\approx\frac{1}{2}\left[n_{F}\left(\alpha+\beta
E_{p}-\frac{\bm{\omega}\cdot{\bm{p}}}{2E_{p}}\right)-n_{F}\left(\alpha+\beta
E_{p}+\frac{\bm{\omega}\cdot{\bm{p}}}{2E_{p}}\right)\right],$ (37)

where $\omega^{i}=\epsilon^{ijk}\omega_{jk}/2$ is the spin potential vector.
The spatial components are
$\bm{S}=\bm{S}_{(\omega)}+\bm{S}_{(\Omega)}+\bm{S}_{(\sigma)}+\bm{S}_{(T)}+\bm{S}_{(\alpha)}$
with

$\displaystyle S^{i}_{(\omega)}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle\left[\frac{\omega^{i}}{2}-\frac{{\bm{p}}^{2}\omega^{i}-\bm{\omega}\cdot{\bm{p}}\,p^{i}}{2E_{p}^{2}}\right]n_{F}(1-n_{F}),$
(38) $\displaystyle S^{i}_{(\Omega)}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle\frac{{\bm{p}}^{2}\Omega^{i}-\bm{\Omega}\cdot{\bm{p}}\,p^{i}}{2T\,E_{p}^{2}}n_{F}(1-n_{F}),$
(39) $\displaystyle S^{i}_{(\sigma)}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle\frac{\epsilon^{ijk}\,p^{j}p^{l}\,\sigma_{kl}}{2T\,E_{p}^{2}}n_{F}(1-n_{F}),$
(40) $\displaystyle S^{i}_{(T)}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle-\frac{\left({\bm{p}}\times\bm{\nabla}T\right)^{i}}{2T^{2}\,E_{p}}\,n_{F}(1-n_{F}),$
(41) $\displaystyle S^{i}_{(\alpha)}(x,{\bm{p}})$ $\displaystyle=$
$\displaystyle\frac{\left({\bm{p}}\times\bm{\nabla}\alpha\right)^{i}}{2E_{p}^{2}}n_{F}(1-n_{F}),$
(42)

where $\bm{\Omega}=(\bm{\nabla}\times{\bm{v}})/2$ is the non-relativistic
vorticity,
$\sigma_{ij}=\left(\partial_{i}v_{j}+\partial_{j}v_{i}+2\delta_{ij}\bm{\nabla}\cdot{\bm{v}}/3\right)/2$
is the non-relativistic shear tensor, $n_{F}=n_{F}(\alpha+\beta E_{p})$, and
we have kept only the leading terms in ${\bm{v}}$.

Equation (36) is particularly useful for the calculation of spin polarization
in transport models because the transport models can provide full phase-space
information about the particles’ production and evolution. But it is less
useful for hydrodynamic modeling of heavy-ion collisions in which the
particles are produced (emitted) on some freeze-out hypersurface $\Xi$. In
this case, we need to average $S^{\mu}(x,{\bm{p}})$ over $\Xi$ to obtain the
mean spin vector 555This definition of the mean spin vector can also be
expressed as
$\bar{S}^{\mu}=-\epsilon^{\mu\nu\rho\sigma}p_{\nu}\bar{S}_{\rho\sigma}/(2E_{p})$
with $\bar{S}_{\rho\sigma}$ the mean spin tensor given by
$\bar{S}^{\mu\nu}=\int_{0}^{\infty}dp_{0}\int
d\Xi_{\lambda}S^{\lambda,\mu\nu}(x,p)/\left[\int_{0}^{\infty}dp_{0}\int
d\Xi_{\lambda}\mathcal{V}^{\lambda}(x,p)\right]$. We note that this is
different from Refs. Becattini et al. (2013, 2021a); Becattini  where the mean
spin vector is defined as
$\bar{S}^{\mu}=-\epsilon^{\mu\nu\rho\sigma}p_{\nu}\bar{S}_{\rho\sigma}/(2m)$
which works only for massive fermions. With this latter definition, finally
$\bar{S}^{\mu}$ is also expressed by Eq. (43) but with $E_{p}$ appearing in
the denominator replaced by $m$.

$\displaystyle\bar{S}^{\mu}({\bm{p}})$ $\displaystyle=$
$\displaystyle\frac{\int d\Xi\cdot p\;S^{\mu}(x,{\bm{p}})}{2\int d\Xi\cdot
p\;n_{F}}$ (43) $\displaystyle=$ $\displaystyle-\frac{\int d\Xi\cdot
p\left\\{\epsilon_{\mu\nu\alpha\beta}p^{\nu}\omega^{\alpha\beta}+4\Sigma^{n}_{\mu\nu}\left[p_{\lambda}(\xi^{\nu\lambda}+\Delta\omega^{\nu\lambda})+\partial^{\nu}\alpha\right]\right\\}n_{F}(1-n_{F})}{8E_{p}\;\int
d\Xi\cdot p\;n_{F}},$

where $p^{\mu}$ on the right-hand side is on-shell. This is a Cooper-Frye type
formula for spin polarization which connects the momentum-space distribution
of mean spin vector of the particles emitted from $\Xi$ with the fluid
properties characterized by $\omega^{\mu\nu}(x)$, $\beta^{\mu}(x)$, and
$\alpha(x)$ on $\Xi$. In the following, we give some comments on formula (43).

(1) At global equilibrium, $\xi^{\mu\nu}=\partial^{\mu}\alpha=0$ and
$\omega^{\mu\nu}=\varpi^{\mu\nu}$ Liu et al. (2019, 2020a); Becattini (2012),
Eq. (43) becomes

$\displaystyle\bar{S}^{\mu}({\bm{p}})$ $\displaystyle=$
$\displaystyle-\frac{\epsilon_{\mu\nu\alpha\beta}\,p^{\nu}\int d\Xi\cdot
p\;\varpi^{\alpha\beta}\;n_{F}(1-n_{F})}{8E_{p}\int d\Xi\cdot p\;n_{F}},$ (44)

which is essentially Eq. (1) when applied to massive fermions.

(2) There is a thermal shear contribution to $\bar{S}^{\mu}({\bm{p}})$ which
may have important implications to phenomenology of spin polarization in
heavy-ion collisions, as having been discussed in Ref. Becattini et al.
(2021a) with $n^{\mu}=\hat{t}^{\mu}$ and in Ref. Liu and Yin (2021) with
$n^{\mu}=u^{\mu}$.

(3) When the system permits finite helicity density at local equilibrium,
there will appear a new contribution to the mean spin vector. This can be
derived, for massless fermions, by relaxing $g_{R/L}$ in Eq. (13) to
$g_{R/L}=\alpha_{R/L}+p\cdot\beta\mp\frac{\hbar}{2}\Sigma_{n}^{\mu\nu}\omega_{\mu\nu}$
where $\alpha_{R}=\alpha+\alpha_{5}$ and $\alpha_{L}=\alpha-\alpha_{5}$ with
$\alpha_{5}=-\beta\mu_{5}$ characterizes the chirality imbalance (which, for
the particle branch of the spectrum, equals to the helicity imbalance).
Counting $\alpha_{5}$ as $O(\hbar)$, we find the leading contribution of
$\alpha_{5}$ to be Liu et al. (2020a)

$\begin{split}\bar{S}_{5}^{\mu}({\bm{p}})=-\frac{p^{\mu}\int d\Xi\cdot
p\;\alpha_{5}\;n_{F}(1-n_{F})}{2E_{p}\int d\Xi\cdot p\;n_{F}},\end{split}$
(45)

where $p^{\mu}$ is on-shell. A similar term can exist for massive fermions as
well, which can be derived by adding a term $-\int
d\Xi_{\mu}(y)J_{5}^{\mu}(y)\alpha_{5}(y)$ to the exponent in Eq. (20) with
$\alpha_{5}(y)$ counted as $O(\partial)$ and following the same procedure in
Sec. III Becattini et al. (2021c). The observable consequences of
$\bar{S}^{\mu}_{5}$ in heavy-ion collisions have been discussed in Refs.
Becattini et al. (2021c); Gao (2021).

(4) The timelike unit vector $n^{\mu}$ does not need to coincide with
$\hat{t}^{\mu}$ (the normal direction of hypersurface $\Xi$) or
$u^{\mu}=T\beta^{\mu}$ (the flow velocity). We can call $n^{\mu}$ the
helicity-frame fixing vector because in its rest frame
$\mathcal{M}^{\mu\nu}\Sigma^{n}_{\mu\nu}$ represents the projection of spin
onto the momentum direction. It is crucial to notice that, because the Wigner
function is Lorentz covariant, the axial current $\mathcal{A}^{\mu}(x,p)$ is
also Lorentz covariant: $\mathcal{A}^{\mu}(\Lambda x,\Lambda
p)=\Lambda^{\mu}_{~{}\nu}\mathcal{A}^{\nu}(x,p)$ for an arbitrary Lorentz
transformation $\Lambda^{\mu}_{~{}\nu}$. This requires $\mathcal{A}^{\mu}$ to
be actually $n^{\mu}$ independent Hidaka et al. (2017); Liu et al. (2020a);
Sheng et al. (2020) and therefore $\omega_{\mu\nu}$ must be $n^{\mu}$
dependent. Under an infinitesimal change $n^{\mu}\rightarrow
n^{\prime\mu}=n^{\mu}+\delta n^{\nu}$, the change in $\omega^{\mu\nu}$,
$\omega_{(n)}^{\mu\nu}\rightarrow\omega_{(n^{\prime})}^{\mu\nu}=\omega_{(n)}^{\mu\nu}+\delta\omega_{(n)}^{\mu\nu}$
is constrained by the following relation 666Note that at global equilibrium
$\omega^{(n)}_{\mu\nu}$ is determined by thermal vorticity and is thus
$n^{\mu}$-independent and Eq. (46) identically holds.:

$\displaystyle\delta\omega^{(n)}_{\mu\nu}\Sigma^{\mu\nu}_{n}=\frac{2}{p\cdot
n}\Sigma_{n}^{\mu\nu}\delta
n_{\nu}\left[\omega^{(n)}_{\mu\lambda}p^{\lambda}+\partial_{\mu}(\alpha+\beta\cdot
p)\right].$ (46)

In particular, for the massive case, we are allowed to choose
$n^{\mu}=p^{\mu}/m$ and the mean spin vector is thus expressed by spin
potential only,

$\displaystyle\bar{S}^{\mu}({\bm{p}})$ $\displaystyle=$
$\displaystyle-\frac{\epsilon_{\mu\nu\alpha\beta}p^{\nu}\int d\Xi\cdot
p\;\omega_{(p)}^{\alpha\beta}\;n_{F}(1-n_{F})}{8E_{p}\;\int d\Xi\cdot
p\;n_{F}}.$ (47)

(5) In heavy-ion collisions, one may be interested in the global spin
polarization $\langle P^{\mu}\rangle=\langle\bar{S}^{\mu}\rangle/s$ ($s=1/2$
for Dirac fermions) where $\langle\bar{S}^{\mu}\rangle$ is the spin vector
averaged over whole phase space,

$\displaystyle\langle\bar{S}^{\mu}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{N}\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}E_{p}}\int
d\Xi\cdot p\;S^{\mu}(x,{\bm{p}}),$ (48)

where $N=2\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}E_{p}}\int d\Xi\cdot p\;n_{F}$ is
the total particle number. Note that $E_{p}$ in the denominator of momentum
integral is necessary to make the dimension of phase space correct. As a case
study, we consider an idealized situation where the normal vector
$\hat{t}^{\mu}$ coincides with fluid velocity $u^{\mu}$ and we obtain

$\displaystyle\langle\bar{S}_{\mu}\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\epsilon_{\mu\nu\alpha\beta}\int d\Xi\cdot
u\,u^{\nu}\left[\omega_{(u)}^{\alpha\beta}F(x)-\frac{2}{3}\Delta\omega_{(u)}^{\alpha\beta}G(x)\right],$
(49)

where $F(x)=N^{-1}\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}}\,n_{F}(1-n_{F})$,
$G(x)=N^{-1}\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}}\frac{{\bm{p}}^{2}}{E^{2}_{p}}\,n_{F}(1-n_{F})$,
and we have taken $n^{\mu}=u^{\mu}$ to perform the calculation and used the
subscript $(u)$ to indicate that $\omega^{\alpha\beta}_{(u)}$ is the spin
potential when $n^{\mu}=u^{\mu}$. Note that, in this case, the gradient of
$\alpha$ and the thermal shear do not contribute to the global spin
polarization.

##  V Conclusion

We have re-visited the spin polarization formula for Dirac fermions at local
thermal equilibrium. Using two different approaches, the chiral/spin kinetic
theory and the method of non-equilibrium density operator, we have derived a
Cooper-Frye type formula that connects the mean spin vector (or equivalently,
spin polarization vector) of Dirac fermions on a hypersurface with the
temperature, chemical potential, flow velocity, and spin potential that
characterize the Gibbs state for local equilibrium. The formula is an
extension of the previously known results valid at global equilibrium, or
under the finite shear tensor, to the situation where a finite independent
spin potential is presented. In such a way, the present formula can apply to
spin hydrodynamics as machinery to convert the spin potential (and
temperature, flow velocity, etc) calculated in spin hydrodynamics to the spin
polarization observable on some freeze-out hypersurface. Application to spin
hydrodynamics will be our future task. Another future task is to extend the
current calculations, which is at $O(\hbar)$ or equivalently at $O(\partial)$
order, to higher orders. Such an extension is straightforward with the
chiral/spin kinetic theory and the method of non-equilibrium density operator
though the calculations would be more involved.

###### Acknowledgements.

We thank Francesco Becattini, Matteo Buzzegoli, Shuai Liu, Shi Pu, Shuzhe Shi,
Qun Wang, Yi Yin for useful discussions. This work is supported by NSFC under
Grant No. 12075061 and Shanghai NSF under Grant No. 20ZR1404100. Y.-C.L is
also supported by China Postdoctoral Science Foundation under Grant No.
2020M681139 and NSFC under Grant No. 12047516. Note added.— During the final
preparation of the present paper, we became aware of a simultaneous work
Buzzegoli (2021) which bears some overlap with our results.

## References

  * Fantoni et al. (2001) S. Fantoni, A. Sarsa, and K. E. Schmidt, Phys. Rev. Lett. 87, 181101 (2001), arXiv:nucl-th/0106026 .
  * Vidana and Bombaci (2002) I. Vidana and I. Bombaci, Phys. Rev. C 66, 045801 (2002), arXiv:nucl-th/0203061 .
  * Isayev and Yang (2004) A. A. Isayev and J. Yang, Phys. Rev. C 69, 025801 (2004), arXiv:nucl-th/0307113 .
  * Stein et al. (2016) M. Stein, A. Sedrakian, X.-G. Huang, and J. W. Clark, Phys. Rev. C 93, 015802 (2016), arXiv:1510.06000 [nucl-th] .
  * Adamczyk et al. (2017) L. Adamczyk et al. (STAR), Nature 548, 62 (2017), arXiv:1701.06657 [nucl-ex] .
  * Adam et al. (2018) J. Adam et al. (STAR), Phys. Rev. C 98, 014910 (2018), arXiv:1805.04400 [nucl-ex] .
  * Abdallah et al. (2021) M. S. Abdallah et al. (STAR), (2021), arXiv:2108.00044 [nucl-ex] .
  * Adam et al. (2021) J. Adam et al. (STAR), Phys. Rev. Lett. 126, 162301 (2021), arXiv:2012.13601 [nucl-ex] .
  * Becattini et al. (2013) F. Becattini, V. Chandra, L. Del Zanna, and E. Grossi, Annals Phys. 338, 32 (2013), arXiv:1303.3431 [nucl-th] .
  * Fang et al. (2016) R.-h. Fang, L.-g. Pang, Q. Wang, and X.-n. Wang, Phys. Rev. C 94, 024904 (2016), arXiv:1604.04036 [nucl-th] .
  * Liu et al. (2020a) Y.-C. Liu, K. Mameda, and X.-G. Huang, Chin. Phys. C 44, 094101 (2020a), [Erratum: Chin.Phys.C 45, 089001 (2021)], arXiv:2002.03753 [hep-ph] .
  * Karpenko and Becattini (2017) I. Karpenko and F. Becattini, Eur. Phys. J. C 77, 213 (2017), arXiv:1610.04717 [nucl-th] .
  * Xie et al. (2016) Y. L. Xie, M. Bleicher, H. Stöcker, D. J. Wang, and L. P. Csernai, Phys. Rev. C 94, 054907 (2016), arXiv:1610.08678 [nucl-th] .
  * Xie et al. (2017) Y. Xie, D. Wang, and L. P. Csernai, Phys. Rev. C 95, 031901 (2017), arXiv:1703.03770 [nucl-th] .
  * Li et al. (2017) H. Li, L.-G. Pang, Q. Wang, and X.-L. Xia, Phys. Rev. C 96, 054908 (2017), arXiv:1704.01507 [nucl-th] .
  * Shi et al. (2019) S. Shi, K. Li, and J. Liao, Phys. Lett. B 788, 409 (2019), arXiv:1712.00878 [nucl-th] .
  * Xia et al. (2018) X.-L. Xia, H. Li, Z.-B. Tang, and Q. Wang, Phys. Rev. C 98, 024905 (2018), arXiv:1803.00867 [nucl-th] .
  * Wei et al. (2019) D.-X. Wei, W.-T. Deng, and X.-G. Huang, Phys. Rev. C 99, 014905 (2019), arXiv:1810.00151 [nucl-th] .
  * Vitiuk et al. (2020) O. Vitiuk, L. V. Bravina, and E. E. Zabrodin, Phys. Lett. B 803, 135298 (2020), arXiv:1910.06292 [hep-ph] .
  * Ivanov et al. (2019) Y. B. Ivanov, V. D. Toneev, and A. A. Soldatov, Phys. Rev. C 100, 014908 (2019), arXiv:1903.05455 [nucl-th] .
  * Fu et al. (2021a) B. Fu, K. Xu, X.-G. Huang, and H. Song, Phys. Rev. C 103, 024903 (2021a), arXiv:2011.03740 [nucl-th] .
  * Guo et al. (2021) Y. Guo, J. Liao, E. Wang, H. Xing, and H. Zhang, Phys. Rev. C 104, L041902 (2021), arXiv:2105.13481 [nucl-th] .
  * Li et al. (2021a) H. Li, X.-L. Xia, X.-G. Huang, and H. Z. Huang, (2021a), arXiv:2106.09443 [nucl-th] .
  * Deng et al. (2021) X.-G. Deng, X.-G. Huang, and Y.-G. Ma, (2021), arXiv:2109.09956 [nucl-th] .
  * Becattini and Karpenko (2018) F. Becattini and I. Karpenko, Phys. Rev. Lett. 120, 012302 (2018), arXiv:1707.07984 [nucl-th] .
  * Florkowski et al. (2019a) W. Florkowski, A. Kumar, R. Ryblewski, and A. Mazeliauskas, Phys. Rev. C 100, 054907 (2019a), arXiv:1904.00002 [nucl-th] .
  * Wu et al. (2019) H.-Z. Wu, L.-G. Pang, X.-G. Huang, and Q. Wang, Phys. Rev. Research. 1, 033058 (2019), arXiv:1906.09385 [nucl-th] .
  * Xie et al. (2020) Y. Xie, D. Wang, and L. P. Csernai, Eur. Phys. J. C 80, 39 (2020), arXiv:1907.00773 [hep-ph] .
  * Xia et al. (2019) X.-L. Xia, H. Li, X.-G. Huang, and H. Z. Huang, Phys. Rev. C 100, 014913 (2019), arXiv:1905.03120 [nucl-th] .
  * Becattini et al. (2019a) F. Becattini, G. Cao, and E. Speranza, Eur. Phys. J. C 79, 741 (2019a), arXiv:1905.03123 [nucl-th] .
  * Liu et al. (2020b) S. Y. F. Liu, Y. Sun, and C. M. Ko, Phys. Rev. Lett. 125, 062301 (2020b), arXiv:1910.06774 [nucl-th] .
  * Adam et al. (2019) J. Adam et al. (STAR), Phys. Rev. Lett. 123, 132301 (2019), arXiv:1905.11917 [nucl-ex] .
  * Acharya et al. (2021) S. Acharya et al. (ALICE), (2021), arXiv:2107.11183 [nucl-ex] .
  * Liu and Huang (2020) Y.-C. Liu and X.-G. Huang, Nucl. Sci. Tech. 31, 56 (2020), arXiv:2003.12482 [nucl-th] .
  * Becattini and Lisa (2020) F. Becattini and M. A. Lisa, Ann. Rev. Nucl. Part. Sci. 70, 395 (2020), arXiv:2003.03640 [nucl-ex] .
  * Gao et al. (2020) J.-H. Gao, G.-L. Ma, S. Pu, and Q. Wang, Nucl. Sci. Tech. 31, 90 (2020), arXiv:2005.10432 [hep-ph] .
  * (37) X.-G. Huang, J. Liao, Q. Wang, and X.-L. Xia, arXiv:2010.08937 [nucl-th] .
  * Gao and Liang (2019) J.-H. Gao and Z.-T. Liang, Phys. Rev. D 100, 056021 (2019), arXiv:1902.06510 [hep-ph] .
  * Weickgenannt et al. (2019) N. Weickgenannt, X.-L. Sheng, E. Speranza, Q. Wang, and D. H. Rischke, Phys. Rev. D 100, 056018 (2019), arXiv:1902.06513 [hep-ph] .
  * Hattori et al. (2019a) K. Hattori, Y. Hidaka, and D.-L. Yang, Phys. Rev. D 100, 096011 (2019a), arXiv:1903.01653 [hep-ph] .
  * Wang et al. (2019) Z. Wang, X. Guo, S. Shi, and P. Zhuang, Phys. Rev. D 100, 014015 (2019), arXiv:1903.03461 [hep-ph] .
  * Sheng et al. (2020) X.-L. Sheng, Q. Wang, and X.-G. Huang, Phys. Rev. D 102, 025019 (2020), arXiv:2005.00204 [hep-ph] .
  * Guo (2020) X. Guo, Chin. Phys. C 44, 104106 (2020), arXiv:2005.00228 [hep-ph] .
  * Weickgenannt et al. (2021a) N. Weickgenannt, E. Speranza, X.-l. Sheng, Q. Wang, and D. H. Rischke, Phys. Rev. Lett. 127, 052301 (2021a), arXiv:2005.01506 [hep-ph] .
  * Yang et al. (2020) D.-L. Yang, K. Hattori, and Y. Hidaka, JHEP 07, 070 (2020), arXiv:2002.02612 [hep-ph] .
  * Manuel and Torres-Rincon (2021) C. Manuel and J. M. Torres-Rincon, Phys. Rev. D 103, 096022 (2021), arXiv:2101.05832 [hep-ph] .
  * Sheng et al. (2021) X.-L. Sheng, N. Weickgenannt, E. Speranza, D. H. Rischke, and Q. Wang, Phys. Rev. D 104, 016029 (2021), arXiv:2103.10636 [nucl-th] .
  * Weickgenannt et al. (2021b) N. Weickgenannt, E. Speranza, X.-l. Sheng, Q. Wang, and D. H. Rischke, Phys. Rev. D 104, 016022 (2021b), arXiv:2103.04896 [nucl-th] .
  * Lin (2021) S. Lin, (2021), arXiv:2109.00184 [hep-ph] .
  * Florkowski et al. (2018) W. Florkowski, B. Friman, A. Jaiswal, and E. Speranza, Phys. Rev. C 97, 041901 (2018), arXiv:1705.00587 [nucl-th] .
  * Hattori et al. (2019b) K. Hattori, M. Hongo, X.-G. Huang, M. Matsuo, and H. Taya, Phys. Lett. B 795, 100 (2019b), arXiv:1901.06615 [hep-th] .
  * Montenegro and Torrieri (2019) D. Montenegro and G. Torrieri, Phys. Rev. D 100, 056011 (2019), arXiv:1807.02796 [hep-th] .
  * Florkowski et al. (2019b) W. Florkowski, A. Kumar, and R. Ryblewski, Prog. Part. Nucl. Phys. 108, 103709 (2019b), arXiv:1811.04409 [nucl-th] .
  * Bhadury et al. (2021a) S. Bhadury, W. Florkowski, A. Jaiswal, A. Kumar, and R. Ryblewski, Phys. Lett. B 814, 136096 (2021a), arXiv:2002.03937 [hep-ph] .
  * Bhadury et al. (2021b) S. Bhadury, W. Florkowski, A. Jaiswal, A. Kumar, and R. Ryblewski, Phys. Rev. D 103, 014030 (2021b), arXiv:2008.10976 [nucl-th] .
  * Garbiso and Kaminski (2020) M. Garbiso and M. Kaminski, JHEP 12, 112 (2020), arXiv:2007.04345 [hep-th] .
  * Gallegos and Gürsoy (2020) A. D. Gallegos and U. Gürsoy, JHEP 11, 151 (2020), arXiv:2004.05148 [hep-th] .
  * Shi et al. (2021) S. Shi, C. Gale, and S. Jeon, Phys. Rev. C 103, 044906 (2021), arXiv:2008.08618 [nucl-th] .
  * Fukushima and Pu (2021) K. Fukushima and S. Pu, Phys. Lett. B 817, 136346 (2021), arXiv:2010.01608 [hep-th] .
  * Li et al. (2021b) S. Li, M. A. Stephanov, and H.-U. Yee, Phys. Rev. Lett. 127, 082302 (2021b), arXiv:2011.12318 [hep-th] .
  * She et al. (2021) D. She, A. Huang, D. Hou, and J. Liao, (2021), arXiv:2105.04060 [nucl-th] .
  * Gallegos et al. (2021) A. D. Gallegos, U. Gürsoy, and A. Yarom, SciPost Phys. 11, 041 (2021), arXiv:2101.04759 [hep-th] .
  * Hongo et al. (2021) M. Hongo, X.-G. Huang, M. Kaminski, M. Stephanov, and H.-U. Yee, JHEP 11, 150 (2021), arXiv:2107.14231 [hep-th] .
  * Wang et al. (2021) D.-L. Wang, S. Fang, and S. Pu, (2021), arXiv:2107.11726 [nucl-th] .
  * Peng et al. (2021) H.-H. Peng, J.-J. Zhang, X.-L. Sheng, and Q. Wang, Chin. Phys. Lett. 38, 116701 (2021), arXiv:2107.00448 [hep-th] .
  * Zubarev et al. (1979) D. N. Zubarev, A. V. Prozorkevich, and S. A. Smolyanskii, Theo. Math. Phys. 40, 821 (1979).
  * van Weert (1982) C. G. van Weert, Ann. Phys. 140, 133 (1982).
  * Hongo (2017) M. Hongo, Annals Phys. 383, 1 (2017), arXiv:1611.07074 [hep-th] .
  * Becattini et al. (2019b) F. Becattini, M. Buzzegoli, and E. Grossi, Particles 2, 197 (2019b), arXiv:1902.01089 [cond-mat.stat-mech] .
  * Becattini et al. (2019c) F. Becattini, W. Florkowski, and E. Speranza, Phys. Lett. B 789, 419 (2019c), arXiv:1807.10994 [hep-th] .
  * Speranza and Weickgenannt (2021) E. Speranza and N. Weickgenannt, Eur. Phys. J. A 57, 155 (2021), arXiv:2007.00138 [nucl-th] .
  * Becattini et al. (2021a) F. Becattini, M. Buzzegoli, and A. Palermo, Phys. Lett. B 820, 136519 (2021a), arXiv:2103.10917 [nucl-th] .
  * Liu and Yin (2021) S. Y. F. Liu and Y. Yin, JHEP 07, 188 (2021), arXiv:2103.09200 [hep-ph] .
  * Becattini et al. (2021b) F. Becattini, M. Buzzegoli, G. Inghirami, I. Karpenko, and A. Palermo, (2021b), arXiv:2103.14621 [nucl-th] .
  * Fu et al. (2021b) B. Fu, S. Y. F. Liu, L. Pang, H. Song, and Y. Yin, Phys. Rev. Lett. 127, 142301 (2021b), arXiv:2103.10403 [hep-ph] .
  * Yi et al. (2021) C. Yi, S. Pu, and D.-L. Yang, Phys. Rev. C 104, 064901 (2021), arXiv:2106.00238 [hep-ph] .
  * Florkowski et al. (2021) W. Florkowski, A. Kumar, A. Mazeliauskas, and R. Ryblewski, (2021), arXiv:2112.02799 [hep-ph] .
  * Tinti and Florkowski (2020) L. Tinti and W. Florkowski, (2020), arXiv:2007.04029 [nucl-th] .
  * Hidaka et al. (2017) Y. Hidaka, S. Pu, and D.-L. Yang, Phys. Rev. D 95, 091901 (2017), arXiv:1612.04630 [hep-th] .
  * Huang et al. (2018) A. Huang, S. Shi, Y. Jiang, J. Liao, and P. Zhuang, Phys. Rev. D 98, 036010 (2018), arXiv:1801.03640 [hep-th] .
  * Liu et al. (2019) Y.-C. Liu, L.-L. Gao, K. Mameda, and X.-G. Huang, Phys. Rev. D 99, 085014 (2019), arXiv:1812.10127 [hep-th] .
  * Gao et al. (2018) J.-H. Gao, Z.-T. Liang, Q. Wang, and X.-N. Wang, Phys. Rev. D 98, 036019 (2018), arXiv:1802.06216 [hep-ph] .
  * (83) F. Becattini, arXiv:2004.04050 [hep-th] .
  * Becattini (2012) F. Becattini, Phys. Rev. Lett. 108, 244502 (2012), arXiv:1201.5278 [gr-qc] .
  * Becattini et al. (2021c) F. Becattini, M. Buzzegoli, A. Palermo, and G. Prokhorov, Phys. Lett. B 822, 136706 (2021c), arXiv:2009.13449 [hep-ph] .
  * Gao (2021) J.-H. Gao, Phys. Rev. D 104, 076016 (2021), arXiv:2105.08293 [hep-ph] .
  * Buzzegoli (2021) M. Buzzegoli, (2021), arXiv:2109.12084 [nucl-th] .


# Fake-Acknowledgment Attack on ACK-based Sensor Power Schedule for Remote
State Estimation

Yuzhe $\mathrm{Li}^{\star}$, Daniel E. $\mathrm{Quevedo}^{{\dagger}}$,
Subhrakanti $\mathrm{Dey}^{\diamond}$, and Ling $\mathrm{Shi}^{\star}$
$\star$: Electronic and Computer Engineering, Hong Kong University of Science
and Technology, Clear Water Bay, Kowloon, Hong Kong. Email: {yliah,
<EMAIL_ADDRESS>Department of Electrical Engineering (EIM-E), The
University of Paderborn, Germany, Email<EMAIL_ADDRESS>Department of Engineering Sciences, Uppsala University, Sweden, Email:
<EMAIL_ADDRESS>
###### Abstract

We consider a class of malicious attacks against remote state estimation. A
sensor with limited resources adopts an acknowledgement (ACK)-based online
power schedule to improve the remote state estimation performance. A malicious
attacker can modify the ACKs from the remote estimator and convey fake
information to the sensor. When the capability of the attacker is limited, we
propose an attack strategy for the attacker and analyze the corresponding
effect on the estimation performance. The possible responses of the sensor are
studied and a condition for the sensor to discard ACKs and switch from online
schedule to offline schedule is provided.

###### Index Terms:

Cyber-physical Systems, Security, Fake-ACK Attack, Remote State Estimation.

##  I Introduction

Cyber-physical systems (CPS) have been a hot topic among both academic and
industrial communities in the past few years. A wide application spectrum of
CPS can be found in areas such as smart grid, intelligent transportation and
environment monitoring [1] with the integration of sensing, control,
communication and computation. In most CPS infrastructures, wireless sensors
are key components with advantages such as low cost, easy installation, self-
power [2], when compared with traditional wired sensors. Therefore, wireless
sensors have been increasingly equipped in CPS to replace wired sensors. The
new issues due to their special characteristics have attracted much attention
in recent years.

The first issue is how to allocate the energy of the sensor efficiently. Since
most wireless sensors use on-board batteries which are difficult to replace,
and are typically expected to work for several years without battery
replacement, energy conservation is critical [3, 4, 5]. In [6], the authors
consider the sensor scheduling problem of whether it should send its data to a
remote estimator or not due to the limited available communication energy. The
optimal offline sensor schedule is derived, where the term “offline” means
that the sensor designs its strategy before the process starting, i.e.,
independent of the state of the process. In [7, 8, 9], the authors extend the
result in [6] and propose a so-called “online” sensor power schedules by
utilizing the real-time information of the process. One typical category of
online information is the ACK from the remote estimator, indicating whether
the data packet from the sensor arrives successfully or not. The ACK-based
online power schedule is proved to improve the estimation performance
significantly [8] compared to the offline schedule which only requires sending
a $1$-bit ACK packet.

ControlCommunicationComputationCyberPhysicalSystemInformation Figure 1:
Architecture of Cyber-Physical Systems.

However, as the communication is through a wireless channel, these sensors are
more vulnerable to cyber security threats. In some safety-critical
infrastructures, the wide use of CPS enlarges the risks and severities of such
attacks. For example, as the largest and most complex CPS in the future, any
severe attack on smart grid may have significant impacts on national economy
and security or even loss of human life [10]. Therefore, the security issue is
of fundamental importance to ensure the safe operation of CPS. There are two
possible types of attacks that are commonly investigated in the literature:
deception (integrity) attacks and denial-of-service (DoS) attacks [11],
corresponding to the two traditional security goals integrity and
availability. The integrity attacks on CPS focus on the integrity of data by
modifying the data packet [12, 13], while the DoS attacks block the exchange
of information including sensor measurements or control inputs between each
part of CPS [14, 15]. In practice, the integrity attacks studied in [12, 13,
16] often require comprehensive information about the system and modifications
of the data. Such information is not needed in DoS attacks. In this work, we
will consider both types of attacks.

Though the online power schedule in [8] improves the estimation performance,
the simple structure of the ACK packet makes it a more reachable (and likely)
alternative for the adversary, under both integrity attacks and DoS attacks.
Elucidating what is the possible attack patterns are and how they will affect
the estimation performance are important to help us improve the design of CPS.
This motivates us to investigate such potential class of malicious attacks
against ACK-based remote state estimation.

The remainder of the work is organized as follows. Section II presents the
system framework and states the main problem of interest. Section III provides
some preliminaries about the optimal offline schedule and an online schedule.
The analysis of our proposed fake-ACK attack pattern is given in Section IV.
Numerical examples and simulations are demonstrated in Section V. Section VI
provides some concluding remarks.

Notations: $\mathbb{Z}$ denotes the set of all integers and $\mathbb{N}$ the
positive integers. $\mathbb{R}$ is the set of real numbers. $\mathbb{R}^{n}$
is the $n$-dimensional Euclidean space. $\mathbb{S}_{+}^{n}$ (and
$\mathbb{S}_{++}^{n}$) is the set of $n$ by $n$ positive semi-definite
matrices (and positive definite matrices). When $X\in\mathbb{S}_{+}^{n}$ (and
$\mathbb{S}_{++}^{n}$), we write $X\geqslant 0$ (and $X>0$) and $X\geqslant Y$
if $X-Y\in\mathbb{S}_{+}^{n}$. The curled inequality symbols $\succeq$ and
$\preceq$ (and their strict forms $\succ$ and $\prec$) are used to denote
generalized componentwise inequalities between vectors: for vectors
$\mathbf{a}=[a_{1},a_{2},...,a_{n}]^{\prime},\mathbf{b}=[b_{1},b_{2},...,b_{n}]^{\prime}$,
we write $\mathbf{a}\succeq\mathbf{b}$ if $a_{i}\geqslant b_{i}$, for
$i=1,2,...,n$. $\mathbf{1}$ denotes vector with all entries one.
$\mathrm{Tr}(\cdot)$ is the trace of a matrix. The superscript ′ stands for
transposition. For functions $g,h$ with appropriate domains, $g\circ h(x)$
stands for the function composition $g\big{(}h(x)\big{)}$, and
$h^{n}(x)\triangleq h\big{(}h^{n-1}(x)\big{)}$, where $n\in\mathbb{N}$ and
with $h^{0}(x)\triangleq x$. $\delta_{ij}$ is the Dirac delta function, i.e.,
$\delta_{ij}$ equals to $1$ when $i=j$, and $0$ otherwise. The notation
$\mathbb{P}[\cdot]$ refers to probability and $\mathbb{E}[\cdot]$ to
expectation.

##  II Problem Setup

###  II-A System Model

Our interest lies in the security of a remote state estimation as depicted in
Fig. 2. Here we consider a general discrete-time linear time-invariant (LTI)
process of the form:

$\displaystyle x_{k+1}$ $\displaystyle=$ $\displaystyle Ax_{k}+w_{k},$ (1)
$\displaystyle y_{k}$ $\displaystyle=$ $\displaystyle Cx_{k}+v_{k},$ (2)

where $k\in\mathbb{N}$, $x_{k}\in\mathbb{R}^{n_{x}}$ is the process state
vector at time $k$, $y_{k}\in\mathbb{R}^{n_{y}}$ is the measurement taken by
the sensor, $w_{k}\in\mathbb{R}^{n_{x}}$ and $v_{k}\in\mathbb{R}^{n_{y}}$ are
zero-mean i.i.d. Gaussian noises with
$\mathbb{E}[w_{k}w_{j}^{\prime}]=\delta_{kj}Q$ ($Q\geqslant 0$),
$\mathbb{E}[v_{k}{v_{j}}^{\prime}]=\delta_{kj}R$ ($R>0$),
$\mathbb{E}[w_{k}{v_{j}}^{\prime}]=0\;\forall j,k\in\mathbb{N}$. The initial
state $x_{0}$ is a zero-mean Gaussian random vector uncorrelated with $w_{k}$
and $v_{k}$ with covariance $\Pi_{0}\geqslant 0$. The pair $(A,C)$ is assumed
to be observable and $(A,Q^{1/2})$ is controllable.

Figure 2: System Architecture.

We consider the so-called “smart sensor” as described in [17], which first
locally estimates the state $x_{k}$ based on all the measurements it has
collected up to time $k$ and then transmits its local estimate to the remote
estimator.

Denote $\hat{x}_{k}^{s}$ and $P_{k}^{s}$ as the sensor’s local minimum mean-
squared error (MMSE) estimate of the state $x_{k}$ and the corresponding error
covariance:

$\displaystyle\hat{x}_{k}^{s}$ $\displaystyle=$
$\displaystyle\mathbb{E}[x_{k}|y_{1},y_{2},...,y_{k}],$ (3)
$\displaystyle\hat{P}_{k}^{s}$ $\displaystyle=$
$\displaystyle\mathbb{E}[(x_{k}-\hat{x}_{k}^{s})(x_{k}-\hat{x}_{k}^{s})^{\prime}|y_{1},y_{2},...,y_{k}],$
(4)

which can be calculated by a standard Kalman filter.

As the estimation error covariance of the Kalman filter converges to a unique
value from any initial condition ([18]), without loss of generality, we assume
that the Kalman filter at the sensor side has entered the steady state and
simplify our subsequent discussion by setting:

$P_{k}^{s}=\overline{P},~{}k\geqslant 1,$ (5)

where $\overline{P}$ is the steady-state error covariance. For notational
ease, we define the Lyapunov and Riccati operators
$h,\tilde{g}:\mathbb{S}_{+}^{n}\rightarrow\mathbb{S}_{+}^{n}$ as:

$\displaystyle h(X)$ $\displaystyle\triangleq$ $\displaystyle AXA^{\prime}+Q,$
$\displaystyle\tilde{g}(X)$ $\displaystyle\triangleq$ $\displaystyle
X-XC^{\prime}[CXC^{\prime}+R]^{-1}CX.$

Then $\overline{P}$ is given by the unique positive semi-definite solution of
$\tilde{g}\circ h(X)=X$ (see [18]).

The error covariance $\overline{P}$ has the following property:

###### Lemma II.1

(Lemma 2.3 in [19]) For $0\leqslant t_{1}\leqslant t_{2}$, the following
inequality holds:

$h^{t_{1}}(\overline{P})\leqslant h^{t_{2}}(\overline{P}).$ (6)

In addition, if $t_{1}<t_{2}$, then

$\mathrm{Tr}\left(h^{t_{1}}(\overline{P})\right)<\mathrm{Tr}\left(h^{t_{2}}(\overline{P})\right).$
(7)

After obtaining $\hat{x}^{s}_{k}$, the sensor will transmit it as a data
packet to the remote estimator. Due to fading and interference, random data
drops will occur. As modelled in [6, 8], we assume that the sensor has two
transmission power levels: when using a higher energy $\Delta$, the packet
will always arrive to the remote estimator; when using a lower energy
$\delta$, the successful arrival rate is $\lambda\in(0,1)$. To simplify the
following discussion, we denote the sensor power schedule as:

$\theta=\\{\gamma_{1},\gamma_{2},...,\gamma_{k},...\\},$ (8)

where $\gamma_{k}=1$ or $0$ represents that the sensor chooses energy $\Delta$
or $\delta$, respectively, at time step $k$.

###  II-B Remote State Estimation

The transmission of $\hat{x}^{s}_{k}$ between the sensor and the remote
estimator can be characterized by a binary random process
$\\{\lambda_{k}\\},k\in\mathbb{N}$:

$\lambda_{k}=\begin{cases}1,&\text{if $\hat{x}^{s}_{k}$ arrives at time
$k$,}\\\ 0,&\text{otherwise (regarded as dropout).}\end{cases}$

Denote $\hat{x}_{k}$ and $P_{k}$ as the remote estimator’s own MMSE state
estimate and the corresponding error covariance based on all the sensor data
packets received up to time step $k$. From [20] and [21], the remote state
estimate $\hat{x}_{k}$ obeys the recursion:

$\hat{x}_{k}=\left\\{\begin{array}[]{ll}\hat{x}_{k}^{s},&\mathrm{if}~{}\lambda_{k}=1,\\\
A\hat{x}_{k-1},&\mathrm{if}~{}\lambda_{k}=0.\end{array}\right.$ (9)

The corresponding state estimation error covariance $P_{k}$ satisfies:

$P_{k}=\left\\{\begin{array}[]{ll}\overline{P},&\mathrm{if}~{}\lambda_{k}=1,\\\
h(P_{k-1}),&\mathrm{if}~{}\lambda_{k}=0.\end{array}\right.$ (10)

Intuitively, higher transmission power leads to better estimation performance.
However, in practice, the sensor has a limited energy budget, which motivates
us to consider the following optimization problem:

###### Problem II.2

$\displaystyle\min_{\theta}~{}~{}~{}J(\theta)\triangleq\limsup_{T\to\infty}\frac{1}{T}\sum_{k=1}^{T}\mathrm{Tr}\left(\mathbb{E}[P_{k}]\right),$
$\displaystyle\mathrm{s.t.}~{}~{}~{}\Phi(\theta)\triangleq\limsup_{T\to\infty}\frac{1}{T}\sum_{k=1}^{T}\mathbb{E}\left[\gamma_{k}\Delta+(1-\gamma_{k})\delta\right]\leqslant\Psi,$

where $\Psi$ is the average energy constraint for the sensor.

###### Remark II.3

We assume that $\Delta$, $\delta$, and $\Psi$ are all positive rational
numbers satisfying $\delta<\Psi<\Delta$. $\square$

##  III Preliminaries

In this section, we will revise some preliminary results about the optimal
offline schedule and an ACK-based online power schedule when there is no
attacker.

###  III-A Optimal Offline Schedule

When the sensor designs its power schedule without utilizing any online
information, the optimal solution is provided by the following theorem:

######  Theorem III.1 (Optimal Offline Schedule Proposed in [6])

Given the average energy constraint $\Psi$, suppose that $p$ and $q$ are two
co-prime integers satisfying $\frac{p}{q}=\frac{\Psi-\delta}{\Delta-\delta}$.
Then the optimal offline power schedule $\theta_{\text{off}}^{\star}$ to
Problem II.2 over a period $q$ can be constructed in the following form:

$\underbrace{(1\underbrace{0\cdots
0}_{s_{0}+1\text{~{}times}})\cdots(1\underbrace{0\cdots
0}_{s_{0}+1\text{times}})}_{n\text{~{}times}}\underbrace{(1\underbrace{0\cdots
0}_{s_{0}\text{~{}times}})\cdots(1\underbrace{0\cdots
0}_{s_{0}\text{times}})}_{m\text{~{}times}}$

where $1$ and $0$ denotes the designed value for $\gamma_{k}$ in (8),
$m=q-p(s_{0}+1)$, $n=p(s_{0}+2)-q$, and $s_{0}$ is the largest integer such
that $s_{0}\leqslant\frac{q}{p}-1$. $\square$

Under the optimal offline schedule $\theta_{\text{off}}^{\star}$, we have the
following closed-form results about the corresponding energy consumption and
estimation performance [6].

###### Proposition III.2

When $\theta_{\text{off}}^{\star}$ is used, the average energy cost is given
by

$\Phi(\theta_{\text{off}}^{\star})=\frac{(m+n)\Delta+(ms_{0}+ns_{0}+m)\delta}{ms_{0}+ns_{0}+2m+n}=\Psi,$

and the trace of the expected average state estimation error covariance is
given by:

$\displaystyle J(\theta_{\text{off}}^{\star})=$
$\displaystyle\frac{1}{m(s_{0}+2)+n(s_{0}+1)}\cdot$
$\displaystyle\mathrm{Tr}\Big{\\{}m\big{[}1+\lambda(s_{0}+1)]\overline{P}+n\big{[}1+\lambda
s_{0})]\overline{P}$
$\displaystyle+m\sum_{i=1}^{s_{0}+1}\big{[}1+\lambda(s_{0}+1-i)\big{]}(1-\lambda)^{i}h^{i}(\overline{P})$
$\displaystyle+n\sum_{i=1}^{s_{0}}\big{[}1+\lambda(s_{0}-i)\big{]}(1-\lambda)^{i}h^{i}(\overline{P})\Big{\\}}.$

$\square$

###  III-B Online Power Schedule

The offline power schedules are designed before the system is running without
using the realtime measurement or state information. In [9], the author showed
that the utilization of the online information can improve the system
performance while require less energy consumption. To further improve the
estimation performance, an online power schedule based on the ACKs from the
remote estimator was proposed in [8] as follows: the remote estimator
generates $1$-bit ACKs ($=\lambda_{k}$) to indicate whether the data packet
arrives successfully or not and an event-detector at the remote estimator’s
side collects and stores the ACKs in an $L$-bit memory. Without loss of
generality, assume the sensor uses $\Delta$ at the first time step. The memory
is set to $11,...,11$ initially and the detector randomly chooses to activate
$z_{0}$-bit memory ($z_{0}<L$) with probability $\mu$, or $(z_{0}+1)$-bit
memory otherwise. At every time step, the memory shifts all bits one bit to
the MSB direction with the existing MSB being deleted and the incoming ACK
stored in the LSB. When the memory becomes $00,...,00$, the detector sends a
flag-ACK to inform the sensor to use high power in the following time step. In
the meanwhile, the memory is set to $11,...,11$ again (see Fig. 3).

Figure 3: System Architecture under Online Sensor Schedule.

By properly designing the parameter $z_{0}$ and $\mu$, the average energy
consumption can meet the constraint in Problem II.2 and the corresponding
estimation performance is given by the following theorem:

######  Theorem III.3 (Theorem 4.5 in [8])

Under the proposed online power schedule $\theta_{\text{on}}$, the trace of
the average expected estimation error covariance is given by:

$\displaystyle J(\theta_{\text{on}})=$
$\displaystyle\frac{\lambda}{1-\mu(1-\lambda)^{z_{0}+1}-\mu(1-\lambda)^{z_{0}+2}}\cdot$
$\displaystyle\mathrm{Tr}\Big{[}\sum_{i=0}^{z_{0}}(1-\lambda)^{i}h^{i}(\overline{P})+\mu(1-\lambda)^{z_{0}+1}h^{z_{0}+1}(\overline{P})\Big{]}.$

$\square$

Compared to the optimal offline schedule, the proposed online schedule
$\theta_{\text{on}}$ has a better estimation performance under the same energy
constraint, which is summarized in the following theorem [8]:

######  Theorem III.4 (Theorem 5.3 in [8])

Under the same energy constraint as in Problem II.2, the proposed online
schedule $\theta_{\text{on}}$ has a better estimation performance than the
optimal offline schedule $\theta_{\text{off}}^{\star}$ provided in Theorem
III.1, i.e.,

$J(\theta_{\text{on}})<J(\theta_{\text{off}}^{\star}).$

##  IV Fake-acknowledgement Attack Against Remote Estimation

As stated before, the ACK-based online power schedule can improve the
estimation performance significantly compared to the offline schedule.
However, the simple structure of the $1$-bit flag-ACK packet makes it a more
reachable (and likely) alternative for the adversary, under both integrity
attacks and DoS attacks. In this section, we propose a possible attack pattern
for the attacker and investigate the corresponding effect on the estimation
performance.

###  IV-A Proposed Fake-ACK Attack Pattern

Based on the communication protocol used by the event-detector, the attacker
may adopt different types of attacks: when the flag-ACK is sent at every time
step with value $1$ or $0$ (Type A), the attacker modifies the bit contained
in the packet to launch the attack, i.e., the integrity attack; when the flag-
ACK is sent only when the event happens (Type B), the attacker can simply
adopt DoS attack to prevent the transmission of the packet. No matter what
types of the communication protocols and attacks are used, the consequent
results are the same: the time steps which would use energy $\Delta$ now use
$\delta$ instead under the attack. Therefore, without loss of generality, we
will assume that the event-detector adopts protocol Type B and the attacker
launches DoS attack in the following discussion (the analysis of the other
case is almost the same).

Assume that the attacker also has an energy constraint and cannot launch
attack at every time step. Otherwise, it is straightforward to see that the
attacker will attack all the flag-ACKs with value $1$ to prevent the
transmission. Let the proportion of the flag-ACK packets that can be prevented
by the attacker among all the flag-ACK packets be constrained to be
$\beta\in(0,1)$ and $t,r(t>r)$ are two co-prime integers satisfying
$\frac{r}{t}=\beta$. In [15], the optimal DoS attack pattern is proved to be
launching the attack consecutively so that the effect on the estimation
performance is maximized. This motivates us to consider the specific attack
pattern of interest described as follows:

  1. 1.

the attacker follows a periodic pattern;

  2. 2.

the attacker maintains a counter and launches DoS attacks to prevent the
transmission of the first $r$ coming flag-ACKs, and leave the rest $t-r$
coming flag-ACK with value $1$ unchanged;

  3. 3.

after every $t$ coming flag-ACKs, the counter is set to $0$ again.

This is summarized in Algorithm 1.

Algorithm 1 Attacking Pattern for The Attacker

1:Process begins;

2:$counter=0$;

3:while 1 do

4: if a flag-ACK packet is sent then

5: if $counter\leqslant r$ then

6: block the transmission of the flag-ACK packet;

7: else

8: do not block;

9: end if

10: $counter=counter+1$;

11: end if

12: if $counter==t$ then

13: $counter=0$;

14: end if

15:end while

###### Example IV.1

We use a simple example to illustrate how this attack pattern works. Assume
that $z_{0}=2$, $\beta=\frac{2}{3}$ (thus $r=2$ and $t=3$), one possible
realization and comparison between the online schedule $\theta_{\text{on}}$
and the one under attack $\tilde{\theta}_{\text{on}}$ are shown in Fig. 4 (the
long arrows represent the arrival packets while the short one is for the
dropped packets; the number between arrows is the value of the counter). In
the realization of $\theta_{\text{on}}$ without attack, the second and the
third flag-ACK are sent due to two ($z_{0}=2$) consecutive packet losses. When
there is an attacker following the pattern in Algorithm 1, it will attack the
first two ($r=2$) flag-ACKs, let the third one transmitted, and reset the
counter for a new period.

Figure 4: Comparison between $\theta_{\text{on}}$ and
$\tilde{\theta}_{\text{on}}$

###  IV-B Performance Analysis

Denote the estimation performance of $\theta_{\text{on}}$ under the proposed
fake-ACK attack as $J(\tilde{\theta}_{\text{on}})$. We will investigate the
analytical form of $J(\tilde{\theta}_{\text{on}})$ in the following part.

Note that due to the recursion of the dynamics in (10), the covariance $P_{k}$
can only take value in the infinitely countable set
$\\{\overline{P},h(\overline{P}),h^{2}(\overline{P}),...\\}$. Denote
$\tau_{k}\in\mathbb{Z}$ as the holding time since the most recent time when
the remote estimator successfully received the data from sensor:

$\tau_{k}\triangleq k-\max_{1\leqslant t\leqslant k}\\{t:\lambda_{t}=1\\}.$
(11)

Then we have

$P_{k}=\left\\{\begin{array}[]{ll}\overline{P},&\mathrm{if}~{}\lambda_{k}=1,\\\
h^{\tau_{k}}(\overline{P}),&\mathrm{if}~{}\lambda_{k}=0.\end{array}\right.$

Define the value of the attacker’s counter at time step $k$ as
$\sigma_{k}\in\\{0,1,2,...,t-1\\}$ and the state of the process at time step
$k$ as:

$S_{k}={(\tau_{k},\sigma_{k})}.$

As the strategy of the attacker guarantees that the flag-ACK will be sent to
the sensor when the value of the counter is larger than $r$, the covariance
$P_{k}$ will be bounded and $\tau_{k}$ only takes value in a finite set.
Therefore, based on the recursion in (10), it is easy to verify that the
process $\\{S_{k}\\}$ constitutes a finite state stationary Markov chain.

###### Remark IV.2

Note that the design of the parameter $\mu$ in the online sensor schedule
$\theta_{\text{on}}$ is to tune the activated memory length between $z_{0}$
and $z_{0}+1$ in the event-detector to satisfy the energy constraint. When the
energy constraint is in some simple form, e.g., $\frac{1}{2}$ or
$\frac{1}{3}$, a fixed length activated memory with $z_{0}$ is sufficient.
Without loss of generality, we assume that the event-detector adopts a fixed
activated memory with length $z_{0}$ to facilitate the discussion and simplify
the notation below. We will later show that our results also apply to the
general case with variant length. $\square$

According to the mechanism of the online sensor schedule $\theta_{\text{on}}$
and the attack strategy, $\tau_{k}$ only takes values in the set
$\\{0,1,2,...,z_{0}(\sigma_{k}+1)\\}$ when $\sigma_{k}\leqslant r$ and in the
set $\\{0,1,2,...,z_{0}\\}$ when $r<\sigma_{k}<t$. Then the state set can be
constructed as:

$\displaystyle\mathcal{S}\triangleq\Big{\\{}$
$\displaystyle(0,0),(1,0),..,(z_{0}-1,0),(z_{0},0),$
$\displaystyle(0,1),(1,1),...,(2z_{0}-1,1),(2z_{0},1),$ $\displaystyle...$
$\displaystyle(0,r),(1,r),...,\big{(}z_{0}(r+1)-1,r\big{)},\big{(}z_{0}(r+1),r\big{)},$
$\displaystyle(0,r+1),(1,r+1),..,(z_{0}-1,r+1),(z_{0},r+1),$
$\displaystyle(0,r+2),(1,r+2),..,(z_{0}-1,r+2),(z_{0},r+2),$
$\displaystyle...$
$\displaystyle(0,t-1),(1,t-1),..,(z_{0}-1,t-1),(z_{0},t-1)\Big{\\}}.$

The total number of states can be calculated as:

$\displaystyle|\mathcal{S}|$
$\displaystyle=\sum_{i=0}^{r}\big{[}z_{0}(i+1)+1\big{]}+\big{[}(t-1)-(r+1)+1\big{]}(z_{0}+1)$
$\displaystyle=\frac{1}{2}r(r+1)z_{0}+t(z_{0}+1).$

It is convenient to label each state $\epsilon_{i}$ in $\mathcal{S}$ in a
sequential order from $i=1$ to $|\mathcal{S}|$. Then we can define the state
transition probability matrix as:

$\mathcal{T}_{s}=\\{T_{ij}\\}_{|\mathcal{S}|\times|\mathcal{S}|},$

where

$T_{ij}\triangleq\mathbb{P}[S_{k+1}=\epsilon_{j}|S_{k}=\epsilon_{i}].$

Based on the mechanism of the online sensor schedule $\theta_{\text{on}}$ and
the attack strategy, we can derive the transition probability from
$S_{k+1}=(\tau_{k+1},\sigma_{k})$ to $S_{k}=(\tau_{k},\sigma_{k})$ as follows:

$1)$ when $\tau_{k}\mod z_{0}\neq 0$ or $\tau_{k}=0$, the event detector will
not send the flag-ACK (and no attack as a consequence). Thus:

$\displaystyle\mathbb{P}\big{[}S_{k+1}|S_{k}\big{]}=\begin{cases}\lambda,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(0,\sigma_{k}),\\\
1-\lambda,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(\tau_{k}+1,\sigma_{k}),\\\
0,&\text{otherwise};\end{cases}$

$2)$ when $\tau_{k}\mod z_{0}=0$, $\tau_{k}\neq 0$ and $\sigma_{k}<r$, the
memory satisfies the triggering condition, the flag-ACK will be sent and the
attack will be launched with the counter increasing by $1$. Thus:

$\displaystyle\mathbb{P}\big{[}S_{k+1}|S_{k}\big{]}=\begin{cases}\lambda,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(0,\sigma_{k}+1),\\\
1-\lambda,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(\tau_{k}+1,\sigma_{k}+1),\\\
0,&\text{otherwise};\end{cases}$

$3)$ when $\tau_{k}\mod z_{0}=0$, $\tau_{k}\neq 0$ and
$r\leqslant\sigma_{k}<t-1$, the memory satisfies the triggering condition, the
flag-ACK will be sent without being attacked and the counter increasing by
$1$. Thus:

$\displaystyle\mathbb{P}\big{[}S_{k+1}|S_{k}\big{]}=\begin{cases}1,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(0,\sigma_{k}+1),\\\
0,&\text{otherwise};\end{cases}$

$4)$ when $\tau_{k}\mod z_{0}=0$, $\tau_{k}\neq 0$ and $\sigma_{k}=t-1$, the
memory satisfies the triggering condition, the flag-ACK will be sent without
being attacked and the counter is reset to $0$. Thus:

$\displaystyle\mathbb{P}\big{[}S_{k+1}|S_{k}\big{]}=\begin{cases}1,&\text{if~{}}(\tau_{k+1},\sigma_{k+1})=(0,0),\\\
0,&\text{otherwise}.\end{cases}$

Based on the above discussion, we can obtain the expression of the transition
matrix $\mathcal{T}_{s}$. Due to space limitation, we provide a simple example
of $\mathcal{T}_{s}$ when $z_{0}=2$, $r=1$ and $t=3$ (therefore
$|\mathcal{S}|=11$) in (12).

$\mathcal{T}_{s}=\left[\begin{array}[]{ccccccccccc}\bm{\lambda}&\bm{\lambda}&\bm{\lambda}&0&0&0&0&0&0&0&\bm{1}\\\
\bm{1-\lambda}&0&0&0&0&0&0&0&0&0&0\\\ 0&\bm{1-\lambda}&0&0&0&0&0&0&0&0&0\\\
0&0&0&\bm{\lambda}&\bm{\lambda}&0&\bm{\lambda}&0&0&0&0\\\
0&0&0&\bm{1-\lambda}&0&0&0&0&0&0&0\\\ 0&0&0&0&\bm{1-\lambda}&0&0&0&0&0&0\\\
0&0&\bm{1-\lambda}&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&\bm{1-\lambda}&0&0&0&0\\\
0&0&0&0&0&\bm{1}&0&\bm{1}&\bm{\lambda}&\bm{\lambda}&0\\\
0&0&0&0&0&0&0&0&\bm{1-\lambda}&0&0\\\
0&0&0&0&0&0&0&0&0&\bm{1-\lambda}&0\end{array}\right]$ (12)

Based on the well-established ergodic Markov chain theory [22], $\\{S_{k}\\}$
has a stationary probability distribution of states, which is denoted as:

$\Pi=\Big{[}\mathbb{P}[\epsilon_{1}],\mathbb{P}[\epsilon_{2}],...,\mathbb{P}[\epsilon_{|\mathcal{S}|}]\Big{]}^{\prime},$

When $\\{S_{k}\\}$ enters the steady state, we have

$\begin{cases}&\mathcal{T}_{s}\Pi=\Pi,\\\
&\sum_{i=1}^{|\mathcal{S}|}\Pi(i)=1,\end{cases}$ (13)

where $\Pi(i)$ is the $i$-th element of $\Pi$.

By solving (13), we can obtain the stationary distribution of state, denoted
as $\Pi^{\star}$. Then we can calculate the estimation performance of
$\theta_{\text{on}}$ under the proposed fake-ACK attack as follows:

$\displaystyle J(\tilde{\theta}_{\text{on}})=$
$\displaystyle\limsup_{T\to\infty}\frac{1}{T}\sum_{k=1}^{T}\mathrm{Tr}\left(\mathbb{E}[P_{k}]\right)$
$\displaystyle=$
$\displaystyle\limsup_{k\to\infty}\mathrm{Tr}\left(\mathbb{E}[P_{k}]\right)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{|\mathcal{S}|}\Pi^{\star}(i)\mathrm{Tr}\Big{(}\mathbb{E}\big{[}P_{k}(\epsilon_{i})\big{]}\Big{)}.$

###  IV-C Discussion

To analyze the effect of the proposed attacking pattern on
$\theta_{\text{on}}$, we first consider two extreme cases.

When the energy constraint for the attacker is $\beta=1$, the attacker can
block the transmissions of all the flag-ACKs, resulting in the sensor using
the low power $\delta$ within the entire time-horizon. The corresponding
estimation performance can be calculated as:

$\displaystyle J(\tilde{\theta}_{\text{on}})_{\text{max}}=$
$\displaystyle\limsup_{T\to\infty}\frac{1}{T}\sum_{k=1}^{T}\mathrm{Tr}\left(\mathbb{E}[P_{k}]\right)$
$\displaystyle=$
$\displaystyle\limsup_{k\to\infty}\mathrm{Tr}\left(\mathbb{E}[P_{k}]\right)$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{+\infty}\lambda(1-\lambda)^{i}\mathrm{Tr}\big{(}h^{i}(\overline{P})\big{)}.$

Note that the necessary and sufficient condition for the stability of the
above case is simply given by $\rho(A)(1-\lambda)<1$, where $\rho(A)$ is the
maximum eigenvalue of $A$. Since, as proved in [20] and [6], more energy
budget is always beneficial, we have:

$J(\tilde{\theta}_{\text{on}})_{\text{max}}=J(\theta_{\text{off}}^{\star})_{(n=0)}\geqslant
J(\theta_{\text{off}}^{\star}),$ (14)

where $n$ is the number of high energy $\delta$ that can be used within each
period of $J(\theta_{\text{off}}^{\star})$ defined in Theorem III.1.

One the other hand, when $\beta=0$, the attacker cannot block any flag-ACK and
the online sensor schedule will operate normally. Therefore, we have:

$J(\tilde{\theta}_{\text{on}})_{\text{min}}=J(\theta_{\text{on}}).$ (15)

Combining (14) and (15), we have:

$J(\theta_{\text{on}})=J(\tilde{\theta}_{\text{on}})_{\text{min}}\leqslant
J(\theta_{\text{off}}^{\star})\leqslant
J(\tilde{\theta}_{\text{on}})_{\text{max}}.$

As a consequence, when the capability of the attacker is limited, keeping the
online schedule will still guarantee a better performance than the offline
schedule $\theta_{\text{off}}^{\star}$. Otherwise, when the attacker has
sufficient energy budget, the estimation performance under the online sensor
schedule $\theta_{\text{on}}$, $J(\tilde{\theta}_{\text{on}})$, will be even
worse than the offline schedule $\theta_{\text{off}}^{\star}$. In such case,
not using ACK is preferable to using fake ACK and the reasonable choice for
the sensor is to adopt the offline schedule rather than keeping the online
schedule. Furthermore, as the performance metric
$J(\tilde{\theta}_{\text{on}})$ is growing with the energy budget $\beta$,
there exists a threshold value $\overline{\beta}$ such that (as later shown in
Fig. 6):

$J(\tilde{\theta}_{\text{on}})_{\beta=\overline{\beta}}=J(\theta_{\text{off}}^{\star}).$

Hence, in practice, based on the information about the attacker’s capability,
the sensor can choose different types of schedules to obtain a better
performance, i.e.,

$\tilde{\theta}_{\text{on}}=\begin{cases}\theta_{\text{on}},&\text{when~{}}\beta<\overline{\beta},\\\
\theta_{\text{off}}^{\star},&\text{when~{}}\beta\geqslant\overline{\beta},\end{cases}$

where $\overline{\beta}$ can be calculated numerically. How to obtain the
value of $\overline{\beta}$ analytically is an interesting problem which will
be considered in our future work.

###### Remark IV.3

Note that the sensor can detect the existence of the attacker by checking the
arrival rate of flag-ACKs based on the theoretical arrival rate provided in
[8]. The observation will enable the sensor to estimate the capability of the
attacker. On the other hand, when the sensor is aware of the existence of the
attacker, the attacker will also need to re-design his strategy to avoid being
detected, which will result in a dynamic decision-making process between both
sides. We will investigate this problem under a game-theoretic framework in
the future. $\square$

##  V Simulation

In this section, we provide numerical examples to illustrate our results in
different situations.

Define

$J_{k}(\theta)=\frac{1}{k}\sum_{i=1}^{k}\mathrm{Tr}\left(\mathbb{E}[P_{i}]\right),$

as the empirical approximation (via 100000 Monte Carlo simulations) of
$J(\theta)$ at every time instant $k$.

Consider a scalar system with parameters $A=1.2$, $C=0.7$, $R=Q=0.8$,
$\lambda=0.5$. Suppose that the energy constraint for the sensor is given by
$\Psi=\frac{1}{7}\Delta+\frac{6}{7}\delta$. It is easy to verify that the
optimal offline schedule is given in a periodic form of
$\\{1000000,1000000,...\\}$ with performance
$J(\theta^{\star}_{\text{off}})=2.0953$. The parameter for the online schedule
is $z_{0}=2$ with performance $J(\theta_{\text{on}})=1.6399$. As shown in Fig.
5, when the energy constraint for the attacker is $\beta=\frac{1}{5}$, the
estimation performance $J(\tilde{\theta}_{\text{on}})$ is better than the
performance using the offline schedule $\theta_{\text{off}}^{\star}$ (thus the
sensor may still choose $\theta_{\text{on}}$). However, when the energy budget
for the attacker is further increased, e.g., $\beta=\frac{2}{3}$,
$J(\tilde{\theta}_{\text{on}})$ is worse even than the offline
$\theta_{\text{off}}^{\star}$ (thus the sensor may switch from the online
schedule to the offline one to obtain better estimation performance). In Fig.
6, the comparison of $J(\tilde{\theta}_{\text{on}})$ as a function of $\beta$
is provided. As we can see, the approximate value of $\overline{\beta}$ is
about $0.28$.

Figure 5: Comparison between $J(\theta^{\star}_{\text{off}})$,
$J(\theta_{\text{on}})$ and $J(\tilde{\theta}_{\text{on}})$. Figure 6:
Comparison of $J(\tilde{\theta}_{\text{on}})$ under different $\beta$.

##  VI Conclusion

In this work, a potential class of malicious attacks against remote state
estimation was studied. We proposed a strategy for the attacker which can
modify the flag-ACKs from the remote estimator and convey fake information to
the sensor. The corresponding effect on the estimation performance is analyzed
based on well-established Markov chain theory. The condition for the sensor to
use an online schedule instead of an offline schedule is provided. Simulations
were provided to illustrate our results. Future works include obtaining the
threshold value of $\overline{\beta}$ analytically, considering other types of
fake-ACK attacks and extending the problem into a game-theoretic framework
where the attacker need to design their pattern without being detected.

## References

  * [1] J. P. Hespanha, P. Naghshtabrizi, and Y. Xu, “A survey of recent results in networked control systems,” _Proceedings of the IEEE_ , vol. 95, no. 1, pp. 138–162, 2007.
  * [2] V. C. Gungor and G. P. Hancke, “Industrial wireless sensor networks: Challenges, design principles, and technical approaches,” _IEEE Transactions on Industrial Electronics_ , vol. 56, no. 10, pp. 4258–4265, 2009\.
  * [3] D. E. Quevedo, A. Ahlén, A. S. Leong, and S. Dey, “On Kalman filtering over fading wireless channels with controlled transmission powers,” _Automatica_ , vol. 48, no. 7, pp. 1306–1316, 2012.
  * [4] O. C. Imer and T. Ba¸sar, “Optimal estimation with limited measurements,” in _Proceedings of the IEEE Conference on Decision and Control and the European Control Conference_ , 2005, pp. 1029–1034.
  * [5] A. Leong and S. Dey, “Power allocation for error covariance minimization in Kalman filtering over packet dropping links,” in _IEEE 51st Annual Conference on Decision and Control (CDC)_ , Dec 2012, pp. 3335–3340.
  * [6] L. Shi, P. Cheng, and J. Chen, “Sensor data scheduling for optimal state estimation with communication energy constraint,” _Automatica_ , vol. 47, no. 8, pp. 1693–1698, 2011.
  * [7] Y. Li, D. E. Quevedo, V. Lau, and L. Shi, “Online sensor transmission power schedule for remote state estimation,” in _IEEE 52nd Aunnal Conference on Decision and Control (CDC), Florence, Italy_ , 2013.
  * [8] D. Han, P. Cheng, J. Chen, and L. Shi, “An online sensor power schedule for remote state estimation with communication energy constraint,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 7, pp. 1942–1947, 2014.
  * [9] J. Wu, K. Johansson, and L. Shi, “An improved hybrid sensor schedule for remote state estimation under limited communication resources,” in _Proceedings of the 51st IEEE Conference on Decision and Control_ , 2012, pp. 3305–3310.
  * [10] Y. Mo, T. Kim, K. Brancik, D. Dickinson, H. Lee, A. Perrig, and B. Sinopoli, “Cyber-physical security of a smart grid infrastructure,” _Proceedings of the IEEE_ , vol. 100, no. 1, pp. 195–209, 2012.
  * [11] A. Cardenas, S. Amin, and S. Sastry, “Secure control: Towards survivable cyber-physical systems,” in _28th International Conference on Distributed Computing Systems Workshops_ , 2008, pp. 495–500.
  * [12] Y. Liu, P. Ning, and M. Reiter, “False data injection attacks against state estimation in electric power grids,” _ACM Transactions on Information and System Security (TISSEC)_ , vol. 14, no. 1, p. 13, 2011.
  * [13] Y. Mo and B. Sinopoli, “Integrity attacks on cyber-physical systems,” in _Proceedings of the 1st international conference on High Confidence Networked Systems_. ACM, 2012, pp. 47–54.
  * [14] Y. Li, L. Shi, P. Cheng, J. Chen, and D. E. Quevedo, “Jamming attack on cyber-physical systems: A game-theoretic approach,” in _IEEE International Conference on CYBER Technology in Automatation, Control, and Intelligent Systems, Nanjing, China_ , 2013.
  * [15] H. Zhang, P. Cheng, L. Shi, and J. Chen, “Optimal dos attack policy against remote state estimation,” in _IEEE 52nd Annual Conference on Decision and Control (CDC)_. IEEE, 2013, pp. 5444–5449.
  * [16] C.-Z. Bai and V. Gupta, “On Kalman filtering in the presence of a compromised sensor: Fundamental performance bounds,” in _American Control Conference_. IEEE, 2014, pp. 3029–3034.
  * [17] P. Hovareshti, V. Gupta, and J. S. Baras, “Sensor scheduling using smart sensors,” in _Proceedings of the 46th IEEE Conference on Decision and Control_ , 2007, pp. 494–499.
  * [18] B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M. I. Jordan, and S. S. Sastry, “Kalman filtering with intermittent observations,” _IEEE Transactions on Automatic Control_ , vol. 49, no. 9, pp. 1453–1464, 2004.
  * [19] L. Shi, K. H. Johansson, and L. Qiu, “Time and event-based sensor scheduling for networks with limited communication resources,” in _World Congress of the International Federation of Automatic Control (IFAC)_ , vol. 18, no. 1, 2011, pp. 13 263–13 268.
  * [20] Y. Li, D. E. Quevedo, V. Lau, and L. Shi, “Optimal periodic transmission power schedules for remote estimation of ARMA processes,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 24, pp. 6164–6174, 2013.
  * [21] L. Shi, M. Epstein, and R. M. Murray, “Kalman filtering over a packet-dropping network: A probabilistic perspective,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 3, pp. 594–604, 2010.
  * [22] J. R. Norris, _Markov chains_. Cambridge University Press, 1998, no. 2.