©2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works. DOI:
10.1109/CSR54599.2022.9850338

# A Bayesian Model Combination based approach to Active Malware Analysis
††thanks: This research is supported by the European Union’s Horizon 2020
Research and Innovation program CONCORDIA under Grant Agreement No. 830927

1st Abhilash Hota  Computer Science
Jacobs University
Bremen, Germany
<EMAIL_ADDRESS>2nd Jürgen Schönwälder  Computer Science
Jacobs University
Bremen, Germany
<EMAIL_ADDRESS>
###### Abstract

Active Malware Analysis involves modeling malware behavior by executing
actions to trigger responses and explore multiple execution paths. One of the
aims is making the action selection more efficient. This paper treats Active
Malware Analysis as a Bayes-Active Markov Decision Process and uses a Bayesian
Model Combination approach to train an analyzer agent. We show an improvement
in performance against other Bayesian and stochastic approaches to Active
Malware Analysis.

###### Index Terms:

Bayesian Model Combination, Bayes-active Markov decision process, active
malware analysis

##  I Introduction

Mobile devices have become increasingly ubiquitous and used in various
critical functions such as banking. The Android OS has a huge share of the
smartphone market [1] and as such, has become a preferred target for malware
authors in recent years. The number of new malwares has been steadily
increasing, with the AV-Test institute registering an average of around 250000
new samples for analysis every month last year [2]. Users of Android OS
devices have access to multiple marketplaces, aside from the official Google
Play store, from where they can download and install applications. Users can
also download and install applications not currently hosted on the
marketplaces. The process of malware analysis generally tends to involve
largely manual work by analysts. However, analyzing the large number of new
applications being released on various Android application marketplaces,
forces malware analysts to rely on various techniques to automate the malware
analysis process.

Malware analysts try to model a malware’s behavior on execution and use this
information to develop signatures that can then be used for malware detection
on a network or device. Approaches to malware analysis can be broadly
classified as static or dynamic. In static analysis [3] an application is
examined without executing the actual instructions. This approach is however
vulnerable to various evasion techniques employed by malware authors like
packing or code obfuscation[4]. Dynamic analysis, on the other hand, involves
executing an application in a sandbox environment and recording the execution
traces. Traditionally, dynamic analysis has been passive in that an
application is simply executed and its behavior observed. A lot of application
behaviors however, especially with mobile applications, are often only
triggered by user interaction.

Active Malware Analysis (AMA) involves modeling a malware sample’s behavior by
executing actions like button clicks, sending messages and microphone
activation in order to trigger responses and explore multiple execution paths.
A major focus in this approach is how to efficiently trigger different
behaviors from an application that might require user interaction in order to
explore the execution paths. Initial efforts in AMA have been based either on
past recorded user activity [5], or a pseudo-random selection of possible
triggering actions [6, 7]. More recent work has focused on employing machine
learning techniques to make smarter choices about which actions to use to
trigger these responses.

Our work models AMA as a Bayes-Adaptive Markov Decision Process (BAMDP) [8]
and applies Bayesian Model Combination (BMC) and $\epsilon$-BMC to the
problem. The intent here is to explore the use of model-free reinforcement
learning using the epsilon-greedy exploration policy to resolve the
exploration-exploitation trade-off in AMA. We evaluate the approach by
training the model on a dataset of 15000 Android malware samples and provide a
comparison against BAMA[9] and Bayes-adaptive Monte Carlo Planning as
approaches to AMA. Comparison is also provided against a stochastic game based
approach using Monte Carlo Tree Search models that have been trained on the
same dataset.

The rest of the paper is structured as follows: Section II provides some
background information and discusses related work. Section III the model
designs implemented in this paper. Section IV explains the experimental setup
and the process followed by the agent for malware analysis. Section V provides
an evaluation of the results and we conclude the paper in Section VI.

##  II Background and Related Work

In this section we provide the relevant background information on Markov
decision processes, Bayesian model combination and AMA.

###  II-A Markov Decision Process

The Markov Decision Process (MDP) is an approach to sequential decision-making
in Markovian dynamical systems [10]. It tracks a system state that can change
according to performed actions and affects the outcome of the system.

An MDP M is a tuple $(S,A,P,P_{0},q)$ where

  * •

S is the set of states,

  * •

A is the set of actions,

  * •

$P(\cdot|s,a)\in P(S)$ is the probability distribution over next states, when
action $a$ is taken in state $s$,

  * •

$P_{0}\in P(S)$ is the probability distribution according to which the initial
state is selected, and

  * •

$R(s,a)\sim q(\cdot|s,a)\in P(R)$ is a random variable representing the reward
obtained when action $a$ is taken in state $s$.

The MDP controlled by a policy $\mu$ induces a Markov chain $M^{\mu}$ with
reward distribution $q^{\mu}(\cdot|s)=q(\cdot|s,\mu(s))$ such that
$R^{\mu}(s)=R(s,\mu(s))\in q^{\mu}(\textperiodcentered|s)$, transition kernel
$P^{\mu}(\cdot|s)=P(\cdot|s,\mu(s))$, and stationary distribution over states
$\pi^{\mu}$. A discount factor $\gamma\in[0,1]$ determines the exponential
devaluation rate of delayed rewards.

The discounted return of a state $s$ is defined as the sum of discounted
rewards that the agent encounters when starting in state $s$ and following
some policy $\mu$ afterwards.

$D^{\mu}(s)=\sum_{t=0}^{\infty}\gamma^{t}R(Z_{t})|Z_{0}=(s,\mu(\cdot|s)),\\\
with\;S_{t+1}\sim P^{\mu}(\cdot|S_{t})$ (1)

The expected value of $D^{\mu}$ gives us the value function of the policy. The
value of a state $s$ under some policy $\mu$ can be represented in terms of
its immediate reward and the values of its successor states using the Bellman
equation as

$V^{\mu}(s)=R^{\mu}(s)+\gamma\int_{S}P^{\mu}(s^{\prime}|s)V^{\mu}(s^{\prime})ds^{\prime}$
(2)

###  II-B Reinforcement Learning

Reinforcement learning (RL) is a field of machine learning that aims at
learning suitable actions to maximize the reward in a particular situation. In
RL the transition and reward functions are unknown. Hence optimal policies are
learned from experience. These are defined as sequences of transitions
$(s_{t},a_{t},r_{t+1},s_{t+1},a_{t+1}),t=0,1,\cdots$ broken up into episodes.
States and rewards are sampled from the environment, and actions follow some
exploration policy $\pi$.

Given an estimate $G_{t}$ of the expected return at time $t$ starting from
state $s$ and taking action $a$, temporal difference learning updates the
expected return as follows:

$Q_{t+1}(s,a)=Q_{t}(s,a)+\eta_{t}(G_{t}-Q_{t}(s,a))$ (3)

where $\eta_{t}\in(0,1]$ is a learning rate parameter and varies by the
application domain. $G_{t}$t is generally bootstrapped from the current
Q-values. Popular bootstrapping algorithms include Q learning and Expected
State–action–reward–state–action (Expected SARSA) .

Q learning is an off-policy algorithm. It evaluates and improves a different
policy than the one being used to select actions. Under certain conditions Q
learning has been shown to converge to the optimal policy with a probability
of 1 . $G_{t}$ is bootstrapped using Q learning as follows:

$G_{t}^{Q}=r_{t+1}+\gamma\max_{a^{\prime}\in A}Q_{t}(s_{t+1},a^{\prime})$ (4)

where $\gamma\in(0,1)$ is a discount factor.

Expected SARSA is an on-policy algorithm. It evaluates and improves the policy
currently being used for action selection. In expected SARSA the uncertainty
of the next action $a_{t+1}$ is averaged out with respect to the policy $\pi$,
resulting in a significant reduction in variance . $G_{t}$ using expected
SARSA is given by

$G_{t}^{ExpSARSA}=r_{t+1}+\gamma
E_{a^{\prime}\sim\pi}[Q_{t}(s_{t+1},a^{\prime})]$ (5)

###  II-C Bayes-Adaptive Markov Decision Process

The Bayes-Adaptive Markov Decision Process (BAMDP) is an extension of the
conventional MDP model. The state space of the BAMDP combines the initial set
of states $S$, with the posterior parameters on the transition function. This
joint space is called a hyper-state. Transitions between hyper-states are
captured in the BAMDP transition model.

A BAMDP $M$ is a tuple
$(S^{\prime},A^{\prime},P^{\prime},P^{\prime}_{0},R^{\prime})$ where

  * •

$S^{\prime}$ is the set of hyper-states,

  * •

$A^{\prime}$ is the set of actions,

  * •

$P^{\prime}(\cdot|s,\phi,a)$ is the transition function between hyper-states,
when action $a$ is taken in hyper-state $(s,\phi)$,

  * •

$P^{\prime}_{0}\in P(S\texttimes\phi)$ combines the initial distribution over
states with the prior over transition functions, and

  * •

$R^{\prime}(s,\phi,a)=R(s,a)$ represents the reward obtained when action $a$
is taken in state $s$.

The value function of the BAMDP can be expressed using the Bellman equation:

$V^{*}(s,\phi)=\max_{a\in
A}[R^{\prime}(s,\phi,a)+\gamma\sum_{(s^{\prime},\phi^{\prime})}{P^{\prime}(s^{\prime},\phi^{\prime}|a,\phi,a)V^{*}(s^{\prime},\phi^{\prime})}]\\\
=\max_{a\in A}[R(s,a)+\gamma\sum_{s^{\prime}\in
S}{\frac{\phi^{a}_{s,s^{\prime}}}{s^{n}\in
S^{\phi^{a}_{s,s^{\prime\prime}}}}V^{*}(s^{\prime},\phi^{\prime})}]$ (6)

The posterior over the transition function is represented by a Dirichlet
distribution which is formulated based on prior knowledge of the application
domain.

###  II-D Exploration vs Exploitation

A major issue in reinforcement learning is exploration, i.e., how the agent
should choose actions while learning about a task. Exploitation, on the other
hand, looks at which actions are to be selected in order to maximize the
expected reward with respect to the current value function estimate.
Exploration here uses epsilon greedy policies defined as

$\pi_{t}(s,a)=\begin{cases}1-\epsilon+\frac{\epsilon}{|A|}&\text{if a =
$argmax_{a^{\prime}Q_{t}(s,a^{\prime})}$}\\\
\frac{\epsilon}{|A|}&\text{otherwise}\end{cases}$ (7)

Either a random action is chosen from the action space $A$ for exploration,
with some probability $\epsilon_{t}\in[0,1]$, or a greedy action is selected
for exploitation according to $Q_{t}$. The optimal value for $\epsilon$ is
generally problem-specific, and determined empirically. $\epsilon$ is
typically annealed in order to favor exploration in the beginning, and
exploitation closer to convergence [Sutton and Barto, 2018].

###  II-E Bayesian Learning

Bayesian learning involves making inferences like point and interval estimates
about some random variable $X$ from a probability distribution over $X$. The
process of inferring $X$ from samples of some random variable $Y$ involves the
following general steps:

  1. 1.

Choose a probability density $P(X)$, called the prior distribution, that
represents beliefs about the random variable $X$ based on prior domain
knowledge

  2. 2.

Select a statistical model $P(Y|X)$ that represents a statistical dependence
between $X$ and $Y$

  3. 3.

Gather data points for $Y$

  4. 4.

Update beliefs about $X$ by computing the posterior distribution using Bayes
rule

$P(X|Y=y)=\frac{P(y|X)P(X)}{\int P(y|X^{\prime})P(X^{\prime})dX^{\prime}}$ (8)

###  II-F Bayesian Model Combination

1:procedure
$\epsilon$-BMC($\mu_{0},\tau_{0},a_{0},b_{0},\hat{\mu}=0,\hat{\sigma}^{2}=\infty,\infty,\beta$)
Initialize $s$ at root

2: for each episode do

3: for each step in the episode do

4: $\epsilon\leftarrow\frac{\alpha}{\alpha+\beta}$

5: choose action a using $\epsilon$-greedy policy $\pi$

6: take action a, observe r and s’

7: $G^{Q}\leftarrow r+\gamma\max_{a^{\prime}}Q(s^{\prime},a^{\prime})$

8: $G^{U}\leftarrow
r+\gamma\frac{1}{|A|}\sum_{a^{\prime}}Q(s^{\prime},a^{\prime})$

9: $G^{E}xpS\leftarrow
r+\gamma\sum_{a^{\prime}}\pi^{\epsilon}(a^{\prime}|s^{\prime})Q(s^{\prime},a^{\prime})$

10: $Q(s,a)\leftarrow Q(s,a)+\eta[G^{E}xpS-Q(s,a)]$

11: Update $\hat{\mu}$ and $\hat{\sigma}^{2}$ using $G^{E}xpS$

12: Compute a, b, $e^{Q}$ and $e^{U}$

13: Update $\alpha$ and $\beta$

14: $s\leftarrow s^{\prime}$

15: end for

16: end for

17:end procedure

Algorithm 1 $\epsilon$ \- Bayesian Model Combination

###  II-G Active Malware Analysis

AMA aims to develop dynamic analysis systems that perform actions in order to
trigger different behaviors of the malware being analyzed. Payload deployment
can often be hidden behind conditional requirements that look for specific
user actions, especially in mobile malware. Initial approaches to AMA have
selected user inputs based on existing data about user behavior patterns or
simply used a pseudo-random selection of possible actions.

Suarez-Tangil et al. [5] propose an analyzer that uses stochastic models
extracted from samples of real recorded user behaviors to reproduce specific
conditions in order to trigger malicious behavior. Bhandari et al. [6] propose
using random triggers in a runtime semantic-aware malware detector that is
hardened against code injection. CANDYMAN [7] introduces a Markov chain model
that expresses the malicious dynamics of the malware being analyzed. The model
is generated using malware behavior observed by executing the malware in a
controlled environment. The features embedded in the Markov chains are then
extracted and used for classification. The goal is to minimize the number of
states composing the Markov chains and consequently the size of the feature
space. Action selection for triggering malware behavior is random.

The random action selection in these initial approaches limited their
effectiveness in terms of execution path exploration since the action
selection could not adapt to the malware behavior. More recent work in AMA
proposes employing more intelligent methods like reinforcement learning to
train autonomous agents for action selection. These take a game-theoretic
approach and model AMA as stochastic or Bayesian games between the malware and
the analyzer agent.

SECUR-AMA [11] treats AMA as a stochastic game and proposes a Monte Carlo Tree
Search based approach to dynamically generate the malware model at runtime.
This approach does not rely on human experts for action selection. Malware
authors often inject noise in their execution traces to defeat detection.
SECUR-AMA remains vulnerable to this issue and proposes some ways to ignore
the possible noise.

BAMA [9] models malware analysis as a Bayesian game. The analyzer uses prior
knowledge about malware families to select triggering actions that reflect the
current belief regarding which family the malware being analyzed belongs to.
The aim is to lower the uncertainty about a malware sample’s type rather than
its behavior.

##  III Model Design

This section describes the design choices made for the malware model and the
Bayesian Model Combination implementations.

###  III-A Malware Model

The action space $A$ is obtained by extracting the list of intents from the
malware manifest file. As such, the action space can potentially change for
each malware sample analyzed. The root node represents the main activity where
each android application is initialized. Further nodes on the tree represent
Android API calls made in response to trigger actions selected by the
analyzer. The final result is a call graph obtained for each malware sample
analyzed. Paths on the model graph are possible execution paths of the malware
and show the probability of reaching a particular terminal state assuming a
specific set of actions chosen by the analyzer.

##  IV Experimental Setup and Analysis Process

The goal of the analyzer is to build the Android API call graph of the malware
sample being analyzed with the minimal number of action selections required.
The choice of the actions executed by the analyzer is based on the Bayesian
Model Combination implementation running on a copy of the current malware
model, which is updated at the end of every round. The action selected by the
analyzer is executed on the running malware sample and the response, in the
form of a sequence of API calls is stored. Each trace is represented as a
path, starting at Init(), and is used to update the malware model graph and
node statistics, including the transition probabilities between consecutive
API calls. The analysis is stopped after a fixed number of trigger actions
have been executed. The model thus generated is given as the final output. A
linear SVM [12] is trained for classification on the call graphs obtained
using different number of analyzer actions and the $F_{1}$ scores are used to
determine the optimal number of analyzer actions required.

The experiments are conducted on an emulated Android 10 image instrumented
using Frida [13] to gather the API calls. The machine used to carry out the
experiments is configured with an AMD Ryzen 5 3600 6-core processor at 3.6
GHz, 32 GB of RAM and an NVIDIA RTX 2070 Super with 8GB of GDDR6 memory. The
analyzers are trained on a dataset containing 15000 malware samples collected
from Androzoo [14]. The set contains 150 sample each from 100 different
malware families.

##  V Evaluation

The average time for the analyzer to select an action and observe the response
was empirically observed to be about 11 seconds. After every observation, the
emulated Android image is reset. In the experiments the time allotted to each
malware sample is fixed at 5 minutes accounting for Android boot time, since a
fresh virtual image is used for every new malware analysis, and the action
selection and trace collection. We do not fix the size of the action space
since a new list is extracted from each new sample. The analyzer thus
considers a fresh action space for each malware and does not consider any
other actions it might have encountered with other malware samples.

TABLE I: Total Time taken per malware sample and optimal number of analyzer actions needed Analyzer Model | Analyzer Actions | Time [seconds]
---|---|---
Bayes-UCT | 8 | 168
BMC-Constant $\epsilon$ | 8 | 168
$\epsilon$-BMC | 7 | 147

Figure 1: F1-Scores for Linear SVM classifier

Table I shows the optimal number of analyzer actions for BMC with a constant
annealing approach to $\epsilon$, $\epsilon$-BMC and BayesUCT [15].
$\epsilon$-BMC gives an improvement compared to constant $\epsilon$ in terms
of the number of actions required. BMC-Constant $\epsilon$ gives similar
performance to BAMA [9] with the number of actions required being 8. The
improvement offered by $\epsilon$-BMC, while small in terms of the number of
actions required, becomes more significant when considering the time taken to
generate the model when the number of actions is limited to the optimal number
found, and scalled up to larger number of samples being analyzed. Fig. 1 shows
the $F_{1}$-score values obtained by a linear SVM trained for classification
on the call graphs generated and shows how performance is affected by the
number of analyzer actions used to obtain the graphs. This allows us to choose
the number of analyzer actions shown in Table I. The Bayesian Model
combination approach described in this paper requires significantly lower
number of analyzer actions to reach the best overall classification scores as
compared to the stochastic game and UCT based approach of SECUR-AMA, even with
different UCB algorithms [16].

##  VI Conclusion

Previous work in AMA has explored some Bayesian approaches by modeling AMA as
a Bayesian game and using Bayesian learning for action selection in dynamic
malware analysis. We consider another model for AMA and examine a Bayesian
Model Combination based approach to the problem. $\epsilon$-BMC shows
significant improvement over stochastic game models and some improvement over
other Bayesian models considered in the literature. There is scope for
improvements by considering Boltzmann exploration and comparing against other
adaptive annealing approaches.

## References

  * [1] “Mobile Operating System Market Share Worldwide”, https://gs.statcounter.com/os-market-share/mobile/worldwide, 2021
  * [2] “AV Test Institute Malware statistics”, https://www.av-test.org/en/statistics/malware, 2021
  * [3] Brian A Wichmann et al.  “Industrial perspective on static analysis”  In _Software Engineering Journal_ 10.2 IET, 1995, pp. 69–75  DOI: 10.1049/sej.1995.0010
  * [4] Andreas Moser, Christopher Kruegel and Engin Kirda  “Limits of static analysis for malware detection”  In _Twenty-Third Annual Computer Security Applications Conference (ACSAC 2007)_ , 2007, pp. 421–430  IEEE  DOI: 10.1109/ACSAC.2007.21
  * [5] Guillermo Suarez-Tangil, Mauro Conti, Juan E. Tapiador and Pedro Peris-Lopez  “Detecting Targeted Smartphone Malware with Behavior-Triggering Stochastic Models”  In _Computer Security - ESORICS 2014_ Cham: Springer International Publishing, 2014, pp. 183–201  DOI: 10.1007/978-3-319-11203-9˙11
  * [6] Shweta Bhandari et al.  “SWORD: Semantic aWare andrOid malwaRe Detector”  In _Journal of Information Security and Applications_ 42, 2018, pp. 46–56  DOI: 10.1016/j.jisa.2018.07.003
  * [7] Alejandro Martín, Víctor Rodríguez-Fernández and David Camacho  “CANDYMAN: Classifying Android malware families by modelling dynamic traces with Markov chains”  In _Engineering Applications of Artificial Intelligence_ 74, 2018, pp. 121–133  DOI: 10.1016/j.engappai.2018.06.006
  * [8] James John Martin  “Bayesian decision problems and Markov chains”  Wiley, 1967
  * [9] Riccardo Sartea, Georgios Chalkiadakis, Alessandro Farinelli and Matteo Murari  “Bayesian Active Malware Analysis”  In _Proceedings of the 19th International Conference on Autonomous Agents and MultiAgent Systems_ , AAMAS ’20  Auckland, New Zealand: International Foundation for Autonomous AgentsMultiagent Systems, 2020, pp. 1206–1214
  * [10] Richard Bellman  “A Markovian decision process”  In _Journal of mathematics and mechanics_ 6.5 JSTOR, 1957, pp. 679–684
  * [11] Riccardo Sartea, Alessandro Farinelli and Matteo Murari  “SECUR-AMA: Active Malware Analysis Based on Monte Carlo Tree Search for Android Systems”  In _Engineering Applications of Artificial Intelligence_ 87, 2020, pp. 103303  DOI: https://doi.org/10.1016/j.engappai.2019.103303
  * [12] Nir Nissim, Robert Moskovitch, Lior Rokach and Yuval Elovici  “Detecting unknown computer worm activity via support vector machines and active learning”  In _Pattern Analysis and Applications_ 15.4 Springer, 2012, pp. 459–475  DOI: 10.1007/s10044-012-0296-4
  * [13] “Frida - A world-class dynamic instrumentation framework”, https://frida.re/, 2021
  * [14] Kevin Allix, Tegawendé F. Bissyandé, Jacques Klein and Yves Le Traon  “AndroZoo: Collecting Millions of Android Apps for the Research Community”  In _Proceedings of the 13th International Conference on Mining Software Repositories_ , MSR ’16  Austin, Texas: ACM, 2016, pp. 468–471  DOI: 10.1145/2901739.2903508
  * [15] Abhilash Hota and Jürgen Schönwälder  “A Bayesian Approach to Monte-Carlo Tree Search in Active Malware Analysis” Submitted to: 37rd International Conference on Information Security and Privacy Protection (IFIP SEC 2022), 2022
  * [16] Abhilash Hota and Jürgen Schönwälder  “A Comparison of Upper Confidence Bounds for Exploration in Active Malware Analysis” Submitted to: 1st ACM Workshop on Robust Malware Analysis (WoRMA 2022), 2022


# Measuring dynamical phase transitions in time series

Bulcsú Sándor<EMAIL_ADDRESS>Department of Physics, Babes-Bolyai
University, Cluj-Napoca, Romania
Transylvanian Institute of Neuroscience, Cluj-Napoca, Romania  András Rusu
Department of Physics, Babes-Bolyai University, Cluj-Napoca, Romania
Transylvanian Institute of Neuroscience, Cluj-Napoca, Romania  Károly Dénes
Department of Physics, Babes-Bolyai University, Cluj-Napoca, Romania
Transylvanian Institute of Neuroscience, Cluj-Napoca, Romania  Mária Ercsey-
Ravasz<EMAIL_ADDRESS>Department of Physics, Babes-Bolyai
University, Cluj-Napoca, Romania
Transylvanian Institute of Neuroscience, Cluj-Napoca, Romania  Zsolt I. Lázár
<EMAIL_ADDRESS>Department of Physics, Babes-Bolyai University, Cluj-
Napoca, Romania
Transylvanian Institute of Neuroscience, Cluj-Napoca, Romania

(August 28, 2024)

###### Abstract

There is a growing interest in methods for detecting and interpreting changes
in experimental time evolution data. Based on measured time series, the
quantitative characterization of dynamical phase transitions at bifurcation
points of the underlying chaotic systems is a notoriously difficult task.
Building on prior theoretical studies that focus on the discontinuities at
$q=1$ in the order-$q$ Rényi-entropy of the trajectory space, we measure the
derivative of the spectrum. We derive within the general context of Markov
processes a computationally efficient closed-form expression for this measure.
We investigate its properties through well-known dynamical systems exploring
its scope and limitations. The proposed mathematical instrument can serve as a
predictor of dynamical phase transitions in time series.

time series analysis, dynamical phase transitions, entropy, chaos

††preprint: APS/123-QED

## Complex systems and time series—

Recently, in the era of complexity science, the significance of efficient
tools for characterizing quantitatively the temporal evolution of real-world
complex systems has surged [1, 2]. In the absence of defining equations of
motion the focus shifts towards analyzing the observed time series [3, 4, 5].
Qualitative changes in the dynamics are important from the point of view of
prediction and control. In biological time series, for example, accurately
detecting epileptic seizures in EEG recordings and identifying cardiac
fibrillation in ECG data are vital [6, 7]. Therefore, the study of critical
transitions and the identification and location of tipping points are
fundamental aspects of complex systems [8, 9, 10]. There are theoretical
methods for describing the critical behavior of dynamical systems allowing for
the definition and estimation of quantities with known analogs in
thermodynamic formalism including entropy, partition function, free energy,
and diffusion coefficient [11, 12]. These are mainly demonstrated in 1D
systems covering both phase space [13, 14] and trajectory (history) space
statistics [15, 16, 17, 18]. However, analytical relations can be exploited
only in a limited number of cases. For a numerical measurement of the entropy,
a symbolic representation of time series is realized in one of two ways:
spatial encoding (threshold-based phase space partitioning) or temporal
encoding (ordinal patterns) [19, 3, 20, 21]. As a result, the problem can be
transferred to the context of symbolic dynamics, setting the stage for a
computational and information-theoretical perspective [22, 23].

In the absence of equations defining the dynamics of the system, such as for
Markov processes or experimental time series it is a challenge to estimate the
relevant statistical quantities related to phase transitions. Here we briefly
review the role of order-$q$ Rényi-entropy in characterizing dynamical
properties as they manifest both in the phase space and trajectory space.
Based on previous studies [15, 16, 17, 18] we propose to measure dynamical
phase transitions via the derivative of the entropy at $q=1$ exhibiting
singularities at tipping points. We derive a computationally efficient closed-
form expression of this measure within the framework of a first-order Markov
process model of the dynamics. Here we present examples for one and two-
dimensional maps, however, the method can be applied in higher dimensions,
too. It can be generalized to continuous time dynamical systems and
multivariate time series by adequate time discretization, for example by
constructing the corresponding Poincare map, as described in Refs. [3, 24].
Our approach bridges a current methodological gap by addressing both the
theoretical and practical aspects of analyzing phase transitions in complex
dynamical systems, with applicability spanning biological, socio-economic,
geophysical, and climate-related time series.

## Generalized entropies.—

The generalized entropy (also known as order-$q$ Rényi entropy) for a discrete
distribution $\bm{p}=(p_{1},p_{2},\dots,p_{N})$ is defined as

$H_{q}(\bm{p})=\frac{1}{1-q}\ln\sum_{i=1}^{N}p_{i}^{q}$ (1)

where $H_{q}\geq H_{q^{\prime}}$ whenever $q^{\prime}>q$ [12]. In the $q=1$
limit, we get the mean

$H_{1}(\bm{p})=-\langle\ln\bm{p}\rangle_{\bm{p}}=-\sum_{i=1}^{N}p_{i}\ln
p_{i}\,,$ (2)

which is the Shannon entropy for the same distribution. The same limit of the
derivative of the Rényi-entropy, with respect to $q$, gives [25]

$-2H_{1}^{\prime}(\bm{p})=\Delta^{2}_{\bm{p}}\ln(\bm{p})=\langle\ln^{2}\bm{p}\rangle_{\bm{p}}-\langle\ln\bm{p}\rangle^{2}_{\bm{p}}\
,$ (3)

hence, it is proportional to the bit-number variance of the distribution,
$\Delta^{2}_{\bm{p}}\ln(\bm{p})$ [13, 12]. Depending on the context, the same
quantity is also referred to as diffusion coefficient [18] or specific heat
[14] depending on whether $\ln\bm{p}$ is associated with displacement during a
random walk or energy of a thermodynamic system.

## Dynamical systems.—

For discrete-time dynamical systems, we need an appropriate partitioning of
the phase space by defining a set of $n$ nonoverlapping cells that cover the
attractor entirely. In the limit of infinitesimal cell sizes, Eq. (1) can be
used to estimate the resolution-independent spectrum of Rényi dimensions. The
Shannon entropy and the bit-number variance can be defined using the coarse-
grained approximation
${\bm{\rho}}^{*}=(\rho_{1}^{*},\rho_{2}^{*},\dots,\rho_{n}^{*})$ of the
corresponding natural distribution on the attractor,

$C_{1}=H_{1}({\bm{\rho}}^{*})\,,\qquad
C_{2}=-2H^{\prime}_{1}({\bm{\rho}}^{*})$ (4)

where the latter exhibits singularities at bifurcation points [13, 14]. While
the natural distribution can be estimated numerically its analytical
expression is only known for some parameters of a limited number of maps. For
time series, in the absence of the map, the quantities $C_{1}$ and $C_{2}$
remain invariant under the reshuffling of the data points, so it is not
suitable for measuring the more relevant trajectory-level information [23]. To
that end a trajectory of the system can be defined as a string of symbols,
viz. coarse-grained states $\mathcal{X}(t)=\\{s_{1},s_{2},\dots,s_{t}\\}$ if
$\mathbf{x}(\tau=1)\in$ cell $s_{1}$, $\mathbf{x}(\tau=2)\in$ cell $s_{2}$,
$\dots$, $\mathbf{x}(\tau=t)\in$ cell $s_{t}$ where $\mathbf{x}(\tau)$ denotes
the state of the dynamical system at time step $\tau$. Let us denote by
$\mathcal{P}(\mathcal{X}(t))$ the joint probability of the trajectory
$\mathcal{X}(t)$ of $t$ consecutive steps. The corresponding dynamic entropy
in the probability space ${\mathbfcal{P}}(t)$ of all possible trajectories
$\mathbfcal{X}(t)$ is defined as [26, 27, 28]:

$H_{q}({\mathbfcal{P}}(t))=\frac{1}{1-q}\ln\sum_{\mathbfcal{X}(t)}\left[\mathcal{P}(\mathcal{X}(t))\right]^{q}$
(5)

While the above expression usually diverges with the length of the
trajectories $t$, the dynamical Rényi entropy rate,

$K_{q}({\mathbfcal{P}})=\lim_{t\rightarrow\infty}\frac{1}{t}H_{q}({\mathbfcal{P}}(t))\,,$
(6)

characterizes the scaling behavior in the limit of infinite time. The above
quantity generally depends on the partitioning. Using the generating partition
the supremum over partitions is reached. As such a partitioning is generally
not known, one may choose instead a grid of cells of equal size. A
sufficiently high resolution is also expected to provide a good estimation of
the supremum [29]. Alternatively, one may also use numerical or machine-
learning algorithms for estimating optimal partitions from observed data (see
e.g. [30, 31, 2]).

The Rényi-entropy at $q=1$, denoted here by $S$, is of special interest and it
is called the Kolmogorov-Sinai (KS) entropy or metric entropy [23]:

$S=K_{1}({\mathbfcal{P}})\,,\qquad\Lambda=-{2}K_{1}^{\prime}({\mathbfcal{P}})\,,$
(7)

as it measures the entropy production rate due to the chaotic nature of the
dynamics. Based on [15, 32] there exists a dynamical phase transition at
borderline cases of chaos, hence the derivative of the Rényi-entropy with
respect to $q$ at $q=1$ (interpreted as diffusion coefficient), quantified by
$\Lambda$ in Eq. (7), diverges at the critical point.

## Markov chains.—

Discrete-time autonomous dynamical systems with a continuous state space can
be approximated by time-homogeneous discrete-time Markov chains with finite
state space. This model requires two elements. Firstly, the states are
discretized by choosing a suitable partitioning as described above. Secondly,
the transition probabilities between states $i$ and $j$ are defined as a
conditional probability $w_{ij}=p(s_{t+1}=j|s_{t}=i)$, where
$i,j=1,\dots,N\leqslant n$, where $n$ is the number of cells and $N$ is the
number of symbols present in the encoded time series. Given such a Markov
chain with the corresponding transition matrix, $(\mathbf{W})_{ij}=w_{ij}$, we
call its network representation a state-transition network (STN) [24]. This
approximation becomes exact if the coarse-graining corresponds to Markov
partitioning, one that is only known exactly for a handful of systems and can
be approximated or optimized numerically in some other cases as well [19].
Additionally, the first-order Markov chain description of the dynamics can be
improved by considering higher-order states corresponding to symbolic
sequences $s_{t-m+1},\dots,s_{t}$ of length $m$ which leads to an $m$-th order
Markov process,
$p(s_{t-m+2},\dots,s_{t+1}|s_{t-m+1},\dots,s_{t})=p(s_{t+1}|s_{t-m+1},\dots,s_{t})$
in terms of the original symbolic states. The number of symbols encoding
$m$-th order states with non-zero probability during the dynamics is typically
$N\ll n^{m}$. The stationary distribution of the Markov chain, ${\bm{\rho}}$,
is the eigenvector of the transition probability matrix $\mathbf{W}$
associated with the unit eigenvalue (see Section B in the Supplemental
Material [33]). In the special case of first-order states, $m=1$, it coincides
with the coarse-grained natural distribution on the attractor,
${\bm{\rho}}={\bm{\rho}}^{*}$.

For estimating the dynamical entropies $K_{q}$, we define the trajectory
$\mathcal{X}(t)$ as a sequence of states/nodes visited during a random walk of
$t$ steps on the network, i.e. a realization of the Markov chain. The
probability of such a trajectory can be given as

$\mathcal{P}(\mathcal{X}(t))=\rho_{s_{1}}\mathcal{P}(\mathcal{X}(t)|s_{1})\,,\quad\mathcal{P}(\mathcal{X}(t)|s_{1})=\prod_{\tau=1}^{t-1}(w_{s_{\tau}s_{\tau+1}})\,.$
(8)

Using the Markov-chain description, the full spectrum of the truncated Rényi-
entropy [34, 16, 32] can be given as (see e.g. [35] and Section C in the
Supplemental Material [33]).

$\tilde{K}_{q}={\ln\alpha_{\text{max}}(q)}/\left({1-q}\right)$ (9)

where $\alpha_{\text{max}}(q)$ is the largest eigenvalue of the
$(\mathbf{W}_{q})_{ij}=w_{ij}^{q}$ matrix (element-wise $q$ power). The
truncated entropy, $\tilde{K}_{q}$, is expected to exceed the exact value,
$K_{q}$, for reasons that are straightforward to show in the $q=1$ case from
the properties of conditional entropies (see also Ref. [34]). Section G in the
Supplemental Material [33] demonstrates Eq. (9) for two prototypical one-
dimensional systems: the piece-wise linear asymmetric tent map [12] for which
the analytically known Rényi-entropy spectrum is reproduced numerically, and
in the case of the critical map [15] exhibiting intermittent dynamics the
discontinuity at $q=1$ is confirmed together with the exact value of the KS
entropy.

For memory-less dynamics the truncated estimation of the entropy rate, $S$,
and diffusion coefficient, $\Lambda$, from Eq. (7) turn into the bit-number
statistics, $C_{1}$ and $C_{2}$, from Eq. (4) (see Section C in the
Supplemental Material [33]).

In the case of the logistic map, $x_{\tau+1}=rx_{\tau}(1-x_{\tau})\,,$ the
truncated entropy exhibits a phase transition-like discontinuity at $q=1$ only
in the vicinity of critical values of the parameter $r$ corresponding to
bifurcation points (e.g. $r=3.82842$, compare the lines denoted by the square-
shaped marker in panels (a) and (c) of Fig. 1). Around these points, the
evolution of the system presents properties characteristic of phase
transitions in physical systems such as critical slowing down. In the case of
limit cycles, the full entropy spectrum is identically zero, for partially
predictable chaos [36], however, it is characterized by relatively low, quasi-
constant values (see the star- and circle-shaped markers in panel (a) of Fig
1, respectively). Finally, for $r=4$ viz for fully developed chaos, all Rényi
entropies are $K_{q}=\ln 2$ (see the dashed line in Fig 1(a) and its numerical
estimate denoted by the yellow line with the diamond-shaped marker) [12],
which is exactly reproduced by its truncated version in Eq. (9) when using the
Markov partitions [0, 0.5) and [0.5, 1) (not shown here).

The same can be done if the map is higher-dimensional. Here, we consider the
well-known two-dimensional Henon map [37],
$x_{\tau+1}=1-ax_{\tau}^{2}+y_{\tau}\,,\,y_{\tau+1}=bx_{\tau}\,,$ with $b=0.3$
and control parameter $a$. A similar phase transition can be observed in the
vicinity of the critical values of the control parameter (e. g. $a\approx
1.2265$, compare panels (b) and (d) of Fig. 1), see also in [24]. Note that
for $a=1.4$ the accepted value of the topological entropy
$h_{\mathrm{t}}\approx 0.465$ is reproduced by $\tilde{K}_{0}\approx
h_{\mathrm{t}}$ [38, 39]. Furthermore, the KS entropy is equal to the maximal
Lyapunov exponent, $\tilde{K}_{1}\approx\lambda\approx 0.42$, as expected by
Pesin’s relation [23].

The critical behaviors illustrated above can be described by the quantity
$\Lambda$ in Eq. (7) corresponding to specific heat in thermodynamics and
associated with the bit number variance of the trajectory probabilities, see
Eq. (3). As it is shown in the next section, this bit number statistics has an
intuitive interpretation and allows for a complete analytical description
within the framework of Markov chains.

Figure 1:  Matching the classical picture of the logistic map (left column)
and the Henon map (right column) to its Markov-chain-based model. For the
logistic map the [0,1] interval, while for the Henon map the phase plane
$[-2,2]\times[-2,2]$ is partitioned into $n=2^{5}$ and $n=2^{5}\times 2^{5}$
equally sized cells, respectively. (a,b) The spectrum of the truncated Rényi
entropy $\tilde{K}_{q}$ computed by using Eq. (9) and states of order $m=12$
for four different parameter values (given in the legends). The location of
markers corresponds to the KS entropy $S=\tilde{K}_{1}$ given by Eq. (11). The
parameter values $r$ and $a$ used as examples in panels (a) and (b),
respectively, are indicated by the corresponding markers and the vertical
dashed lines in the panels below. (a) For the logistic map, the thin dashed
line denotes the theoretical Rényi entropy, $K_{q}=\ln 2$, for $r=4$ [12]. (b)
For the Henon map the topological entropy $h_{\mathrm{t}}\approx
0.465\approx\tilde{K}_{0}$ and the maximal Lyapunov exponent $\lambda\approx
0.42\approx S$ are shown for $a=1.4$ for comparison [38]. (c,d) The
bifurcation diagrams of the above systems (see Sec. G in the Suppl. Mat. [33]
for a broader range of the control parameters). (e,f) The KS entropy, $S(m)$,
and (g,h) the Lyapunov measure, $\Lambda(m)$, from Eqs. (11) and (14),
respectively. Lines in increasingly darker shades of gray correspond to Markov
chains using states of order $m=1,2,4,12$, respectively. The analytical
results (gray lines) given in Eq. (14) are compared to the statistics of
$10^{3}$ simulated random walks (gray circular markers) on the corresponding
STNs in the insets for order $m=12$. In panels (e,f) the largest Lyapunov
exponent, $\lambda$, is also shown by the thin black line for comparison. Note
that $S(12)\approx\lambda$ for $\lambda>0$.

## Random-walks statistics.—

The two measures from Eq. (7) can be expressed as the mean and variance of

$\mathcal{L}(t)\equiv-\ln\mathcal{P}(\mathcal{X}(t)|s_{1})$ (10)

over all trajectories $\mathcal{X}(t)$ and initial states $s_{1}$ (see Section
D in Supplemental Material [33]). They can be estimated within the framework
of first-order Markovian processes discussed above by using the truncated
entropy $\tilde{K}_{q}$. According to Eq. (8) the quantity $\mathcal{L}(t)$
can be interpreted as the total length of a random walk trajectory,
$\mathcal{X}(t)$, made up of steps of length $l_{\tau}=-\ln
w_{s_{\tau}s_{\tau+1}}$ with $\mathcal{L}(t)=\sum_{\tau=1}^{t}l_{\tau}$.

Due to the Markov chain central limit theorem [40], the distribution of total
path lengths, $\mathbfcal{L}(t)$, over all possible realizations of random
walks, $\mathbfcal{X}(t)$, approaches a normal distribution in the limit of
$t\rightarrow\infty$, which can be fully characterized by its mean,
$\langle\mathbfcal{L}(t)\rangle_{\mathbfcal{P}}$, and variance,
$\Delta_{\mathbfcal{P}}^{2}\mathbfcal{L}(t)$ (see Section G in the
Supplemental Material [33]). These two quantities describing the statistics of
random-walk path lengths turn out to be closely related to the phase-
transition measures proposed in Eq. (7) as discussed below.

The truncated estimation of the entropy rate from Eq. (7) can be given as the
time-normalized mean path length,
$\langle\mathbfcal{L}(t)\rangle_{\mathbfcal{P}}/t$ (see Sections D and F in
Supplemental Material [33]), which reduces to the Kolmogorov-Sinai entropy
rate of the underlying Markov process [24]:

$S=\langle l\rangle=-\sum_{i,j=1}^{N}\rho_{i}w_{ij}\ln
w_{ij}={\bm{\rho}}^{\top}\mathbf{L}_{1}\mathbf{1}\,,$ (11)

where $\langle\cdot\rangle$ denotes the average over the ensemble of random
walkers, ${\bm{\rho}}^{\top}$ and $\mathbf{1}=(1,1,\dots,1)$ are the left and
right eigenvectors of the stochastic matrix $\mathbf{W}$ with corresponding
eigenvalue $\alpha_{1}=1$, and $(\mathbf{L}_{1})_{ij}=-w_{ij}\ln w_{ij}$ is
the probability-weighted edge-length matrix. In the context of STNs, this
entropy is equivalent to the average of the entropies associated with
individual nodes (see Section E in Supplemental Material [33]).

Similarly, the time-normalized variance,
$\Delta_{\mathbfcal{P}}^{2}\mathbfcal{L}(t)/t$, of random-walk path lengths on
the same STN is given by the Lyapunov measure (see Sections D and F in
Supplemental Material [33]),

$\Lambda=\sigma_{l}^{2}+2\sum_{\tau=1}^{\infty}\Bigl{(}\langle
l_{1}l_{1+\tau}\rangle-\langle l_{1}\rangle\langle
l_{1+\tau}\rangle\Bigr{)}\,,$ (12)

where the first term denotes the variance of the individual edge lengths,

$\sigma_{l}^{2}=\langle l^{2}\rangle-\langle l\rangle^{2}\,,\qquad\langle
l^{2}\rangle=\sum_{i,j=1}^{N}\rho_{i}w_{ij}\ln^{2}w_{ij}\,,$ (13)

sampled according to the Markov process, while the second term is the sum of
auto-covariances for all possible time lags $\tau$, with $l_{1}$ being the
edge length in the first step [41]. Using the transition probability matrix
$\mathbf{W}$, the Lyapunov measure can be expressed in finite form (see
Section F in Supplemental Material [33]) as

$\Lambda={\bm{\rho}}^{\top}\mathbf{{L}}_{2}\mathbf{1}-({\bm{\rho}}^{\top}\mathbf{L}_{1}\mathbf{1})^{2}+2{\bm{\rho}}^{\top}\mathbf{L}_{1}\mathbf{B}\mathbf{L}_{1}\mathbf{1}\,,$
(14)

where $(\mathbf{{L}}_{2})_{ij}=w_{ij}\ln^{2}w_{ij}$,
$\mathbf{B}\equiv(\mathbf{I}-\mathbf{W}+\mathbf{1}{\bm{\rho}}^{\top})^{-1}-\mathbf{1}{\bm{\rho}}^{\top}$.

The equivalence of the closed-form expressions of the two measures in Eqs.
(11), (14) and the corresponding path length statistics obtained from the
direct simulation of the underlying random walk is demonstrated for the
logistic and the Henon maps in the insets of panels (e), (g) and (f), (h),
respectively, of Fig. 1. As expected based on Pesin’s relation the Kolmogorov-
Sinai entropy (gray lines in Fig. 1(e), (f)) matches the Lyapunov exponent,
$\lambda$ (black line), computed for all control parameter values using
Benettin’s method [42, 43]. The two get closer as the order, $m$, of the
states increases (darker shades of gray), $S(m)\rightarrow\lambda$. However,
the Lyapunov measure tends to diverge in the neighborhood of certain
bifurcation points irrespective of the details of the partitioning and the
used order, $m$, for symbolic encoding (see the red peaks in the insets of
panels Fig. 1(g), (h). On the other hand, the width and height of the peaks
can be controlled via the applied resolution, $n$, (not shown here) and order,
$m$. Systems undergoing parameter drift are expected to produce similar
divergences during their quasi-stationary evolution prior to a dynamical phase
transition. Hence an abrupt increase in the Lyapunov measure, for example
preceding the collapse of chaos to periodic behavior, can serve as a precursor
of such tipping points.

The proposed method for measuring dynamical phase transitions expressed by
Eqs. (11) and (14) is based on common matrix algebra at a computational cost
of $\mathcal{O}(N^{3})$, and storage cost of $\mathcal{O}(N^{2})$ where $N$ is
the number of symbols. For large $N$ (high phase space resolution, $n$, or
high-order representations, $m$) the transition matrix, $\mathbf{W}$, is
typically sparse allowing for low complexity iterative numerical methods (see
Sections H and I for the iterative method and code availability, respectively,
of the Supplemental Material [33]).

## Conclusions.—

The proposed Lyapunov measure based on the first-order Markov process model of
dynamical systems can detect dynamical phase transitions manifesting as
discontinuities in the Rényi-entropy (e.g. at weak intermittency). For memory-
less dynamics, the KS entropy and the Lyapunov measure reduce to the bit-
number statistics $C_{1}$ and $C_{2}$ from Eq. (4) which rely solely on state
distribution. However, for general dynamics, the former two measures
calculated from trajectories convey additional information on the underlying
behavior of the system telling apart meaningful dynamics from random noise
with the same statistics [44, 4].

The formalism developed here allows for any symbolic representation of the
time series including ordinal partition networks [20]. The accuracy of the
numerical estimates of entropy-related quantities is solely dependent on the
intrinsic limitations of the employed representation, compare, for example,
the results on topological entropy using the ordinal patterns representations
of the Henon-map [39]. The proposed method can be employed for characterizing
dynamical systems and time series even at high-resolution coarse graining and
high-order symbolic representation allowing cost-effective numerical
estimation of the two statistical measures. Our transition matrix-based
procedure can be applied in higher dimensions and it can be generalized to
continuous time dynamical systems and multivariate time series.

In light of the recent surge in experimental data, there is a growing need for
tools that can extract useful information. We have outlined the mathematical
foundations of a method that can serve as an effective instrument in this
regard, enabling the measurement of dynamical phase transitions in time series
and facilitating the detection of early warning signals of critical
transitions in chaotic dynamics with immediate applications in areas such as
finances, climate research and neuroscience.

###### Acknowledgements.

This work was supported by the grant of the Romanian Ministry of Research,
Innovation and Digitization, CNCS - UEFISCDI, project number PN-
III-P4-PCE-2021-0408 (BS, KD, MER, ZL), PN-III-P4-ID-PCE-2020-0647 (BS, ZL),
ERANET-FLAG-ERA-ModelDXConsciousness (BS, MER, ZL), ERANET-NEURON-2-RESIST-D
(MER, BS, AR, ZL), ERANET-NEURON-2-IBRAA (MER), ERANET-FLAG-ERA-JTC2023-MONAD
(MER), ERANET-NEURON-UnscrAMBLY (MER) within PNCDI III, and SRG-UBB
32993/23.06.2023 (BS) within UBB Starting Research Grants of the Babes-Bolyai
University, and by the Collegium Talentum Programme of Hungary (AR). The
authors also acknowledge the useful discussion with Tamás Tél.

## References

  * Crutchfield [2012] J. P. Crutchfield, Between order and chaos, Nature Physics 8, 17 (2012).
  * Murphy and Bassett [2024] K. A. Murphy and D. S. Bassett, Machine-learning optimized measurements of chaotic dynamical systems via the information bottleneck, Physical Review Letters 132, 197201 (2024).
  * Sakellariou  _et al._ [2019] K. Sakellariou, T. Stemler, and M. Small, Markov modeling via ordinal partitions: An alternative paradigm for network-based time-series analysis, Physical Review E 100, 62307 (2019).
  * Boaretto  _et al._ [2021] B. R. R. Boaretto, R. C. Budzinski, K. L. Rossi, T. L. Prado, S. R. Lopes, and C. Masoller, Discriminating chaotic and stochastic time series using permutation entropy and artificial neural networks, Scientific Reports 11, 15789 (2021).
  * Bartlett  _et al._ [2022] S. Bartlett, J. Li, L. Gu, L. Sinapayen, S. Fan, V. Natraj, J. H. Jiang, D. Crisp, and Y. L. Yung, Assessing planetary complexity and potential agnostic biosignatures using epsilon machines, Nature Astronomy 6, 387 (2022).
  * Lehnertz  _et al._ [2023] K. Lehnertz, T. Bröhl, and R. von Wrede, Epileptic-network-based prediction and control of seizures in humans, Neurobiology of Disease 181, 106098 (2023).
  * Toker  _et al._ [2020] D. Toker, F. T. Sommer, and M. D’Esposito, A simple method for detecting chaos in nature, Communications Biology 3, 11 (2020).
  * Sándor  _et al._ [2013] B. Sándor, F. Járai-Szabó, T. Tél, and Z. Néda, Chaos on the conveyor belt, Phys. Rev. E 87, 42920 (2013).
  * Köglmayr and Räth [2024] D. Köglmayr and C. Räth, Extrapolating tipping points and simulating non-stationary dynamics of complex systems using efficient machine learning, Scientific Reports 14, 10.1038/s41598-023-50726-9 (2024).
  * Morr and Boers [2024] A. Morr and N. Boers, Detection of approaching critical transitions in natural systems driven by red noise, Physical Review X 14, 21037 (2024).
  * Tél [1988] T. Tél, Fractals, multifractals, and thermodynamics, An Introductory Review 43, 1154 (1988).
  * Beck and Schögl [1993] C. Beck and F. Schögl,  _Thermodynamics of Chaotic Systems: An Introduction_, Cambridge Nonlinear Science Series (Cambridge University Press, 1993).
  * Schlögl [1974] F. Schlögl, Specific heat as a general statistical measure, Zeitschrift für Physik 267, 77 (1974).
  * Schlögl and Schöll [1988] F. Schlögl and E. Schöll, Generalized specific heat as a characteristic measure in chaos, Zeitschrift für Physik B Condensed Matter 71, 231 (1988).
  * Szépfalusy  _et al._ [1987] P. Szépfalusy, T. Tél, A. Csordás, and Z. Kovács, Phase transitions associated with dynamical properties of chaotic systems, Physical Review A 36, 3525 (1987).
  * P Szépfalusy [1989] P Szépfalusy, Characterization of Chaos and Complexity by Properties of Dynamical Entropies, Physica Scripta 1989, 226 (1989).
  * Csordás and Szépfalusy [1989] A. Csordás and P. Szépfalusy, Singularities in rényi information as phase transitions in chaotic states, Phys. Rev. A 39, 4767 (1989).
  * Fujisaka [1983] H. Fujisaka, Statistical dynamics generated by fluctuations of local lyapunov exponents, Progress of Theoretical Physics 70, 1264–1275 (1983).
  * Rubido  _et al._ [2018] N. Rubido, C. Grebogi, and M. S. Baptista, Entropy-based generating Markov partitions for complex systems, Chaos 28, 10.1063/1.5002097 (2018).
  * Zou  _et al._ [2019] Y. Zou, R. V. Donner, N. Marwan, J. F. Donges, and J. Kurths, Complex network approaches to nonlinear time series analysis, Physics Reports 787, 1 (2019).
  * Zanin and Olivares [2021] M. Zanin and F. Olivares, Ordinal patterns-based methodologies for distinguishing chaos from noise in discrete time series, Communications Physics 4, 190 (2021).
  * Crutchfield and Young [1990] J. P. Crutchfield and K. Young, Computation at the onset of chaos, Entropy, complexity, and the physics of information 8, 223 (1990).
  * Boffetta  _et al._ [2002] G. Boffetta, M. Cencini, M. Falcioni, and A. Vulpiani, Predictability: a way to characterize complexity, Physics Reports 356, 367 (2002).
  * Sándor  _et al._ [2021] B. Sándor, B. Schneider, Z. I. Lázár, and M. Ercsey-Ravasz, A Novel Measure Inspired by Lyapunov Exponents for the Characterization of Dynamics in State-Transition Networks, Entropy 23, 1 (2021).
  * Buryak and Mishura [2021] F. Buryak and Y. Mishura, Convexity and robustness of the Rényi entropy, Modern Stochastics: Theory and Applications 8, 387 (2021).
  * Grassberger and Procaccia [1983a] P. Grassberger and I. Procaccia, Measuring the strangeness of strange attractors, Physica D: Nonlinear Phenomena 9, 189 (1983a).
  * Grassberger and Procaccia [1983b] P. Grassberger and I. Procaccia, Characterization of strange attractors, Physical Review Letters 50, 346 (1983b).
  * Grassberger and Procaccia [1984] P. Grassberger and I. Procaccia, Dimensions and entropies of strange attractors from a fluctuating dynamics approach, Physica D: Nonlinear Phenomena 13, 34 (1984).
  * Cohen and Procaccia [1985] A. Cohen and I. Procaccia, Computing the kolmogorov entropy from time signals of dissipative and conservative dynamical systems, Physical Review A 31, 1872 (1985).
  * Kennel and Buhl [2003] M. B. Kennel and M. Buhl, Estimating good discrete partitions from observed data: Symbolic false nearest neighbors, Physical Review Letters 91, 84102 (2003).
  * Hirata  _et al._ [2004] Y. Hirata, K. Judd, and D. Kilminster, Estimating a generating partition from observed time series: Symbolic shadowing, Physical Review E 70, 16215 (2004).
  * Kaufmann and Szépfalusy [1989] Z. Kaufmann and P. Szépfalusy, Properties of the entropies at weak intermittent states of Lorenz-type systems, Physical Review A 40, 2615 (1989).
  * [33] See Supplemental Material at [URL will be inserted by publisher] for mathematical derivations, for figures with broader range of the control parameters, further examples and code availability.
  * Szépfalusy and Györgyi [1986] P. Szépfalusy and G. Györgyi, Entropy decay as a measure of stochasticity in chaotic systems, Physical Review A 33, 2852 (1986).
  * Wu  _et al._ [2023] C. Wu, Y. Li, L. Xu, and G. Han, Rényi Entropy Rate of Stationary Ergodic Processes, IEEE Transactions on Information Theory , 1 (2023).
  * Wernecke  _et al._ [2017] H. Wernecke, B. Sándor, and C. Gros, How to test for partially predictable chaos, Scientific Reports 7, 1087 (2017).
  * All [1996] Chaos in two-dimensional maps, in  _Chaos: An Introduction to Dynamical Systems_ (Springer New York, New York, NY, 1996) pp. 193–230.
  * D’Alessandro  _et al._ [1990] G. D’Alessandro, P. Grassberger, S. Isola, and A. Politi, On the topology of the Henon map, Journal of Physics A: Mathematical and General 23, 5285 (1990).
  * Sakellariou  _et al._ [2021] K. Sakellariou, T. Stemler, and M. Small, Estimating topological entropy using ordinal partition networks, Phys. Rev. E 103, 022214 (2021).
  * Jones [2004] G. L. Jones, On the markov chain central limit theorem, Probability Surveys 1, 299 (2004).
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of Markov Chain Monte Carlo_ (CRC press, 2011).
  * Benettin  _et al._ [1976] G. Benettin, L. Galgani, and J.-M. Strelcyn, Kolmogorov entropy and numerical experiments, Physical Review A 14, 2338 (1976).
  * Datseris [2018] G. Datseris, DynamicalSystems.jl: A Julia software library for chaos and nonlinear dynamics, Journal of Open Source Software 3, 598 (2018).
  * Rosso  _et al._ [2007] O. A. Rosso, H. A. Larrondo, M. T. Martin, A. Plastino, and M. A. Fuentes, Distinguishing noise from chaos, Physical Review Letters 99, 154102 (2007).


[a,c]Qun Wang

# Global polarization and spin alignment in heavy-ion collisions: past,
present and future

Wen-Bo Dong  Xin-Li Sheng  Yi-Liang Yin

###### Abstract

We give a brief overview on global polarization and spin alignment in heavy
ion collisions. The current theoretical understandings on global polarization
of hyperons and the global spin alignment of vector mesons are summarized.

##  1 Introduction

Non-central heavy-ion collisions induce a huge global orbital angular momentum
(OAM) perpendicular to the reaction plane. Due to the spin-orbit coupling, the
initial OAM is partially converted to the spin angular momentum of particles
in the quark-gluon plasma, in a similar way as the Barnett effect [1] and the
Einstein-de Hass effect [2]. This in turn results in phenomena such as
polarizations of $\Lambda$ hyperons [3, 4, 5] and vector mesons [6, 7, 8, 9],
which have been measured in recent experiments. In this proceeding, we
reviewed recent progresses on the polarization phenomena in heavy-ion
collisions and discussed our present understanding about the vector meson’s
polarization. We also proposed questions for possible future studies. One can
also read [10, 11, 12, 13, 14, 15] for recent reviews on polarizations in
heavy-ion collisions.

##  2 Global polarization

The idea was first proposed by Z.-T. Liang and X.-N. Wang in 2004 that the
huge OAM in non-central heavy-ion collisions can polarize quarks due to the
spin-orbit coupling [16, 17]. They studied the parton scatterings within the
screened potential model. Under the small angle approximation, the cross
section can be decomposed into a spin-dependent part and a spin-independent
part, while the former part is proportional to the spin-orbit coupling ${\bf
n}\cdot({\bf x}_{T}\times{\bf p})$, where $\bf n$ is the spin quantization
direction, ${\bf x}_{T}$ is the impact parameter relative to the center of the
static potential, and ${\bf p}$ is the momentum of incoming parton. Such a
coupling indicates an unpolarized quark will be polarized along the direction
of local OAM, i.e., the direction of ${\bf x}_{T}\times{\bf p}$, after
scattering with the static potential. In non-central heavy-ion collisions, the
initial geometry leads to the fact that the average longitudinal momentum per
parton depends on the position in $x$-direction, described by a nonzero
$dp_{z}/dx$. Microscopically, the local OAM in parton collisions, given by
$dp_{z}dx$, has a preferred direction along $-y$, which is the same as the
direction of the global OAM. Numerical similations based on Woods-Saxon model
and hard-sphere model [17, 18] show that the global OAM in Au+Au collisions at
200 GeV can reach $10^{5}\hbar$, while the local OAM can be as large as $1\sim
4\,\hbar$, leading to the global quark polarization perpendicular to the
reaction plane.

Instead of treating one quark as an effective potential, the polarized quark
scatterings are studied in [18] in a more realistic way by treating two quarks
as plane waves and consider the scattering through one-gluon exchange. In
order to recover the impact parameter, they made a two dimensional Fourier
transformation with respect to the transverse momentum transfer. The result
shows that the spin-dependent part of the cross section is proportional to the
spin-orbit coupling, which qualitatively agrees with the result obtained by
the potential model [17]. The model was further improved in Ref. [19] by
parameterizing the incident particles as wave packets with certain momenta and
center positions. The initial nonvanishing OAM is included by a finite
transverse distance (impact parameter) between centers of two wave packets.
For simplicity, the outgoing particles are treated as plane waves with zero
OAM. Therefore the initial OAM is fully converted to the spin polarization of
outgoing particles. The authors in [19] also imposed the causality condition
in the center of mass frame for each parton collision, which declared that the
collision happens at the same time, but the incident particles are displaced
by the impact parameter. The cross section contains a 16-dimensional integral,
which is challenging for numerical calculations. The authors developed a
Monte-Carlo integration package running on multi-GPUs [20], which makes the
integral possible within acceptable computing time. Another challenge in the
model is that the polarized squared amplitude contains more than 5000 terms,
making it hard to see the physics behind the formula. Fortunately, the
polarization production rate has a concise form as $d^{4}{\bf
P}_{q}/d^{4}x\propto\boldsymbol{\nabla}\times(\beta{\bf u})$ when the partons’
momentum is in local equilibrium while their spin is not. Here $\beta$ is the
inverse of temperature and $\bf u$ is the local fluid velocity. This work [19]
therefore provides a microscopic description for the generation of global spin
polarization, which is based on a comprehensive first principle calculation
including all possible parton scattering channels.

The kinetics of spin 1/2 particles can be analytically described by the Wigner
function approach [21, 22, 23, 24, 25, 26, 27, 28, 29, 30], see Ref. [31, 12,
15] for recent reviews. The axial vector component of the Wigner function
denotes the quark polarization density. In general, the polarization at local
equilibrium contains various sources at linear order in gradient, including
the thermal vorticity, thermal shear, fluid acceleration, and electromagnetic
fields [32, 33]. Among these sources, the thermal vorticity is dominant for
$\Lambda$’s global polarization, while the thermal shear plays the most
important role for the $\Lambda$’s longitudinal polarization [34, 35, 36, 37].
However, how the quark-gluon plasma reaches local equilibrium with spin still
remains a puzzle. For massive fermions, the kinetic theory becomes messy
because of increasing degrees of freedom (DOF). Scalar particles are described
by one Boltzmann equation, while chiral fermions are described by the chiral
kinetic theory, containing two equations for left and right-handed fermions,
respectively. However, the kinetics for massive fermions are described by the
Kadanoff-Baym equation [38, 39] for the Wigner function, which contains 16
DOF. Fortunately, the semi-classical expansion approach provides a systematic
way to derive the Wigner function order by order, which reduces the
independent DOF to 4 at each order. In this case, particles are characterized
by the matrix-valued spin dependent distributions (MVSD),

$f^{q}_{rs}(x,p)\equiv\int\frac{d^{4}q}{2(2\pi)^{3}}e^{-iq\cdot
x}\delta(p\cdot q)\left\langle a^{\dagger}(s,p_{2})a(r,p_{1})\right\rangle\,,$
(1)

where $a^{\dagger},a$ are creation and annihilation operators, respectively,
and $p_{1,2}\equiv p\pm q/2$. The trace in spin space of the MVSD is
interpreted as the unpolarized distributions, while its projection onto Pauli
matrices are related to the spin polarization as

$f^{q}_{rs}(x,p)=\frac{1}{2}f_{q}(x,p)\left[\delta_{rs}-P^{q}_{\mu}(x,p)n_{j}^{(+)\mu}(x,p)\tau_{rs}^{j}\right]\,,$
(2)

where $P^{q}_{\mu}$ is the polarization vector, $n_{j}^{(+)\mu}$ is the
Lorentz transformed four-vector of the basis direction ${\bf n}_{j}$ in the
rest frame of particle, and $\tau^{j}$ is the Pauli matrix corresponding to
the direction ${\bf n}_{j}$. In terms of the MVSD, the spin Boltzmann
equations (SBEs) can be derived from the Kadanoff-Baym equation. At zeroth
order in $\hbar$, the SBEs read [30]

$\displaystyle\frac{1}{E_{p}}p\cdot\partial_{x}\text{tr}\left[f^{(0)}(x,p)\right]$
$\displaystyle=$ $\displaystyle\mathscr{C}_{\text{scalar}}[f^{(0)}]\,,$
$\displaystyle\frac{1}{E_{p}}p\cdot\partial_{x}\text{tr}\left[n_{j}^{(+)\mu}\tau_{j}f^{(0)}(x,p)\right]$
$\displaystyle=$ $\displaystyle\mathscr{C}_{\text{pol}}[f^{(0)}]\,,$ (3)

in which only leading order MVSD $f^{(0)}$ appears. Explicit expressions for
collision terms $\mathscr{C}_{\text{scalar}}$ and $\mathscr{C}_{\text{pol}}$
can be found in Ref. [40]. At first order in $\hbar$, the SBEs are written as
[40]

$\displaystyle\frac{1}{E_{p}}p\cdot\partial_{x}\text{tr}\left[f^{(1)}(x,p)\right]$
$\displaystyle=$ $\displaystyle\mathscr{C}_{\text{scalar}}[f^{(1)},\partial
f^{(0)},f^{(0)}]\,,$
$\displaystyle\frac{1}{E_{p}}p\cdot\partial_{x}\text{tr}\left[n_{j}^{(+)\mu}\tau_{j}f^{(1)}(x,p)\right]$
$\displaystyle=$ $\displaystyle\mathscr{C}_{\text{pol}}[f^{(1)},\partial
f^{(0)},f^{(0)}]\,.$ (4)

The collision terms on the right-hand side in general depend not only on the
MVSD, but also on the space-time derivate of MVSD, $\partial^{\mu}f^{(0)}$.
Microscopically, this reflects the spin-orbit coupling during particle
scatterings, which is the source of spin polarization. Equations 3 and 4
provide a possible way to numerically simulate the evolution of spin
polarization. One can also express the kinetic equations in terms of the
distribution function in extended phase space by treating spin as an
additional phase space parameter [28, 41], or in terms of the axial current
[29].

##  3 Spin alignment

For vector mesons, studies usually focus on the spin alignment instead of the
spin polarization. This is because vector mesons mainly decay through strong
decays or dilepton decays, which preserve the parity symmetry and make it
difficult to measure the spin polarization. One of the experimental
observables for vector meson is the $00$-element of the spin density matrix,
i.e. $\rho_{00}$, known as the spin alignment of vector mesons. As spin-1
particles, vector meson’s spin can be $0$ and $\pm 1$, thus $\rho_{00}=1/3$ if
the spin does not have any preferred direction. A positive (negative)
deviation from 1/3 means the vector meson’s spin vector is preferred to align
parallel (perpendicular) to the direction of spin quantization direction. It
was first proposed by Z.-T. Liang and X.-N. Wang in 2004 that the spin
alignment is related to the polarization of constituent quark and antiquark
[16]. Since quarks can be globally polarized by the OAM in non-central
collisions, it is expected that the spin alignment measured along the
direction of global OAM (which is referred as the global spin alignment) also
have nonvanishing deviation from 1/3. Recently, the STAR collaboration
measured the spin alignment of $\phi$ and $K^{*0}$ mesons at RHIC energies
[8]. The result is shown in Fig. (1), where the spin alignment of $K^{*0}$ is
consistent with 1/3, while the spin alignment of $\phi$ shows a significant
positive derivation from 1/3. Experiments at LHC also observed global OAM for
$\phi$ , $K^{*0}$, and $J/\psi$ in Pb-Pb collisions at
$\sqrt{s_{\text{NN}}}=2.76$ TeV [6, 7, 9].

Figure 1: The STAR’s measurement on vector mesons’ $\rho_{00}$ with respect to
the second-order event plane in heavy ion collision [8] The stars and circles
correspond to $\phi$ and $K^{*0}$ mesons, respectively. The red solid line is
the fit curve using the model in Ref. [42].

It was recently found that the ordinary mechanism for quark spin polarization
can not explain the significant spin alignment for $\phi$ meson [42]. For
example, the thermal vorticity field in heavy-ion collisions is of order
$10^{-2}$, therefore a quick estimation shows that the induced spin alignment,
which is $\propto\omega^{2}$, is $\sim 10^{-4}$ and is much smaller than
$\rho_{00}$ observed in experiments. A new mechanism of strong force field
(effective vector field) effect is then introduced in [42, 43, 44, 40] to
solve this puzzle. This field refers to the strong interaction between
constituent quark (antiquark) inside the meson and the quark (antiquark) in
the surrounding medium. The field is dominated by fluctuations and therefore
has vanishing mean value. However, it has finite contribution to the spin
alignment of a flavorless vector meson, because the spin polarization for
constituent quarks and antiquarks can have strong correlation. The interplay
between quark and the strong force field can be effectively described by the
chiral quark model, with the Lagrangian having $SU_{L}(3)\times SU_{R}(3)$
symmetry [45, 46]

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\overline{\psi}\left[i\gamma_{\mu}\left(\partial^{\mu}+igG^{\mu}\right)+g_{V}\gamma_{\mu}V^{\mu}\right]\psi+g_{A}\overline{\psi}\gamma_{\mu}A^{\mu}\psi$
(5)
$\displaystyle+\frac{1}{4}f^{2}\text{Tr}\left(\partial^{\mu}\Sigma^{\dagger}\partial_{\mu}\Sigma\right)-\frac{1}{2}\text{Tr}F_{\mu\nu}F^{\mu\nu}\,,$

where $\psi$ is the quark field and $V^{\mu}$ is the effective vector fields
induced by quark currents. The interaction Lagrangian for an $s/\overline{s}$
pair coupled with the $\phi$ field is
$g_{V}\overline{s}\gamma^{\mu}sV^{\phi}_{\mu}$. The vector $\phi$ field
$V_{\phi}^{\mu}$ can polarize $s/\overline{s}$ though the magnetic or electric
components of the field strength tensor
$F_{\rho\sigma}^{\phi}=\partial_{\rho}V_{\sigma}^{\phi}-\partial_{\sigma}V_{\rho}^{\phi}$,
in a similar way as the classical electromagnetic field. The corresponding
$\phi$ meson’s spin alignment in the rest frame is given by [43]

$\displaystyle\rho_{00}\left(x,\mathbf{0}\right)$ $\displaystyle\approx$
$\displaystyle\frac{1}{3}+C_{1}\left[\frac{1}{3}\omega^{\prime}\cdot\omega^{\prime}-\left(\epsilon_{0}\cdot\omega^{\prime}\right)^{2}\right]$
(6)
$\displaystyle+C_{2}\left[\frac{1}{3}\varepsilon^{\prime}\cdot\varepsilon^{\prime}-\left(\epsilon_{0}\cdot\varepsilon^{\prime}\right)^{2}\right]$
$\displaystyle-\frac{4g_{\phi}^{2}}{m_{\phi}^{2}T_{\text{eff}}^{2}}C_{1}\left[\frac{1}{3}\mathbf{B}_{\phi}^{\prime}\cdot\mathbf{B}_{\phi}^{\prime}-\left(\epsilon_{0}\cdot\mathbf{B}_{\phi}^{\prime}\right)^{2}\right]$
$\displaystyle-\frac{4g_{\phi}^{2}}{m_{\phi}^{2}T_{\text{eff}}^{2}}C_{2}\left[\frac{1}{3}\mathbf{E}_{\phi}^{\prime}\cdot\mathbf{E}_{\phi}^{\prime}-\left(\epsilon_{0}\cdot\mathbf{E}_{\phi}^{\prime}\right)^{2}\right]\,,$

where $\epsilon_{0}$ denotes the spin quantization direction,
$\mathbf{E}_{\phi}$ and $\mathbf{B}_{\phi}$ denote electric and magnetic part
of $\phi$ field strength tensor $F_{\phi}^{\mu\nu}$, respectively.
Contributions from thermal vorticity $\omega$ and acceleration
$\boldsymbol{\varepsilon}$ are also listed in Eq. (6) for a comparison. Here
all fields with primes are defined in the meson’s rest frame. The coefficients
$C_{1}$ and $C_{2}$ are constants that depend only on the quark mass and the
$\phi$ meson’s mass [43]

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\frac{8m_{s}^{4}+16m_{s}^{2}m_{\phi}^{2}+3m_{\phi}^{4}}{120m_{s}^{2}\left(m_{\phi}^{2}+2m_{s}^{2}\right)}\,,$
$\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle\frac{8m_{s}^{4}-14m_{s}^{2}m_{\phi}^{2}+3m_{\phi}^{4}}{120m_{s}^{2}\left(m_{\phi}^{2}+2m_{s}^{2}\right)}\,.$
(7)

In order to compare with the experimental data, the results in Eq. (6) should
be transformed back to the lab frame. Then the spin alignment in the lab frame
are expressed in terms of fields in the lab frame, while the transformation
depends on momenta of $\phi$ mesons. By taking average in the phase space, the
mean value of spin alignment in the lab frame reads [44]

$\displaystyle\left\langle\rho_{00}^{\phi}\left(x,\mathbf{p}\right)\right\rangle_{x,\mathbf{p}}$
$\displaystyle\approx$
$\displaystyle\frac{1}{3}+\frac{1}{3}\underset{i=1,2,3}{\Sigma}\left\langle
I_{B,i}(\mathbf{p})\right\rangle_{p}\frac{1}{m_{\phi}^{2}}\left[\left\langle\omega_{i}^{2}\right\rangle_{x}-\frac{4g_{\phi}^{2}}{m_{\phi}^{2}T_{\text{eff}}^{2}}\left\langle\left(\mathbf{B}_{i}^{\phi}\right)^{2}\right\rangle_{x}\right]$
(8) $\displaystyle+\frac{1}{3}\underset{i=1,2,3}{\Sigma}\left\langle
I_{E,i}(\mathbf{p})\right\rangle_{p}\frac{1}{m_{\phi}^{2}}\left[\left\langle\varepsilon_{i}^{2}\right\rangle_{x}-\frac{4g_{\phi}^{2}}{m_{\phi}^{2}T_{\text{eff}}^{2}}\left\langle\left(\mathbf{E}_{i}^{\phi}\right)^{2}\right\rangle_{x}\right]\,,$

where $\left\langle\cdots\right\rangle_{x}$,
$\left\langle\cdots\right\rangle_{p}$ denote average in space-time and in
momentum space, respectively. The momentum-dependent functions $I_{E/B,i}(p)$
can be found in [44]. Equation (8) has a perfect factorization of $x$ and
$\mathbf{p}$ dependence. Meanwhile, all field components appear in squares,
indicating that the spin alignment measures the fluctuations of the fields
instead of their mean values.

Due to the lack of theoretical inputs for the fluctuations, the authors in
[43] choose to extract the transverse fluctuation
$F_{T}^{2}=\left\langle\mathbf{E}_{x,y}^{2}\right\rangle=\left\langle\mathbf{B}_{x,y}^{2}\right\rangle$
and longitudinal fluctuation
$F_{z}^{2}=\left\langle\mathbf{E}_{z}^{2}\right\rangle=\left\langle\mathbf{B}_{z}^{2}\right\rangle$
as functions of collision energies from experimental data for the $\phi$
meson’s spin alignment in in-plain and out-of-plain directions [8]. With these
parameters, the transverse momentum spectra of $\phi$ meson’s spin alignment
is calculated in [43]. The results are shown in the left panel of Fig. 2,
where the experiment datas are plotted by red stars with error bars, while the
model predictions are solid lines with shaded areas for uncertainties. The
rapidity dependence of $\rho_{00}$ for $\phi$ mesons has also been predicted
in [40], with the major results shown in the right panel of Fig. 2.

Figure 2: Left panel: model predictions and uncertainties [43] (solid line
with shaded areas) for the spin alignment $\rho_{00}$ of $\phi$ meson with
respect to the event plane. The red stars with error bars are experiment
results by the STAR collaboration [8]. Right panel: model predictions [40] and
uncertainties (solid line with shaded areas) for the spin alignment
$\rho_{00}$ of $\phi$ meson as functions of rapidity.

##  4 Summary

We briefly introduced the history of studies on polarizations in heavy ion
collisions. The polarization for quarks is induced by nonlocal parton
scatterings through the spin-orbit coupling, which can be analytically
described by spin kinetic equations based on the Wigner function approach. At
the hadronization stage, polarized quarks form hadrons and mesons. During this
quark recombination process, the angular momentum is converted to the
polarization of hadrons or vector mesons. For $\Lambda$ hyperons, the spin
polarization is mainly carried by the constituent $s$ quark, whose global
polarization is dominated by the vorticity field and the local polarization is
dominated by the thermal shear tensor. On the other hand, the spin alignment
of vector mesons may be mainly induced by fluctuating strong force field,
which is generated by currents of pesudo-Goldstone bosons at the hadronization
stage.

However, there are still some open questions waiting to be answered in the
future. For the global or local polarization of hyperons, we are still lack of
a comprehensive simulation based on spin kinetic theories or spin
hydrodynamics which includes non-equilibrium effects. For vector mesons, we
are wondering about possible connection between the spin alignment with the
gluon condensates or glasma fields in the initial stage of heavy-ion
collisions [47]. The contribution from strong force field [43] and some other
contributions from hydro quantities [48, 49, 50], to the spin alignment also
need detailed and comprehensive quantative studies in the future.

## References

  * [1] S.J. Barnett, _Magnetization by rotation_ , _Physical review_ 6 (1915) 239.
  * [2] A. Einstein and W. De Haas, _Experimental proof of the existence of ampère’s molecular currents_ , in _Proc. KNAW_ , vol. 181, p. 696, 1915.
  * [3] STAR collaboration, _Global $\Lambda$ hyperon polarization in nuclear collisions: evidence for the most vortical fluid_, _Nature_ 548 (2017) 62 [1701.06657].
  * [4] STAR collaboration, _Global polarization of $\Lambda$ hyperons in Au+Au collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV_, _Phys. Rev. C_ 98 (2018) 014910 [1805.04400].
  * [5] STAR collaboration, _Polarization of $\Lambda$ ($\bar{\Lambda}$) hyperons along the beam direction in Au+Au collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV_, _Phys. Rev. Lett._ 123 (2019) 132301 [1905.11917].
  * [6] ALICE collaboration, _Evidence of Spin-Orbital Angular Momentum Interactions in Relativistic Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 125 (2020) 012301 [1910.14408].
  * [7] ALICE collaboration, _First measurement of quarkonium polarization in nuclear collisions at the LHC_ , _Phys. Lett. B_ 815 (2021) 136146 [2005.11128].
  * [8] STAR collaboration, _Pattern of global spin alignment of $\phi$ and K∗0 mesons in heavy-ion collisions_, _Nature_ 614 (2023) 244 [2204.02302].
  * [9] ALICE collaboration, _Measurement of the J/ $\psi$ Polarization with Respect to the Event Plane in Pb-Pb Collisions at the LHC_, _Phys. Rev. Lett._ 131 (2023) 042303 [2204.10171].
  * [10] Q. Wang, _Global and local spin polarization in heavy ion collisions: a brief overview_ , _Nucl. Phys. A_ 967 (2017) 225 [1704.04022].
  * [11] F. Becattini and M.A. Lisa, _Polarization and Vorticity in the Quark–Gluon Plasma_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 395 [2003.03640].
  * [12] X.-G. Huang, J. Liao, Q. Wang and X.-L. Xia, _Vorticity and Spin Polarization in Heavy Ion Collisions: Transport Models_ , _Lect. Notes Phys._ 987 (2021) 281 [2010.08937].
  * [13] J.-H. Gao, Z.-T. Liang, Q. Wang and X.-N. Wang, _Global polarization effect and spin-orbit coupling in strong interaction_ , _Lect. Notes Phys._ 987 (2021) 195 [2009.04803].
  * [14] F. Becattini, _Spin and polarization: a new direction in relativistic heavy ion physics_ , _Rept. Prog. Phys._ 85 (2022) 122301 [2204.01144].
  * [15] Y. Hidaka, S. Pu, Q. Wang and D.-L. Yang, _Foundations and applications of quantum kinetic theory_ , _Prog. Part. Nucl. Phys._ 127 (2022) 103989 [2201.07644].
  * [16] Z.-T. Liang and X.-N. Wang, _Spin alignment of vector mesons in non-central A+A collisions_ , _Phys. Lett. B_ 629 (2005) 20 [nucl-th/0411101].
  * [17] Z.-T. Liang and X.-N. Wang, _Globally polarized quark-gluon plasma in non-central A+A collisions_ , _Phys. Rev. Lett._ 94 (2005) 102301 [nucl-th/0410079].
  * [18] J.-H. Gao, S.-W. Chen, W.-t. Deng, Z.-T. Liang, Q. Wang and X.-N. Wang, _Global quark polarization in non-central A+A collisions_ , _Phys. Rev. C_ 77 (2008) 044902 [0710.2943].
  * [19] J.-j. Zhang, R.-h. Fang, Q. Wang and X.-N. Wang, _A microscopic description for polarization in particle scatterings_ , _Phys. Rev. C_ 100 (2019) 064904 [1904.09152].
  * [20] H.-Z. Wu, J.-J. Zhang, L.-G. Pang and Q. Wang, _ZMCintegral: a Package for Multi-Dimensional Monte Carlo Integration on Multi-GPUs_ , _Comput. Phys. Commun._ 248 (2020) 106962 [1902.07916].
  * [21] H.T. Elze, M. Gyulassy and D. Vasak, _Transport Equations for the QCD Quark Wigner Operator_ , _Nucl. Phys. B_ 276 (1986) 706.
  * [22] D. Vasak, M. Gyulassy and H.T. Elze, _Quantum Transport Theory for Abelian Plasmas_ , _Annals Phys._ 173 (1987) 462.
  * [23] Y. Hidaka, S. Pu and D.-L. Yang, _Relativistic Chiral Kinetic Theory from Quantum Field Theories_ , _Phys. Rev. D_ 95 (2017) 091901 [1612.04630].
  * [24] N. Weickgenannt, X.-L. Sheng, E. Speranza, Q. Wang and D.H. Rischke, _Kinetic theory for massive spin-1/2 particles from the Wigner-function formalism_ , _Phys. Rev. D_ 100 (2019) 056018 [1902.06513].
  * [25] K. Hattori, Y. Hidaka and D.-L. Yang, _Axial Kinetic Theory and Spin Transport for Fermions with Arbitrary Mass_ , _Phys. Rev. D_ 100 (2019) 096011 [1903.01653].
  * [26] J.-H. Gao and Z.-T. Liang, _Relativistic Quantum Kinetic Theory for Massive Fermions and Spin Effects_ , _Phys. Rev. D_ 100 (2019) 056021 [1902.06510].
  * [27] Z. Wang, X. Guo, S. Shi and P. Zhuang, _Mass Correction to Chiral Kinetic Equations_ , _Phys. Rev. D_ 100 (2019) 014015 [1903.03461].
  * [28] N. Weickgenannt, E. Speranza, X.-l. Sheng, Q. Wang and D.H. Rischke, _Generating Spin Polarization from Vorticity through Nonlocal Collisions_ , _Phys. Rev. Lett._ 127 (2021) 052301 [2005.01506].
  * [29] D.-L. Yang, K. Hattori and Y. Hidaka, _Effective quantum kinetic theory for spin transport of fermions with collsional effects_ , _JHEP_ 07 (2020) 070 [2002.02612].
  * [30] X.-L. Sheng, N. Weickgenannt, E. Speranza, D.H. Rischke and Q. Wang, _From Kadanoff-Baym to Boltzmann equations for massive spin-1/2 fermions_ , _Phys. Rev. D_ 104 (2021) 016029 [2103.10636].
  * [31] J.-H. Gao, G.-L. Ma, S. Pu and Q. Wang, _Recent developments in chiral and spin polarization effects in heavy-ion collisions_ , _Nucl. Sci. Tech._ 31 (2020) 90 [2005.10432].
  * [32] Y. Hidaka, S. Pu and D.-L. Yang, _Nonlinear Responses of Chiral Fluids from Kinetic Theory_ , _Phys. Rev. D_ 97 (2018) 016004 [1710.00278].
  * [33] C. Yi, S. Pu and D.-L. Yang, _Reexamination of local spin polarization beyond global equilibrium in relativistic heavy ion collisions_ , _Phys. Rev. C_ 104 (2021) 064901 [2106.00238].
  * [34] F. Becattini, M. Buzzegoli and A. Palermo, _Spin-thermal shear coupling in a relativistic fluid_ , _Phys. Lett. B_ 820 (2021) 136519 [2103.10917].
  * [35] F. Becattini, M. Buzzegoli, G. Inghirami, I. Karpenko and A. Palermo, _Local Polarization and Isothermal Local Equilibrium in Relativistic Heavy Ion Collisions_ , _Phys. Rev. Lett._ 127 (2021) 272302 [2103.14621].
  * [36] S.Y.F. Liu and Y. Yin, _Spin polarization induced by the hydrodynamic gradients_ , _JHEP_ 07 (2021) 188 [2103.09200].
  * [37] B. Fu, S.Y.F. Liu, L. Pang, H. Song and Y. Yin, _Shear-Induced Spin Polarization in Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 127 (2021) 142301 [2103.10403].
  * [38] L.P. Kadanoff, _Quantum statistical mechanics_ , CRC Press (2018).
  * [39] S. Mrowczynski and U.W. Heinz, _Towards a relativistic transport theory of nuclear matter_ , _Annals Phys._ 229 (1994) 1.
  * [40] X.-L. Sheng, S. Pu and Q. Wang, _Momentum dependence of the spin alignment of the $\phi$ meson_, _Phys. Rev. C_ 108 (2023) 054902 [2308.14038].
  * [41] N. Weickgenannt, E. Speranza, X.-l. Sheng, Q. Wang and D.H. Rischke, _Derivation of the nonlocal collision term in the relativistic Boltzmann equation for massive spin-1/2 particles from quantum field theory_ , _Phys. Rev. D_ 104 (2021) 016022 [2103.04896].
  * [42] X.-L. Sheng, L. Oliva and Q. Wang, _What can we learn from the global spin alignment of $\phi$ mesons in heavy-ion collisions?_, _Phys. Rev. D_ 101 (2020) 096005 [1910.13684].
  * [43] X.-L. Sheng, L. Oliva, Z.-T. Liang, Q. Wang and X.-N. Wang, _Spin Alignment of Vector Mesons in Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 131 (2023) 042304 [2205.15689].
  * [44] X.-L. Sheng, L. Oliva, Z.-T. Liang, Q. Wang and X.-N. Wang, _Relativistic spin dynamics for vector mesons_ , 2206.05868.
  * [45] A. Manohar and H. Georgi, _Chiral Quarks and the Nonrelativistic Quark Model_ , _Nucl. Phys. B_ 234 (1984) 189.
  * [46] F. Fernandez, A. Valcarce, U. Straub and A. Faessler, _The Nucleon-nucleon interaction in terms of quark degrees of freedom_ , _J. Phys. G_ 19 (1993) 2013.
  * [47] A. Kumar, B. Müller and D.-L. Yang, _Spin alignment of vector mesons by glasma fields_ , _Phys. Rev. D_ 108 (2023) 016020 [2304.04181].
  * [48] F. Li and S.Y.F. Liu, _Tensor Polarization and Spectral Properties of Vector Meson in QCD Medium_ , 2206.11890.
  * [49] W.-B. Dong, Y.-L. Yin, X.-L. Sheng, S.-Z. Yang and Q. Wang, _Linear response theory for spin alignment of vector mesons in thermal media_ , 2311.18400.
  * [50] D. Wagner, N. Weickgenannt and E. Speranza, _Generating tensor polarization from shear stress_ , _Phys. Rev. Res._ 5 (2023) 013187 [2207.01111].