# Identifying Lightning Processes in ERA5 Soundings with Deep Learning

###### Abstract

Atmospheric environments favorable for lightning and convection are commonly
represented by proxies or parameterizations based on expert knowledge such as
CAPE, wind shears, charge separation, or combinations thereof. Recent
developments in the field of machine learning, high resolution reanalyses, and
accurate lightning observations open possibilities for identifying tailored
proxies without prior expert knowledge.

To identify vertical profiles favorable for lightning, a deep neural network
links ERA5 vertical profiles of cloud physics, mass field variables and wind
to lightning location data from the Austrian Lightning Detection & Information
System (ALDIS), which has been transformed to a binary target variable
labelling the ERA5 cells as lightning and no lightning cells. The ERA5
parameters are taken on model levels beyond the tropopause forming an input
layer of approx. 670 features. The data of 2010–2018 serve as
training/validation.

On independent test data, 2019, the deep network outperforms a reference with
features based on meteorological expertise. Shapley values highlight the
atmospheric processes learned by the network which identifies cloud ice and
snow content in the upper and mid-troposphere as relevant features. As these
patterns correspond to the separation of charge in thunderstorm cloud, the
deep learning model can serve as physically meaningful description of
lightning.

Journal of Advances in Modeling Earth Systems (JAMES)

Department of Mathematics, University of Innsbruck, Austria Department of
Atmospheric and Cryospheric Sciences, University of Innsbruck, Austria
Department of Statistics, University of Innsbruck, Austria

Tobias<EMAIL_ADDRESS>
Without prior knowledge of atmospheric environments favorable for lightning,
deep learning is capable of distinguishing between cells where lightning _or_
no lightning has been observed based on the large-scale meteorological state
of the atmosphere captured by ERA5 model level data.

The deep learning model identifies profiles associated with the separation of
charge within thunderstorm clouds from the temporally and vertically well
resolved ERA5 cloud properties of snow and ice.

In addition, the deep learning approach outperforms a reference model based on
expert variables for the domain of the European eastern Alps.

## Plain Language Summary

Lightning can have a damaging effect on infrastructure and people. As
lightning is a short and local event it is not explicitly resolved in
atmospheric models. Instead, auxiliary descriptions based on meteorological
expert knowledge, so called proxies, are used to assess lightning. We used AI
that successfully discovered on its own the ingredients that experts know to
be essential for lightning in the well-studied region of the Alps. The same
method can then be applied to regions about which much less is known.

##  1 Introduction

Lightning affects many fields of our everyday’s life. Cloud-to-ground flashes
might hit infrastructure such as wind turbines [Becerra . (2018)] and power
lines [Cummins . (1998)] and thus cause power outages. Humans might get
injured [Ritenour . (2008)] or even die [Holle (2016)] after being hit by
lightning. Ignited wildfires [Reineking . (2010)] release carbon dioxide into
the climate system, and thus limit the biosphere’s capacity to store carbon
dioxide. Lightning also affects the climate system by producing nitrogen
oxides which play a key role in ozone conversion and acid rain production
[DeCaria . (2005)]. Ozone is an important greenhouse gas and changes in
concentration can lead to warming or cooling of the atmosphere. Thus,
understanding of lightning is also an important factor in climate change
research [Finney . (2018)].

Given lightning’s impact and that an average of 46 flashes are occurring
around the globe every second [Cecil . (2014)] it is desirable to have
numerical models of the atmosphere capable to simulate lightning and its
underlying dynamic processes down to the resolved scales of the numeric model.
Beyond the resolved scales one relies on so called proxies or parameterization
to further describe lightning. The term proxy is commonly used for quantities
derived from atmospheric model output after numeric computations.
Parameterizations mean the description of lightning while numeric computations
of the atmosphere model.

Proxies are frequently applied to assess historic and future behavior of
convection and lightning. Popular proxies are cloud top height [Price Rind
(1992)], cloud ice flux [Finney . (2014)], CAPE times precipitation [Romps .
(2018)], or the lightning potential index [Brisson . (2021)]. Though, these
proxies perform reasonably good [Tippett . (2019)], there is a need for more
complex or holistic proxies, as the behavior of lightning in a changing
climate is still uncertain [Murray (2018)]. Another application that makes
clear that more research on the description of lightning is needed in the
field of operational weather forecasting. Experience shows, for instance, that
CAPE needs to be adapted to local conditions in order to perform well
[Groenemeijer . (2019)].

Parameterizations are an internal part of numeric models, as they emulate sub-
scale processes that cannot be resolved due the discretization of governing
equations. Therefore, the emulated processes give feedback to the other
processes, also on larger scales, within the atmospheric model. For instance,
tost2007 showed that modelled nitrogen oxide is sensitive to lightning
parameterizations in numerical models. Next to the classic description of
lightning using cloud top height [Price Rind (1992)], parameterizations have
been developed using polynomial regression [Allen Pickering (2002)] and
schemes based on hydrometeors in the mixed-phase region which is important for
cloud-resolving models [McCaul . (2009)]. A comparison of several
parameterizations using a superparameterized model is given by charn2021.
Recently, the ECMWF launched a product for total lightning densities expressed
as a function of hydrometeors contents, CAPE, and (convective) cloud-base
height output by the convective parameterization [Lopez (2016)].

In recent years machine learning (ML) approaches have been proposed to
describe convection and lightning. Forty preselected single level parameters
from ERA5 were processed by artificial neural networks and gradient boosting
machines for lightning in parts of Europe and Sri Lanka [Ukkonen . (2017),
Ukkonen Mäkelä (2019)]. The authors also bring up the idea to feed ERA5’s
model level parameters directly to an appropriate ML tool, i.e. neural
network. Other studies tested random forests for Hubei Province in China [Shi
. (2022)] and generalized additive models (GAM) for the european Alps [Simon .
(2022)]. All these studies confirm that the use of ML approaches for the
description of lightning is promising. In concurrent research, also
explainable artificial intelligence (XAI) techniques are used to move towards
understanding the underlying reasoning of complex AI models and shows
encouraging results in Earth System Sciences applications [Barnes . (2020),
Dutta Pal (2022), Hilburn . (2021), Mayer Barnes (2021), Stirnberg . (2021),
Toms . (2021)]. silva2022 use XGBoost classification trees to explore when the
NASA Goddard Earth Observing System model of lightning flash occurrence shows
weaknesses and apply shapley additive explanations (SHAP) to describe which
meteorological drivers are related to the model errors. They found that these
errors are strongly related to convection in the atmosphere and certain
characteristics of the land surface.

We build upon these studies and aim at finding a _holistic_ description of
lightning. Supervised deep learning harvests temporally and vertically well
resolved ERA5 soundings of atmospheric dynamics and cloud physics to explain
observations from the Austrian Lightning Information & Detection System
(ALDIS). The pattern found in ERA5 serve as proxy, but could also guide
towards a parameterization of lightning. Using ERA5 on model levels comes with
the benefit that a complete picture of the atmosphere is considered to find
patterns explaining lightning. As the approach sees the _raw model atmosphere_
, no expert parameters diagnosed from the model levels are used as inputs, the
study also answers whether deep learning can identify physically meaningful
patterns within the ERA5 sounding to describe lightning processes.

The region of interest are the eastern Alps which are characterized by complex
terrain. Atmospheric dynamics on a gamut of scales interacting with
topography, which lead to various meso-scale [Feldmann . (2021)] and local
processes [Houze (2012)] that can trigger convection and lightning. We focus
on lightning during the peak phase of the warm season (June, July, August)
which differs fundamentally in the underlying dynamic processes to lightning
during the cold season [Morgenstern . (2022)].

The paper is structured as follows. Section 2 presents both the lightning
detection data and the atmospheric reanalyses, which both enter the supervised
deep learning approach and a reference model (Sect. 3). Additionally, section
3 illustrates the methods to analyze performance and explainability of the
deep learning approach. The results of these analyses are given in section 4.
Section 5 discusses the physical patterns identified by the method,
scrutinizes further applications and research that is made possible with the
novel insights and finally concludes the study.

##  2 Data

Two data sets build the foundation for our supervised machine learning task.
First, the observational data from the lightning location system ALDIS (Sect.
2.1) is used to derive the labels distinguishing cells with and without
lightning activity. Second, pseudo soundings from ERA5 (Sect. 2.2) serve as
input for the deep learning approach.

###  2.1 Lightning Detection Data

The Austrian Lightning Detection & Information System (ALDIS) is part of the
European effort EUCLID [Schulz . (2016)]. Cloud-to-ground flashes with a
current of $>\,15\,\operatorname{kA}$ or $<\,-2\,\operatorname{kA}$ (Fig. 1)
are aggregated to the spatio-temporal grid cells of ERA5 (Sect. 2.2). Each
cell has a horizontal extent of approx.
$30\,\operatorname{km}\,\times\,30\,\operatorname{km}$ and temporally of one
hour. If at least one flash has been detected in such a grid cell, then the
cell is labelled as lightning cell. Otherwise, if not a single flash has been
detected, the cell is labelled as no lightning cell. Spatially, the grid
centers range from $8.25^{\circ}E$ to $16.75^{\circ}E$ and from
$45.25^{\circ}N$ to $49.75^{\circ}N$. Temporally, data for the meteorological
summers (June, July, August) from 2010 to 2019 are available.

Figure 1: Left: ALDIS visualization for June 20, 2019, all observed cloud to
ground flashes between 17 UTC and 18 UTC. Right: Gray crosses show the
$0.25^{\circ}$ grid of ERA5. The blue lines indicate two arbitrary samples of
vertical profiles.

###  2.2 Atmospheric Reanalysis

ECMWF’s fifth reanalyses, ERA5 [Hersbach . (2020)], has a horizontal
resolution of $0.25^{\circ}$ (approx. $30\,\operatorname{km}$) and temporally
of 1 hour. Vertically it consists of 137 hybrid model levels that align with
topography near ground and approach isobars in the upper atmosphere111See
https://confluence.ecmwf.int/display/UDOC/L137+model+level+definitions. (Fig.
1). On these model levels nine parameters (Tab. 1) are available to describe
the state of the atmosphere. Next to classical parameters, temperature,
specific humidity and three-dimensional winds, ERA5 provides a description of
liquid and solid water particles in clouds, i.e. the specific content of ice,
snow (including graupel), liquid water, and rain. For our study we use these
nine parameters (Tab. 1) on the lowest 74 model levels, spanning from level
number 64 (approx. $15\,000\,\operatorname{m}$ geopotential height) to level
number 137 ($10\,\operatorname{m}$ above ground).

Table 1: ERA5 parameters on model levels. Name | Short Name | Units | Parameter ID
---|---|---|---
Temperature | $\mathtt{t}$ | $K$ | 130
Specific humidity | $\mathtt{q}$ | $kg\,kg^{-1}$ | 133
U component of wind | $\mathtt{u}$ | $m\,s^{-1}$ | 131
V component of wind | $\mathtt{v}$ | $m\,s^{-1}$ | 132
Vertical velocity | $\mathtt{w}$ | $Pa\,s^{-1}$ | 135
Specific rain water content | $\mathtt{crwc}$ | $kg\,kg^{-1}$ | 75
Specific snow water content | $\mathtt{cswc}$ | $kg\,kg^{-1}$ | 76
Specific cloud liquid water content | $\mathtt{clwc}$ | $kg\,kg^{-1}$ | 246
Specific cloud ice water content | $\mathtt{ciwc}$ | $kg\,kg^{-1}$ | 247

###  2.3 Composition of Datasets

The two data sets are merged in order to obtain a tabular data shape. Each row
of this tabular data refers to a spatio-temporal grid cell. Thus, it can be
indexed by the longitude and latitude of its center as well as its hourly time
stamp. Each row is either labelled as lightning or no lightning cell. The nine
ERA5 parameters (Tab. 1) on their 74 model levels enter the tabular data such
that each resulting column refers to an individual parameter on an individual
level, making up a total of $9\times 74=666$ ERA5 feature columns. Further,
each row is complemented with the information of the hour of the day and day
of the season to account for diurnal and seasonal variations, respectively.
Finally, the model topography is added as another column.

##  3 Methods

To avoid incorporating expert knowledge by using specialized deep learning
architectures, a classical fully connected deep neural network (Sect. 3.1) is
used to fit a model which capable of distinguishing whether a specific spatio-
temporal grid cell corresponds to a lightning cell. To make sure that the
neural network can model lightning sufficiently well, we compare the resulting
outputs on unseen test data to those of a state-of-the-art reference model
(Sect. 3.2). Finally, we gain insights into the patterns exploited by the
trained model by applying shapley additive explanations (Sect. 3.3).

###  3.1 Deep Learning Approach

We designed a rather general fully connected neural network consisting of
eight hidden layers with $512\times 512\times 512\times 512\times 128\times
128\times 128\times 16$ nodes and use the leaky rectified linear unit (leaky
ReLU) for activation. The input dimension is predetermined by the number of
input features and thus equals $671$ (nine atmospheric variables on 74 levels,
longitude, latitude, hour of the day, day of the season, and topography). The
dimension of the output layer equals one, since we classify whether the cell
is a lightning cell and the output is activated with the sigmoid function. The
input features get standardized222The nine atmospheric variables are
standardized by considering the 74 levels altogether. prior training. To
prevent the model from overfitting, we apply dropout [Srivastava . (2014)]
with a value of $0.15$ and early stopping with a patience of ten epochs. We
use binary cross entropy as loss function with a weight of approximately $41$
for positive events (flash occurrences) to address for the highly imbalanced
data set.

###  3.2 Reference Model

For reference we used a generalized additive model [<]GAM,¿wood2017generalized
which was trained on longitude, latitude, hour of the day, day of the season,
topography and the following ten variables which were derived from ERA5
soundings on meteorological expertise [Simon . (2022)]:

  * •

$\mathtt{cape}$: Convective available potential energy,

  * •

$\mathtt{cloud\\_exists}$: Binary indicator whether cloud is present,

  * •

$\mathtt{cp}$: Convective precipitation,

  * •

$\mathtt{cswc2040}$: Mass of specific snow water content between the
$-20^{\circ}C$ and $-40^{\circ}C$ isotherms,

  * •

$\mathtt{cth}$: Cloud top height in height above ground,

  * •

$\mathtt{ishf}$: Instantaneous surface sensible heat flux,

  * •

$\mathtt{mcc}$: Medium cloud cover

  * •

$\mathtt{tcslw}$: Total column supercooled liquid ,

  * •

$\mathtt{wvc1020}$: Mass of water vapor between the $-10^{\circ}C$ and
$-20^{\circ}C$ isotherms,

  * •

$\mathtt{2t}$: Two meter temperature.

Thus, the input dimension for the reference model is only $15$. The GAM is
fitted using an algorithm tailored for gigadata [Wood . (2017)].

###  3.3 Explainability

While generalized additive models are interpretable by users [Lou . (2012)],
interpretability research of deep neural networks still suffers many gaps
[Zhang . (2021)]. We use Deep SHAP [Lundberg Lee (2017)] to gain insights into
the patterns exploited by our neural network from Section 5 which attribute to
a spatio-temporal cell being classified as a lightning cell. SHAP is a game
theoretic approach to explain the relation of input and output of any machine
learning model. It uses the concept of shapley values [Shapley (1952)] to
provide local interpretability by computing feature attributions which led to
the model’s output for a given input. Unfortunately, the computation time for
calculating the exact shapley values grows exponentially with the number of
input features. Common implementations for computing shapley values use
simplifications to ensure computational feasibility.

In this work we use Deep SHAP which is a model agnostic method that leverages
extra knowledge about the nature of deep neural networks to approximate
shapley values more efficiently but also assumes independence among the input
features. As in various applications, this property obviously is not fulfilled
in our data set but using a more accurate approximation as described in
aas2021 is not feasible with the large number of input features. Despite this
formal prerequisite, shapley values are successfully utilized in a variety of
imaging tasks (e.g. in the medical field [van der Velden . (2022)]) although
the independence assumption is also heavily violated when images are used as
input and every pixel serves as an individual input feature. Understanding the
raw ERA5 vertical profiles as a collection of 1D-images it is safe to assume
that Deep SHAP leads to sufficiently good approximations to the precise
shapley values and can therefore be applied to gain insights into the neural
network trained in this study.

##  4 Results

In this section we first investigate the performance of the deep learning
approach by comparing its output on unseen test data against observations and
the output of the reference model (Sect. 4.1). Further, the application of
SHAP allows us to gain insights about the vertical profiles exploited by the
neural network which indicate the occurrence of lightning Sect. 4.2.

###  4.1 Performance of the deep learning approach

We first trained the neural network as described in Sect. 3.1 to distinguish
whether a given spatio-temporal cell is a lightning cell or not. To map the
output of the model to a binary category, we needed to define a threshold. The
threshold was determined by applying a grid search on the F1-scores of the
validation data. Then we fitted the reference model as described in Sect. 3.2
and computed the threshold following the same procedure. The resulting
confusion matrices are displayed in Tab. 2. We observe that the neural network
slightly outperforms the reference model in every category of the confusion
matrix on previously unseen test data (year 2019). This can also be seen by
comparing the Matthew correlation coefficients (mcc) of the two models, where
$+1$ represents a perfect match between model output and observations, and $0$
no better than random guessing. The deep learning model has an mcc of
approximately $0.278$ and the reference model $0.237$.

Previous studies have shown that descriptions of lightning based on ML models
reproduce the observed diurnal cycle more realistically than simple proxies
such as CAPE [<]Fig 4. in¿ukkonen2019. To investigate the ability of our model
to reproduce the diurnal cycle of lightning we averaged the binary model
output and the binary observations of the test year 2019 for each hour of the
day and for four different small subdomains (Fig.2). The threshold to
transform the model output into a binary quantity has been chosen in such a
way that within the validation data the mean modelled occurrence of lightning
equals the mean observed occurrence.

The comparison reveals a good match of the shapes of the modelled and observed
diurnal cycles. In particular the transition from the low values in the
morning to the peak in the afternoons are well reproduced. For all four
subdomains the model slightly overestimates the observed occurrence
probabilities, which is a result of different mean occurrences of lightning in
the validation (used for finding the threshold) and test data (plotted). The
only deviation that cannot be explained by the differences between validation
and test data can be found during the late afternoon in the high Alps where
the model overestimates the observed diurnal cycle quite strongly.

Table 2: Confusion matrices of the neural network model (left) and the
reference model (right) on test year 2019.

|  | observed
---|---|---
|  | yes | no
modelled | yes | 14 370 | 61 446
no | 15 768 | 1 374 741

|  | observed
---|---|---
|  | yes | no
| yes | 12 654 | 65 176
no | 17 484 | 1 371 011

Figure 2: Diurnal cycles of the probability of a lightning occurring in a
spatio-temporal cell. Both binary model output and binary observations have
been averaged over four subdomains.

###  4.2 Identifying patterns exploited by the deep learning model

The good performance of the deep learning approach motivates us to take a
closer look at what patterns the model has learned in order to distinguish
between lightning and no lightning. Shapley values (Sect. 3.3) enable us to
get an idea about which inputs the neural network is particularly interested
in. Since we want to find patterns which are valid throughout the full region
used for training, we computed the shapley values separately333We used the
full number of samples of the corresponding spatio-temporal cell as background
data for SHAP’s DeepExplainer. for each spatio-temporal cell and aggregated
the results as visualized in Fig. 3. We observe that cloud ice (ciwc) and snow
water content (cwsc) contribute the most to the network’s output. Also note
that ciwc with its lighter-weighted ice particles is particularly interesting
at a geopotential height of approx. $9000$ to $11000\,\operatorname{m}$ and
cswc with its heavier particles at approx. $3000$ to
$10000\,\operatorname{m}$.

Figure 3: Shapley values on correctly modeled lightning events (true
positives). The two colors represent the confidence (stratified by median) of
the network in its output. The dark green color summarizes the events where
the network is very confident that a lightning event occurred. The light green
color summarizes the events where the network still modeled correctly, but
with less confidence. The solid line shows the median of all observations and
the colored areas highlight the 50% quantiles.

Having a closer look (Fig. 4) at the actual cloud ice water content and snow
water content at these heights we observe that the model is more confident
when the ciwc and cswc are significantly higher. As can be seen in Fig. 5 the
model tends to produce false positives when ciwc and cswc are high and false
negatives when ciwc and cswc are low compared to correctly classified
lightning events.

Figure 4: Cloud ice and cloud snow water content on true positives.

Figure 5: Cloud ice and cloud snow water content on less confident true
positives and false positives (top) and on less confident true positives and
false negatives (bottom). Note the different scaling compared to Fig. 4.

###  4.3 Sample case study

A sample case on unseen test data illustrates how the model from the deep
learning approach _sees_ a specific weather event. In the afternoon of June
20, 2019, a weak upper level trough embedded in southwesterly flow passed over
the Alps, whereas below crest height the flow was predominantly around the
Alps. At lower levels, the _patch_ shown in Fig. 6 was still part of the zone
with the highest equivalent potential temperature in a large warm sector,
whose delineating cold front had just reached the western European coast.

In the observations this _patch_ with lightning occurrences was propagated
with the background flow (Fig. 1). The binary model output identifies this
propagating patch of lightning, which is in the eastern part of the domain
around 18:00 UTC (Fig. 6), to a reasonable degree. The center of the patch is
well identified, whereas at its bounds misclassified cells both false positive
and false negative are present. This emphasizes the complexity of the
underlying binary classification task.

Figure 6: The map shows true positives (green diamonds), false negative (red
diamonds) and false positive (dots) for the June 20, 2019, at 18:00 UTC which
is a case of the unseen test data. The size of the green diamonds indicates
whether it is a _very_ or _less_ confident true positive. Low saturation of
the red diamonds indicates that the output of the network was close to
labelling the cell as _lightning_ cell. Figure 7: Cloud top height of very
confident true positives, less confident true positives, false negatives and
false positives. Clearly, the model is more confidence when the value of cloud
top height becomes larger.

##  5 Discussion and Conclusions

In this study we trained a neural network on the vertical columns of raw ERA5
data without inducing any further expert knowledge about atmospheric processes
to classify whether there was a lightning event or not. We then used shapley
values to explain which variables and vertical levels attribute the most to
correct classifications of lightning cells. As shown in Sect. 4.2, the
specific snow water and ice water content by far attracts the most attention
and the interest peaks at a geopotential height of approx.
$7000\,\operatorname{m}$ and at $10000\,\operatorname{m}$ respectively. The
neural network discovered by itself the essential ingredient for lightning,
namely charge separation. It occurs when ice cyrstals (ciwc) and larger frozen
particles (graupel, cswc) are present in the convective updraft. Once the
graupel is sufficiently heavy, its velocity is smaller than the velocity of
the rising ice crystals, and the collisions between ice crystals and graupel
result in oppositely charged particles [Reynolds . (1957), Saunders . (2006)].
¡Fig. 1¿lopez2016 shows the typical distribution of charges in a mature
thunderstorm cloud. However, it is noteworthy that the model seems to be
particularly interested in the cloud ice water content at a height of $9000$
to $11000\,\operatorname{m}$ while recent literature usually looks at the
cloud ice water content at $440\,\operatorname{hPa}$ (typically about
$6000\,\operatorname{m}$) [Finney . (2014), Finney . (2018), Silva . (2022)].
While important proxies such as cloud top height were not explicitly given to
the neural network, Fig. 7 illustrates that it was able to extract that
information from the full vertical profiles and utilize this feature. Focusing
on the region close to the tropopause between 9000 and 11000 m means that it
is crucial to vent ice particles all the way up to the tropopause and form
anvils as is typical of the thunderstorm clouds.

Our results suggest promising future applications. Being able to train a
neural network directly on atmospheric soundings with good ability to
distinguish lighting and non lightning, and then opening the black box may
enable researchers to gain a better understanding of atmospheric processes in
regions like e.g. equatorial Africa where ample studies are scarce
[Chakraborty . (2022)]. The launch of the first MGT-I satellite is expected at
the end of 2022 and will provide a lightning imager [Holmlund . (2021)] which
appears to be a promising source for the target variable. Furthermore, many
existing models come with two very different parameterizations for ocean and
land [Finney . (2014)] and this inevitably leads to discontinuities in coastal
areas. Also the reasons for the much lower lightning frequency over ocean are
not as well understood yet. Explainable AI might be a valuable building block
in moving towards a more holistic understanding of the underlying atmospheric
processes.

Future work might improve the results presented in this study. Convection and
cloud processes are not purely vertical processes and thus ML parameterization
greatly benefits from using multiple neighboring vertical atmospheric columns
instead of a single column. wang2022 work with $192\,\operatorname{km}\times
192\,\operatorname{km}$ grid cells to model, among others, subgrid zonal and
meridional momentum flux due to vertical advection and suggest that a $3\times
3$ subgrid could further improve the performance of the deep learning
approach. We used a simple fully connected neural network and therefore lost
information about the connectivity of the values along the levels of the
vertical profiles. Using convolutional layers to process the profiles would,
most likely, further improve the results. However, the goal of this work was
to use a very simple machine learning approach to detect lightning cells and
then to break up the model and understand which atmospheric conditions the
model has found to be typical for lightning. The input data was preprocessed
with only very little meteorological expertise to ensure that the methodology
is easily transferable to other regions of the earth where the understanding
of lightning related atmospheric processes is still scarce.

##  6 Open Research

The Python and R code used to produce the results and plots in this study are
available at https://github.com/noxthot/xai_lightningprocesses. The source
code relies on two data sources. ERA5 [Hersbach . (2020)] data are available
via the Climate Data Store (https://cds.climate.copernicus.eu/). Scripts for
sending the retrievals are included in the data-preprocessing directory of the
GitHub repo. The ALDIS data [Schulz . (2016)], which are the second important
source of data, cannot be made available to the public. However, ALDIS data
are available on request from ALDIS<EMAIL_ADDRESS>– fees may be charged.

###### Acknowledgements.

This work was funded by the Austrian Science Fund (FWF, grant no. P 31836) and
the Austrian Research Promotion Agency (FFG, grant no. 872656). We are
grateful for data support by Gerhard Diendorfer and Wolfgang Schulz from OVE-
ALDIS. We thank Deborah Morgenstern and Johannes Horak for their script to
compute geopotential height on ERA5 model levels. Also, we thank Johanna
Rissbacher for contributing Fig. 6 and the corresponding code.

## References

  * Aas . (2021) aas2021Aas, K., Jullum, M. Løland, A.  2021\.  Explaining Individual Predictions when Features are Dependent: More Accurate Approximations to Shapley Values Explaining individual predictions when features are dependent: More accurate approximations to Shapley values.  Artificial Intelligence298103502.  10.1016/j.artint.2021.103502
  * Allen Pickering (2002) allen2002Allen, DJ. Pickering, KE.  2002\.  Evaluation of Lightning Flash Rate Parameterizations for Use in a Global Chemical Transport Model Evaluation of lightning flash rate parameterizations for use in a global chemical transport model.  Journal of Geophysical Research: Atmospheres107D23ACH 15-1-ACH 15-21.  10.1029/2002JD002066
  * Barnes . (2020) barnes2020Barnes, EA., Toms, B., Hurrell, JW., Ebert-Uphoff, I., Anderson, C. Anderson, D.  2020\.  Indicator Patterns of Forced Change Learned by an Artificial Neural Network Indicator patterns of forced change learned by an artificial neural network.  Journal of Advances in Modeling Earth Systems129e2020MS002195.  10.1029/2020MS002195
  * Becerra . (2018) becerra2018Becerra, M., Long, M., Schulz, W. Thottappillil, R.  2018\.  On the Estimation of the Lightning Incidence to Offshore Wind Farms On the estimation of the lightning incidence to offshore wind farms.  Electric Power Systems Research157211–226.  10.1016/j.epsr.2017.12.008
  * Brisson . (2021) brisson2021Brisson, E., Blahak, U., Lucas-Picher, P., Purr, C. Ahrens, B.  2021\.  Contrasting Lightning Projection Using the Lightning Potential Index Adapted in a Convection-Permitting Regional Climate Model Contrasting lightning projection using the lightning potential index adapted in a convection-permitting regional climate model.  Climate Dynamics5772037–2051.  10.1007/s00382-021-05791-z
  * Cecil . (2014) cecil2014Cecil, DJ., Buechler, DE. Blakeslee, RJ.  2014\.  Gridded Lightning Climatology from TRMM-LIS and OTD: Dataset Description Gridded lightning climatology from TRMM-LIS and OTD: Dataset description.  Atmospheric Research135404–414.  10.1016/j.atmosres.2012.06.028
  * Chakraborty . (2022) chakraborty2022Chakraborty, R., Menghal, P., Harshitha, M. Sodunke, M.  2022\.  Climatology of Lightning Activities Across the Equatorial African Region Climatology of lightning activities across the Equatorial African region.  2022 3rd URSI Atlantic and Asia Pacific Radio Science Meeting (AT-AP-RASC) 2022 3rd ursi atlantic and asia pacific radio science meeting (at-ap-rasc) ( 1–4).  10.23919/AT-AP-RASC54737.2022.9814276
  * Charn Parishani (2021) charn2021Charn, AB. Parishani, H.  2021\.  Predictive Proxies of Present and Future Lightning in a Superparameterized Model Predictive proxies of present and future lightning in a superparameterized model.  Journal of Geophysical Research: Atmospheres12617.  10.1029/2021JD035461
  * Cummins . (1998) cummins1998Cummins, K., Krider, E. Malone, M.  1998\.  The US National Lightning Detection Network and Applications of Cloud-to-Ground Lightning Data by Electric Power Utilities The US National Lightning Detection Network and applications of cloud-to-ground lightning data by electric power utilities.  IEEE Transactions on Electromagnetic Compatibility404465–480.  10.1109/15.736207
  * DeCaria . (2005) decaria2005DeCaria, AJ., Pickering, KE., Stenchikov, GL. Ott, LE.  2005\.  Lightning-Generated $\mathrm{NO}_{x}$ and its Impact on Tropospheric Ozone Production: A Three-Dimensional Modeling Study of a Stratosphere-Troposphere Experiment: Radiation, Aerosols and Ozone (STERAO-A) Thunderstorm Lightning-generated $\mathrm{NO}_{x}$ and its impact on tropospheric ozone production: A three-dimensional modeling study of a Stratosphere-Troposphere Experiment: Radiation, Aerosols and Ozone (STERAO-A) thunderstorm.  Journal of Geophysical Research: Atmospheres110D14.  10.1029/2004JD005556
  * Dutta Pal (2022) dutta2022Dutta, D. Pal, SK.  2022\.  Interpretation of Black Box for Short-Term Predictions of Pre-Monsoon Cumulonimbus Cloud Events over Kolkata Interpretation of black box for short-term predictions of pre-monsoon cumulonimbus cloud events over Kolkata.  Journal of Data, Information and Management42167–183.  10.1007/s42488-022-00071-9
  * Feldmann . (2021) feldmann2021Feldmann, M., Germann, U., Gabella, M. Berne, A.  2021\.  A Characterisation of Alpine Mesocyclone Occurrence A characterisation of Alpine mesocyclone occurrence.  Weather and Climate Dynamics241225–1244.  10.5194/wcd-2-1225-2021
  * Finney . (2014) finney2014Finney, DL., Doherty, RM., Wild, O., Huntrieser, H., Pumphrey, HC. Blyth, AM.  2014\.  Using Cloud Ice Flux to Parametrise Large-Scale Lightning Using cloud ice flux to parametrise large-scale lightning.  Atmospheric Chemistry and Physics142312665–12682.  10.5194/acp-14-12665-2014
  * Finney . (2018) finney2018Finney, DL., Doherty, RM., Wild, O., Stevenson, DS., MacKenzie, IA. Blyth, AM.  2018\.  A Projected Decrease in Lightning under Climate Change A projected decrease in lightning under climate change.  Nature Climate Change83210–213.  10.1038/s41558-018-0072-6
  * Groenemeijer . (2019) groenemeijer2019ecmwfmemoGroenemeijer, P., Púcik, T., Tsonevsky, I. Bechtold, P.  201911\.  An Overview of Convective Available Potential Energy and Convective Inhibition provided by NWP models for operational forecasting An overview of convective available potential energy and convective inhibition provided by NWP models for operational forecasting.  852\.  10.21957/q392hofrl
  * Hersbach . (2020) hersbach2020Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J.Thépaut, JN.  2020\.  The ERA5 Global Reanalysis The ERA5 global reanalysis.  Quarterly Journal of the Royal Meteorological Society1467301999–2049.  10.1002/qj.3803
  * Hilburn . (2021) hilburn2021Hilburn, KA., Ebert-Uphoff, I. Miller, SD.  2021\.  Development and Interpretation of a Neural-Network-Based Synthetic Radar Reflectivity Estimator Using GOES-R Satellite Observations Development and interpretation of a neural-network-based synthetic radar reflectivity estimator using GOES-R satellite observations.  Journal of Applied Meteorology and Climatology6013–21.  10.1175/JAMC-D-20-0084.1
  * Holle (2016) holle2016Holle, RL.  2016\.  A Summary of Recent National-Scale Lightning Fatality Studies A summary of recent national-scale lightning fatality studies.  Weather, Climate, and Society8135–42.  10.1175/WCAS-D-15-0032.1
  * Holmlund . (2021) holmlund2021Holmlund, K., Grandell, J., Schmetz, J., Stuhlmann, R., Bojkov, B., Munro, R.Blythe, P.  2021\.  Meteosat Third Generation (MTG): Continuation and Innovation of Observations from Geostationary Orbit Meteosat Third Generation (MTG): Continuation and innovation of observations from geostationary orbit.  Bulletin of the American Meteorological Society1025990–1015.  10.1175/BAMS-D-19-0304.1
  * Houze (2012) houze2012Houze, RA.  2012\.  Orographic Effects on Precipitating Clouds Orographic effects on precipitating clouds.  Reviews of Geophysics5011–47.  10.1029/2011RG000365
  * Lopez (2016) lopez2016Lopez, P.  2016\.  A Lightning Parameterization for the ECMWF Integrated Forecasting System A lightning parameterization for the ECMWF integrated forecasting system.  Monthly Weather Review14493057–3075.  10.1175/MWR-D-16-0026.1
  * Lou . (2012) lou2017Lou, Y., Caruana, R. Gehrke, J.  2012\.  Intelligible Models for Classification and Regression Intelligible models for classification and regression.  Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining Proceedings of the 18th acm sigkdd international conference on knowledge discovery and data mining ( 150–158).  New York, NY, USAAssociation for Computing Machinery.  10.1145/2339530.2339556
  * Lundberg Lee (2017) lundberg2017Lundberg, SM. Lee, SI.  2017\.  A Unified Approach to Interpreting Model Predictions A unified approach to interpreting model predictions.  I. Guyon . (), Advances in Neural Information Processing Systems Advances in neural information processing systems ( 30).  Curran Associates, Inc.  https://proceedings.neurips.cc/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf
  * Mayer Barnes (2021) mayer2021Mayer, KJ. Barnes, EA.  2021\.  Subseasonal Forecasts of Opportunity Identified by an Explainable Neural Network Subseasonal forecasts of opportunity identified by an explainable neural network.  Geophysical Research Letters4810e2020GL092092.  10.1029/2020GL092092
  * McCaul . (2009) mccaul2009McCaul, EW., Goodman, SJ., LaCasse, KM. Cecil, DJ.  2009\.  Forecasting Lightning Threat Using Cloud-Resolving Model Simulations Forecasting lightning threat using cloud-resolving model simulations.  Weather and Forecasting243709–729.  10.1175/2008WAF2222152.1
  * Morgenstern . (2022) morgenstern2022Morgenstern, D., Stucke, I., Simon, T., Mayr, GJ. Zeileis, A.  2022\.  Differentiating Lightning in Winter and Summer with Characteristics of the Wind Field and Mass Field Differentiating lightning in winter and summer with characteristics of the wind field and mass field.  Weather and Climate Dynamics31361–375.  10.5194/wcd-3-361-2022
  * Murray (2018) murray2018Murray, LT.  2018\.  An Uncertain Future for Lightning An uncertain future for lightning.  Nature Climate Change83191–192.  10.1038/s41558-018-0094-0
  * Price Rind (1992) price1992Price, C. Rind, D.  1992\.  A Simple Lightning Parameterization for Calculating Global Lightning Distributions A simple lightning parameterization for calculating global lightning distributions.  Journal of Geophysical Research: Atmospheres97D99919–9933.  10.1029/92JD00719
  * Reineking . (2010) reineking2010Reineking, B., Weibel, P., Conedera, M. Bugmann, H.  2010\.  Environmental Determinants of Lightning- v. Human-Induced Forest Fire Ignitions Differ in a Temperate Mountain Region of Switzerland Environmental determinants of lightning- v. human-induced forest fire ignitions differ in a temperate mountain region of Switzerland.  International Journal of Wildland Fire195541–557.  10.1071/WF08206
  * Reynolds . (1957) reynolds1957Reynolds, S., Brook, M. Gourley, MF.  1957\.  Thunderstorm Charge Separation Thunderstorm charge separation.  Journal of Atmospheric Sciences145426–436.  10.1175/1520-0469(1957)014¡0426:TCS¿2.0.CO;2
  * Ritenour . (2008) ritenour2008Ritenour, AE., Morton, MJ., McManus, JG., Barillo, DJ. Cancio, LC.  2008\.  Lightning Injury: A Review Lightning injury: A review.  Burns345585–594.  10.1016/j.burns.2007.11.006
  * Romps . (2018) romps2018Romps, DM., Charn, AB., Holzworth, RH., Lawrence, WE., Molinari, J. Vollaro, D.  2018\.  CAPE Times P Explains Lightning Over Land But Not the Land-Ocean Contrast CAPE times P explains lightning over land but not the land-ocean contrast.  Geophysical Research Letters452212,623–12,630.  10.1029/2018GL080267
  * Saunders . (2006) saunders2006Saunders, CPR., Bax-norman, H., Emersic, C., Avila, EE. Castellano, NE.  2006\.  Laboratory Studies of the Effect of Cloud Conditions on Graupel/Crystal Charge Transfer in Thunderstorm Electrification Laboratory studies of the effect of cloud conditions on graupel/crystal charge transfer in thunderstorm electrification.  Quarterly Journal of the Royal Meteorological Society1326212653-2673.  10.1256/qj.05.218
  * Schulz . (2016) schulz2016Schulz, W., Diendorfer, G., Pedeboy, S. Poelman, DR.  2016\.  The European Lightning Location System EUCLID Part 1: Performance Analysis and Validation The European lightning location system EUCLID part 1: Performance analysis and validation.  Natural Hazards and Earth System Sciences162595–605.  10.5194/nhess-16-595-2016
  * Shapley (1952) shapley1952Shapley, LS.  1952\.  A Value for N-Person Games A value for n-person games.  Santa Monica, CARAND Corporation.  10.7249/P0295
  * Shi . (2022) shi2022Shi, M., Zhang, W., Fan, P., Chen, Q., Liu, Z., Li, Q. Liu, X.  2022\.  Modelling Deep Convective Activity Using Lightning Clusters and Machine Learning Modelling deep convective activity using lightning clusters and machine learning.  International Journal of Climatology422952-973.  10.1002/joc.7282
  * Silva . (2022) silva2022Silva, SJ., Keller, CA. Hardin, J.  2022\.  Using an Explainable Machine Learning Approach to Characterize Earth System Model Errors: Application of SHAP Analysis to Modeling Lightning Flash Occurrence Using an explainable machine learning approach to characterize earth system model errors: Application of SHAP analysis to modeling lightning flash occurrence.  Journal of Advances in Modeling Earth Systems144e2021MS002881.  10.1029/2021MS002881
  * Simon . (2022) simon2022Simon, T., Mayr, GJ., Morgenstern, D., Umlauf, N. Zeileis, A.  2022\.  Amplification of Annual and Diurnal Cycles of Alpine Lightning Amplification of annual and diurnal cycles of alpine lightning [Preprint].  10.21203/rs.3.rs-965951/v3
  * Srivastava . (2014) srivastava2014Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. Salakhutdinov, R.  2014\.  Dropout: A Simple Way to Prevent Neural Networks from Overfitting Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research15561929–1958.  http://jmlr.org/papers/v15/srivastava14a.html
  * Stirnberg . (2021) stirnberg2021Stirnberg, R., Cermak, J., Kotthaus, S., Haeffelin, M., Andersen, H., Fuchs, J.Favez, O.  2021\.  Meteorology-Driven Variability of Air Pollution (PM 1) Revealed with Explainable Machine Learning Meteorology-driven variability of air pollution (PM 1) revealed with explainable machine learning.  Atmospheric Chemistry and Physics2153919–3948.  10.5194/acp-21-3919-2021
  * Tippett . (2019) tippett2019Tippett, MK., Lepore, C., Koshak, WJ., Chronis, T. Vant-Hull, B.  2019\.  Performance of a Simple Reanalysis Proxy for U.S. Cloud-to-Ground Lightning Performance of a simple reanalysis proxy for U.S. cloud-to-ground lightning.  International Journal of Climatology39103932–3946.  10.1002/joc.6049
  * Toms . (2021) toms2021Toms, BA., Barnes, EA. Hurrell, JW.  2021\.  Assessing Decadal Predictability in an Earth-System Model Using Explainable Neural Networks Assessing decadal predictability in an earth-system model using explainable neural networks.  Geophysical Research Letters4812e2021GL093842.  10.1029/2021GL093842
  * Tost . (2007) tost2007Tost, H., Jöckel, P. Lelieveld, J.  2007\.  Lightning and Convection Parameterisations – Uncertainties in Global Modelling Lightning and convection parameterisations – uncertainties in global modelling.  Atmospheric Chemistry and Physics7174553–4568.  10.5194/acp-7-4553-2007
  * Ukkonen Mäkelä (2019) ukkonen2019Ukkonen, P. Mäkelä, A.  2019\.  Evaluation of Machine Learning Classifiers for Predicting Deep Convection Evaluation of machine learning classifiers for predicting deep convection.  J. Adv. Model. Earth Sy.1161784–1802.  10.1029/2018MS001561
  * Ukkonen . (2017) ukkonen2017Ukkonen, P., Manzato, A. Mäkelä, A.  2017\.  Evaluation of Thunderstorm Predictors for Finland Using Reanalyses and Neural Networks Evaluation of thunderstorm predictors for Finland using reanalyses and neural networks.  ournal of Applied Meteorology and Climatology5682335–2352.  10.1175/JAMC-D-16-0361.1
  * van der Velden . (2022) vandervelden2022van der Velden, BH., Kuijf, HJ., Gilhuijs, KG. Viergever, MA.  2022\.  Explainable artificial intelligence (XAI) in deep learning-based medical image analysis Explainable artificial intelligence (xai) in deep learning-based medical image analysis.  Medical Image Analysis79102470.  10.1016/j.media.2022.102470
  * Wang . (2022) wang2022Wang, P., Yuval, J. O’Gorman, PA.  2022\.  Non-local parameterization of atmospheric subgrid processes with neural networks Non-local parameterization of atmospheric subgrid processes with neural networks.  Journal of Advances in Modeling Earth Systemse2022MS002984.  10.1029/2022MS002984
  * Wood (2017) wood2017generalizedWood, SN.  2017\.  Generalized Additive Models: An Introduction with R Generalized additive models: An introduction with R (2nd ).  Boca RatonChapman & Hall/CRC.  10.1201/9781420010404
  * Wood . (2017) wood2017gigadataWood, SN., Li, Z., Shaddick, G. Augustin, NH.  2017\.  Generalized Additive Models for Gigadata: Modeling the U.K. Black Smoke Network Daily Data Generalized additive models for gigadata: Modeling the U.K. black smoke network daily data.  Journal of the American Statistical Association1125191199–1210.  10.1080/01621459.2016.1195744
  * Zhang . (2021) zhang2021Zhang, Y., Tiňo, P., Leonardis, A. Tang, K.  2021\.  A Survey on Neural Network Interpretability A survey on neural network interpretability.  IEEE Transactions on Emerging Topics in Computational Intelligence55726-742.  10.1109/TETCI.2021.3100641


2022 121 07 June

Flow angle and flow magnitude fluctuations in Pb–Pb collisions at the LHC
ALICE CollaborationALICE Collaboration

This Letter reports on the first measurements of transverse momentum dependent
flow angle $\Psi_{n}$ and flow magnitude $v_{n}$ fluctuations, determined
using new four-particle correlators. The measurements are performed for
various centralities in Pb–Pb collisions at a centre-of-mass energy of
$\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV with ALICE at the LHC. Both flow angle and
flow magnitude fluctuations are observed in the presented centrality ranges
and are strongest in the most central collisions and for a transverse momentum
$p_{\rm T}>2$ GeV/c. Comparison with theoretical models, including iEBE-
VISHNU, MUSIC, and AMPT, show that the measurements bring novel insights into
the fluctuating initial conditions that are not well known. In addition, these
new results exhibit unique sensitivities to the specific shear viscosity,
$\eta/s$, of the quark–gluon plasma (QGP) and to the initial state of the
heavy-ion collisions. As such fluctuations are getting stronger with
increasing $p_{\rm T}$, a re-examination of existing models is needed to have
a more precise and unbiased extraction of properties of the QGP.

In ultrarelativistic collisions of heavy ions, such as those at the
Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC), a
deconfined state of strongly interacting matter, commonly referred to as
quark–gluon plasma (QGP), is predicted to be created under extreme conditions
of temperature and energy densities [1, 2]. Many experimental results indicate
that a strongly-coupled QGP is formed in heavy-ion collisions [3, 4, 5, 6, 7].
Initial anisotropies of the geometric overlap of the colliding nuclei and
spatial inhomogeneities in the energy density drive the collective expansion
of the QGP and are transformed through the evolution of the QGP into a
momentum anisotropy in the final state [8, 9, 10]. This momentum anisotropy is
characterised by a Fourier expansion of the distribution of the azimuthal
angle, $\varphi$, of emitted particles [11]

$\displaystyle\frac{\mathrm{d}^{2}N}{\mathrm{d}p_{\rm
T}\mathrm{d}\varphi}=\frac{\mathrm{d}N}{2\pi\mathrm{d}p_{\rm
T}}\left(1+2\sum_{n=1}^{\infty}v_{n}(p_{\rm T})\cos[n(\varphi-\Psi_{n}(p_{\rm
T}))]\right),$ (1)

where $v_{n}(p_{\rm T})$ and $\Psi_{n}(p_{\rm T})$ correspond to the magnitude
and angle, respectively, of the $n^{\mathrm{th}}$-order harmonic flow vector
$\vec{V}_{n}(p_{\rm T})=v_{n}(p_{\rm T})e^{in\Psi_{n}(p_{\rm T})}$. Here the
transverse momentum, $p_{\rm T}$, dependence of both the flow magnitude and
flow angle has been made explicit. The flow vector quantifies the orientation
and magnitude of the anisotropic flow, and the flow angle $\Psi_{n}$ is the
angle of the symmetry plane of the $n^{\mathrm{th}}$-order flow vector.

Typically, the largest flow coefficient is the elliptic flow $v_{2}$, since it
is largely determined by the geometrical overlap of nuclei colliding at a
finite impact parameter, the distance between the nuclei centres in the
transverse plane. However, fluctuations in the position of the colliding
nuclei and in the position of nucleons within the nuclei can induce more
complex geometrical shapes, which will give rise to non-zero flow coefficients
with $n>2$ [12, 13, 14, 15], such as the triangular flow $v_{3}$. Both
elliptic and triangular flow coefficients have been measured extensively at
RHIC [16, 17, 18, 19] and at the LHC [20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30]. The comparison of measurements of $v_{n}$ at low $p_{\rm T}$ with
hydrodynamical calculations has been used to constrain the initial conditions
of the heavy-ion collisions and the transport properties of the QGP, such as
the specific shear viscosity $\eta/s$. These comparisons indicate that the
system behaves as a strongly-coupled low-viscosity fluid [10, 31, 32, 33].
Studies probing the collective expansion of the system using hard probes, such
as measurements of the anisotropic flow of charged and identified hadrons at
high $p_{\rm T}$ [34, 35, 36, 37, 38, 39], or of hadrons containing heavy
quarks [40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52], rely on the
comparison of theoretical model calculations to experimental data to extract
the QGP properties, such as the path-length dependence of the heavy-quark
energy loss [53, 54]. In this Letter, the observation of separate flow angle
and flow magnitude fluctuations in heavy-ion data is presented. It is argued
that the notion of a common symmetry plane, widely used in various flow
measurements, does not hold in measurements that correlate particles at high
$p_{\rm T}$ (> 2 GeV) with the common symmetry plane determined mainly by low
$p_{\rm T}$ particles.

Fluctuations of both the flow angle $\Psi_{n}$ and the flow magnitude $v_{n}$
have been shown to be present in hydrodynamical models [55, 56]. The flow
angle fluctuations are the fluctuations of the symmetry plane determined by a
subset of the particles, in this case, those at a specific $p_{\rm T}$, with
respect to the symmetry plane determined by the total set of particles. In
heavy-ion collisions, the produced particle yield is typically dominated by
low $p_{\rm T}$ particles. The flow magnitude fluctuations can be understood
as the fluctuations of the $v_{n}$ coefficients at different transverse
momentum ranges. Measurements at the LHC have reported the existence of
$p_{\rm T}$-dependent flow vector fluctuations [57, 58, 59, 60]. However,
those measurements rely on observables constructed from two-particle
correlations, which intrinsically contain contributions from both the flow
angle and magnitude fluctuations with no way to separate the two
experimentally. A Principal Component Analysis (PCA) has been used for the
extraction of leading and subleading modes of flow fluctuations in heavy-ion
collisions [61, 62, 63, 64]. However, these analyses were also based on
information retrieved from two-particle correlations and cannot be used to
separate the individual components of the flow vector fluctuations. In this
Letter, we propose new four-particle correlation functions to separate the
$p_{\rm T}$-dependent fluctuations of the flow angle and flow magnitude. The
flow angle fluctuations can be quantified by

$\displaystyle A^{\mathrm{f}}_{n}=\frac{\langle\langle\cos
n[\varphi_{1}^{\mathrm{a}}+\varphi_{2}^{\mathrm{a}}-\varphi_{3}-\varphi_{4}]\rangle\rangle}{\langle\langle\cos
n[\varphi_{1}^{\mathrm{a}}+\varphi_{2}-\varphi_{3}^{\mathrm{a}}-\varphi_{4}]\rangle\rangle}$
$\displaystyle=\frac{\langle
v_{n}^{2}(p_{\mathrm{T}}^{\mathrm{a}})~{}v_{n}^{2}\cos
2n[\Psi_{n}(p_{\mathrm{T}}^{\mathrm{a}})-\Psi_{n}]\rangle}{\langle
v_{n}^{2}(p_{\mathrm{T}}^{\mathrm{a}})v_{n}^{2}\rangle}$
$\displaystyle\simeq\langle\cos 2n\left[\Psi_{n}(p_{\rm T}^{\rm
a})-\Psi_{n}\right]\rangle_{\mathrm{w}},$ (2)

where the $a$ superscript refers to an associate particle selected from a
narrow transverse momentum range $p_{\rm T}^{\rm a}$ and the $w$ means that
$A^{\rm f}_{n}$ is averaged over the event ensemble with each event having a
weight of the fourth power of $v_{n}$ [65]. The double brackets refer to an
average over all particles and all events, and the single brackets refer to an
average over all events. A value of $A^{\rm f}_{n}<1$ indicates the presence
of $p_{\rm T}$-dependent flow angle fluctuations. A larger deviation from
unity suggests that the symmetry plane at a specific $p_{\rm T}$ is more
decorrelated with respect to the common symmetry plane. On the other hand, the
flow magnitude fluctuations can be studied with the double ratio

$\displaystyle M^{\rm f}_{n}$ $\displaystyle=\frac{\langle\langle\cos
n[\varphi_{1}^{\mathrm{a}}+\varphi_{2}-\varphi_{3}^{\mathrm{a}}-\varphi_{4}]\rangle\rangle/\left(\langle\langle\cos
n[\varphi_{1}^{\mathrm{a}}-\varphi_{3}^{\mathrm{a}}]\rangle\rangle\langle\langle\cos
n[\varphi_{2}-\varphi_{4}]\rangle\rangle\right)}{\langle\langle\cos
n[\varphi_{1}+\varphi_{2}-\varphi_{3}-\varphi_{4}]\rangle\rangle/\langle\langle\cos
n[\varphi_{1}-\varphi_{2}]\rangle\rangle^{2}}=\frac{\langle v_{n}^{2}(p_{\rm
T}^{\rm a})v_{n}^{2}\rangle/\langle v_{n}^{2}(p_{\rm T}^{\rm a})\rangle\langle
v_{n}^{2}\rangle}{\langle v_{n}^{4}\rangle/\langle v_{n}^{2}\rangle^{2}}.$ (3)

A deviation of $M^{\rm f}_{n}$ from unity indicates the presence of $p_{\rm
T}$-dependent flow magnitude fluctuations. The magnitude of the deviation will
show how strongly the flow magnitude is decorrelated at a specific $p_{\rm T}$
range with respect to the integrated flow. These new correlators probe a
higher moment of the distribution of flow fluctuations compared to correlators
traditionally used in studies of flow and flow fluctuations using two-particle
techniques [57, 58, 59]. The lower-order moments of the flow angle and flow
magnitude fluctuations cannot be measured separately in experiments but can be
approximated by constructing the lower and upper limits of the first moment of
flow angle and magnitude fluctuations, respectively. The lower limit of the
first-moment flow angle fluctuations $\langle\cos n\left[\Psi_{n}(p_{\rm
T}^{\rm a})-\Psi_{n}\right]\rangle$ can be calculated as

$\displaystyle\sqrt{\frac{A^{\rm f}_{n}+1}{2}}\geq\langle\cos
n[\Psi_{n}(p_{\rm T}^{\rm a})-\Psi_{n}]\rangle.$ (4)

The overall flow vector fluctuations can be probed by the ratio of the $p_{\rm
T}$-differential flow coefficient $v_{n}\\{2\\}=\langle v_{n}(p_{\rm T}^{\rm
a})v_{n}\cos n[\Psi_{n}(p_{\rm T}^{\rm a})-\Psi_{n}]\rangle/\sqrt{\langle
v_{n}^{2}\rangle}$ and the $p_{\rm T}$-integrated flow coefficient in a narrow
$p_{\rm T}$ interval $v_{n}[2]=\sqrt{\langle v_{n}^{2}(p_{\rm T}^{\rm
a})\rangle}$ [55]. Taking the ratio of Eq. (4) and the flow vector
fluctuations

$\displaystyle v_{n}\\{2\\}/v_{n}[2]=\frac{\langle v_{n}(p_{\rm T}^{\rm
a})~{}v_{n}\cos n\left[\Psi_{n}(p_{\rm T}^{\rm
a})-\Psi_{n}\right]\rangle}{\sqrt{\langle v_{n}^{2}(p_{\rm T}^{\rm
a})\rangle)}\sqrt{\langle v_{n}^{2}\rangle}},$ (5)

enables the calculation of the upper limit of the first-order flow magnitude
fluctuations

$\displaystyle\frac{v_{n}\\{2\\}/v_{n}[2]}{\sqrt{(A^{\rm
f}_{n}+1)/2}}\leq\frac{\langle v_{n}(p_{\rm T}^{\rm
a})~{}v_{n}\rangle}{\sqrt{\langle v_{n}^{2}(p_{\rm T}^{\rm
a})\rangle}\sqrt{\langle v_{n}^{2}\rangle}}.$ (6)

The limits on the first-moment flow angle and magnitude fluctuations connect
the study of the separated fluctuations with the previous studies of flow
vector fluctuations using only two-particle correlations [57, 58, 59]. All the
above correlators are calculated with the generic framework [66, 67], which
corrects the non-uniformities in the acceptance of the detector. As the
correlations in Eqs. (2) and (3) are themselves correlated, the statistical
uncertainty is estimated by the bootstrap method of random sampling with
replacement.

The correlators $A^{\rm f}_{2}$ and $M^{\rm f}_{2}$ are measured in
approximately 54 million Pb–Pb collisions collected in 2015 at a centre-of-
mass energy per nucleon pair of $\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV.
Experimentally, events are selected based on a minimum bias trigger achieved
by requiring a coincidence of signals in the two V0 scintillator arrays, V0A
with a pseudorapidity range $2.8<\eta<5.1$, and V0C with a pesudorapidity
range $-3.7<\eta<-1.7$. Additionally, a reconstructed primary vertex within
$\pm 10$ cm of the nominal interaction point along the beam axis is required.
Rejection of pileup events, which are events close in time to the triggered
event due to the high luminosity, is based on the correlation between tracks
measured with different detectors. A variation of the criteria for pileup
rejection is considered for the systematic uncertainties [68]. The centrality
of the events is measured using a fit of the V0A and V0C signal amplitude sum,
with a two-component model for particle production and a Glauber Monte Carlo
model of the collision geometry [69]. Charged tracks are reconstructed using
the Inner Tracking System (ITS) [70] and the Time Projection Chamber (TPC)
[71]. Tracks are selected with at least 70 TPC space points out of 159 points,
and a $\chi^{2}$ per TPC cluster less than $4$ [72]. To reduce contamination
from secondary particles, tracks are selected with a distance of closest
approach to the primary vertex of less than 2 cm in the longitudinal direction
and a $p_{\rm T}$-dependent selection ranging from 0.2 cm at $p_{\rm T}$ = 0.2
GeV/c to 0.018 cm at $p_{\rm T}$ = 4.0 GeV/c in the transverse direction.
Additionally, charged particle tracks are selected within the kinematic range
$|\eta|<0.8$ and $0.2<p_{\rm T}<5$ GeV/c. In order to suppress non-flow
correlation contributions, which are not related to collective behaviour, such
as resonance decays and jets, a gap in pseudorapidity is introduced between
the correlated particles. An $\eta$-gap of $|\Delta\eta|>0.8$ is used for two-
particle correlations and $|\Delta\eta|>0$ for four-particle correlations in
order to ensure optimal balance between the statistical precision and non-flow
suppression. To further investigate the non-flow suppression, the analysis was
also performed with only particles of same-sign charge and with increasing
pseudorapidity gaps between the subevents of correlated particles. Additional
Monte Carlo studies using HIJING [73], a model that does not feature
collective effects, but involves particle correlations arising from other
sources, indicate that non-flow is sufficiently suppressed with the applied
pseudorapidity gaps. Based on these model studies, the non-flow correlation
contributions to the correlators are less than 1% of the measured value across
the presented centralities.

Systematic uncertainties are evaluated by varying the event and track
selection criteria. Since the systematic uncertainty may depend on the
collision centrality and the $p_{\rm T}$ of particles, only the largest
contribution from each source is mentioned below. The systematic uncertainty
associated to the event selection criteria is evaluated by varying the
selection on the vertex position along the beam direction (from 10 cm to 7, 8,
or 9 cm), the magnetic field polarity, and the criteria for selecting pileup
events, and is below 1%. Possible biases from the centrality determination
were investigated by repeating the analysis using the centrality estimated at
midrapidity from hits in the Silicon Pixel Detector (SPD) instead of the V0,
resulting in a negligible difference. Uncertainties related to track selection
are estimated by considering different track reconstruction and track quality
selection criteria. Changing the track reconstruction to include tracks
without hits in the SPD leads to a variation of at most $1.7\%$. The quality
of the reconstructed tracks is varied by changing the minimum number of TPC
space points to $80$ and $90$, which leads to a difference of $1.5\%$ on the
measured correlators. Uncertainties related to the variations in the distance
of closest approach in both longitudinal and transverse directions are
negligible, indicating that the effect of contaminations from secondary
particles on the measurements is insignificant. Finally, tightening the
$\chi^{2}$ per number of TPC clusters from 4 to 2.5 gives an uncertainty of at
most $3\%$. The total systematic uncertainty is calculated as the quadratic
sum of the individual sources that have a statistically significant
contribution according to a statistical test [74].

The $A^{\rm f}_{2}$ measurements as a function of the associate particle
transverse momentum $p_{\rm T}^{\rm a}$ in centrality classes 0–5%, 10–20%,
and 30–40% are shown in Fig. 1. The results are presented from 0.2 GeV/c up to
4 GeV/c, since the requirement of two particles at high $p_{\rm T}^{\rm a}$
for the four-particle correlations limits the available data sample at higher
transverse momentum. In the most central collisions 0–5%, a large deviation
from unity is observed starting from $p_{\rm T}^{\rm a}$ $\approx$ 2.5 GeV/c,
which increases towards higher $p_{\rm T}^{\rm a}$ range. As previously
mentioned, this deviation cannot be attributed to non-flow effects, whose
contributions are negligible. With more than $5\sigma$ significance of the
deviation at high $p_{\rm T}^{\rm a}$ across the presented centralities, these
measurements provide the first observation of $p_{\rm T}$-dependent flow angle
fluctuations in Pb–Pb collisions at $\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV. In
centralities 10–20% and 30–40%, the fluctuation weakens and reaches around
5–7% deviation from unity at high $p_{\rm T}^{\rm a}$ $\approx 3.5$ GeV/c in
comparison to 0–5% most central collisions. The increasing deviation from
unity with $p_{\rm T}^{\rm a}$ above 3 GeV/c observed in data suggests that
the elliptic flow at large transverse momentum ($p_{\rm T}>2.5$ GeV/c) may not
be correlated with a common symmetry plane. This will affect the comparison of
measurements relying on a common symmetry plane between particles at high and
low $p_{\rm T}$ with theoretical models that do not feature flow angle
fluctuations, e.g. measurements of the flow of charmed hadrons, used to
constrain the path-length dependence of in-medium energy loss of the heavy-
quarks [53, 54] and their transport properties [75].

Figure 1: The flow angle fluctuation $A^{\rm f}_{2}$ in Pb–Pb collisions at
$\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV as a function of the associate particle
transverse momentum, $p_{\rm T}^{\rm a}$, in centrality classes 0–5% (left),
10–20% (middle), and 30–40% (right). Comparison with iEBE-VISHNU with $\rm
T_{R}ENTo$ initial conditions and $\eta/s(T)$ [76], iEBE-VISHNU with AMPT
initial conditions and $\eta/s=0.08$ [76] and $\eta/s=0.20$, CoLBT [77], and
AMPT [78] are shown as coloured bands.

Theoretical calculations from iEBE-VISHNU [79], MUSIC [80, 81] as well as AMPT
[82, 78] models are, when available, compared to the data in Fig. 1. The iEBE-
VISHNU model is a (2+1)D viscous hydrodynamical model coupled to the hadronic
cascade model UrQMD [83]. This model has been successful in describing
collective phenomena in different systems and event-by-event fluctuations [76,
84]. The iEBE-VISHNU model is used with $\rm T_{R}ENTo$ initial conditions and
has a temperature-dependent specific shear viscosity $\eta/s(T)$ [85]. The
model is limited to $p_{\rm T}$ ranges below 3 GeV/c. The MUSIC model is an
event-by-event (3+1D) viscous hydrodynamic model and is used with both Glauber
[86] and $\rm T_{R}ENTo$ [87] initial conditions. Different values of the
specific shear viscosity of 0.08, 0.12 and 0.16 are used with $\rm T_{R}ENTo$
initial conditions [88]. The AMPT transport model uses partonic interactions
with string melting, and a quark coalescence model is used to simulate the
formation of hadrons, which are then transported through a hadronic cascade
model [89]. It is tuned to measurements of $\mathrm{d}N/\mathrm{d}\eta$,
$p_{\rm T}$-spectra and elliptic flow of charged pions, kaons and protons from
the ALICE experiment [90, 78]. The iEBE-VISHNU calculation underestimates the
deviation of $A^{\rm f}_{2}$ from unity at above 2.5 GeV/c across the
presented centralities with the largest difference in the 0–5% most central
collisions. The iEBE-VISHNU model with $\rm T_{R}ENTo$ initial conditions uses
parameters extracted from a Bayesian analysis [85] in contrast to the MUSIC
models which use standard $\rm T_{R}ENTo$ initial conditions [87]. The
Bayesian $\rm T_{R}ENTo$ represents the current best understanding of the
initial conditions and QGP transport properties, so these measurements bring
further constraints on the state-of-the-art models. The comparison of the
MUSIC calculations with Glauber and $\rm T_{R}ENTo$ initial conditions shows
that $A^{\rm f}_{2}$ is sensitive to the fluctuations in the initial state
with little sensitivity to the different values of specific shear viscosity.
The AMPT calculation describes the data well in the 0–5% most central
collisions and captures the deviation from unity in the highest $p_{\rm
T}^{\rm a}$ bin. At higher centralities the AMPT calculation overestimates the
deviation from unity at high $p_{\rm T}^{\rm a}$.

Figure 2: The flow magnitude fluctuation $M^{\rm f}_{2}$ in Pb–Pb collisions
at $\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV as a function of the associate particle
transverse momentum, $p_{\rm T}^{\rm a}$, in centrality classes 0–5% (left),
10–20% (middle), and 30–40% (right). Comparison with iEBE-VISHNU with $\rm
T_{R}ENTo$ initial conditions and $\eta/s(T)$ [76], MUSIC with Glauber initial
conditions and $\eta/s=0.08$ [88], MUSIC with $\rm T_{R}ENTo$ initial
conditions and $\eta/s=0.08,0.12,0.16$ [88], and AMPT [78] are shown as
coloured bands.

The measurements of the $p_{\rm T}$-dependent fluctuations of the flow
magnitude, $M^{\rm f}_{2}$, as a function of the transverse momentum, $p_{\rm
T}^{\rm a}$, in centrality classes 0–5%, 10–20%, and 30–40% are shown in Fig.
2. A substantial deviation of $M^{\rm f}_{2}$ from unity is observed in the
0–5% most central collisions. This deviation by more than $5\sigma$
significance in the highest bin of $p_{\rm T}^{\rm a}$ reveals the first
observation of $p_{\rm T}$-dependent flow magnitude fluctuations in Pb–Pb
collisions at $\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV. In centrality 0–5%, where
the initial state fluctuations are most significant, the flow magnitude
deviates from unity for $p_{\rm T}^{\rm a}$ above 2 GeV/c and the deviation
increases with increasing $p_{\rm T}^{\rm a}$. The deviation from unity is
more pronounced in the central collisions 0–5% compared to 10–20% and 30–40%.
By construction, $M^{\rm f}_{2}$ is not restricted to be below unity as seen
in Fig. 2 at 30–40% centrality. This is in contrast to the flow vector
fluctuations measured with two-particle correlations and the flow angle
fluctuations measured with $A^{\rm f}_{2}$, which can only be larger than
unity due to non-flow effects [56]. Model studies with HIJING and the studies
of the correlators with same-sign charge and larger pseudorapidity gaps show
that non-flow correlations are negligible for $M^{\rm f}_{2}$, so the
deviation from unity must be due to the fluctuations of the flow magnitude.
The iEBE-VISHNU calculations underestimate the effect in 0–5% centrality
showing almost no $p_{\rm T}$ dependence. In the 10–20% and 30–40% centrality
intervals, the model captures the increasing trend with $p_{\rm T}^{\rm a}$
and, consequently, describes the data. The MUSIC models show strong
sensitivity to the specific shear viscosity in 0–5% most central collisions,
and also a sensitivity to the different initial conditions. In the 10–20% and
30–40% centrality intervals, $M^{\rm f}_{2}$ shows no sensitivity to the
specific shear viscosity, but is still affected by the different initial
conditions. The MUSIC models overestimate the deviation of $M^{\rm f}_{2}$
from unity. The comparison of the measurements with these models reveals that
the flow magnitude fluctuations are driven by initial state fluctuations, and
that in most central collisions they are also affected by the transport
properties of the QGP. The AMPT transport model calculation succeeds in
describing the flow magnitude fluctuations in the most central collisions,
where it also describes the flow angle fluctuations (Fig. 1 left). At higher
centralities, the AMPT model significantly overestimates the data even at low
$p_{\rm T}^{\rm a}$, but it qualitatively captures the increasing trend of
$M^{\rm f}_{2}$ in the data in the 30–40% centrality interval.

Figure 3: The lower limit of the first-order flow angle fluctuations, upper
limit of the first-order flow magnitude fluctuations, and the flow vector
fluctuations as a function of centrality for the $3.0<p_{\rm T}^{\rm a}<4.0$
GeV/c. The lower and upper limits are denoted by the arrows. The top panel
shows the absolute contribution and the bottom panel the contribution relative
to the overall flow vector fluctuations. Comparisons to AMPT [78] are shown as
coloured bands.

In order to connect the four-particle correlation observables to the
previously measured flow vector fluctuations, the limits from Eqs. (4) and (6)
are compared to the flow vector fluctuations probed by Eq. (5). Thus the lower
moments of the flow angle fluctuations and flow magnitude fluctuations, which
cannot be directly accessed in experiments, can be explored. The lower limit
on the first-order flow angle fluctuations, the upper limit on the first-order
flow magnitude fluctuations, and the total flow vector fluctuations are shown
as a function of centrality in the $3<p_{\rm T}^{\rm a}<4$ GeV/c range in Fig.
3. In central collisions, the upper limit on the flow magnitude fluctuations
is higher than the lower limit of the flow angle fluctuations up to 10%
centrality. For the 10–30% centrality interval, the limits are similar, and
above 30% centrality, the flow magnitude upper limits are much smaller, and
the flow angle fluctuations dominate the overall flow vector fluctuations
completely. This is consistent with the measurements shown in Figs. 1 and 2,
where $M^{\rm f}_{2}$ approaches unity faster with increasing centrality than
$A^{\rm f}_{2}$ and even goes above unity in semicentral collisions. While
$p_{\rm T}$-dependent flow angle and flow magnitude fluctuations are both
present in central Pb–Pb collisions, the effects of flow angle fluctuations
are consistently present across the shown centralities compared to the much
smaller flow magnitude fluctuations in non-central collisions. The AMPT
transport model calculations overestimate the flow vector fluctuations
$v_{2}\\{2\\}/v_{2}[2]$ as well as the limits of the first-order flow angle
and flow magnitude fluctuations in the central collisions. At higher
centralities, AMPT describes the lower and upper limits well, whereas the flow
vector fluctuations are still overestimated by the model as in central
collisions.

In summary, the $p_{\rm T}$-dependent flow angle fluctuations and flow
magnitude fluctuations of the second-order flow vector $\vec{V}_{2}$ are
measured in Pb–Pb collisions at $\sqrt{s_{{}_{\rm NN}}}$ = 5.02 TeV with the
correlators $A^{\rm f}_{2}$ and $M^{\rm f}_{2}$ as function of the associate
particle momentum in centrality intervals 0–5%, 10–20%, and 30–40%. Large
deviations from unity of both $A^{\rm f}_{2}$ ($\sim$20%) and $M^{\rm f}_{2}$
($\sim$10%) are observed in central collisions, where the event-by-event
fluctuations of the position of the colliding nucleons dominate over the
geometric response. In semicentral collisions, both flow angle and flow
magnitude fluctuations decrease. The flow angle fluctuations reach 5–7% and
the flow magnitude fluctuations reach around 2% and are above unity in 30–40%
centrality. The flow magnitude fluctuations decrease faster than the flow
angle fluctuations up to 40% centrality, where the flow vector fluctuations
are almost solely due to flow angle fluctuations. The proposed correlators
decompose the flow vector fluctuations into constituent flow angle and flow
magnitude fluctuations. The observation of flow angle and flow magnitude
fluctuations gives further insight into the nature of event-by-event
fluctuations in the initial state of heavy-ion collisions and comparison with
theoretical models provides further constraints on the initial conditions and
even the QGP transport coefficient $\eta/s$ in central collisions. The flow
angle fluctuations and flow magnitude fluctuations are most pronounced at high
$p_{\rm T}$ and will affect measurements that assume a common symmetry plane,
such as high-$p_{\rm T}$ azimuthal correlations, which is the typical method
for high-$p_{\rm T}$ flow measurements. The measurements provide further
constraints on existing models, which are necessary to avoid biasing the
extraction of QGP transport properties, such as the jet transport coefficient,
$\hat{q}$, with flow measurements.

## References

  * [1] E. V. Shuryak, “Quark-Gluon Plasma and Hadronic Production of Leptons, Photons and Psions”, Phys. Lett. B78 (1978) 150. [Yad. Fiz.28,796(1978)].
  * [2] E. V. Shuryak, “Quantum Chromodynamics and the Theory of Superdense Matter”, Phys. Rept. 61 (1980) 71–158.
  * [3] BRAHMS Collaboration, I. Arsene et al., “Quark gluon plasma and color glass condensate at RHIC? The Perspective from the BRAHMS experiment”, Nucl. Phys. A 757 (2005) 1–27, arXiv:nucl-ex/0410020.
  * [4] STAR Collaboration, J. Adams et al., “Experimental and theoretical challenges in the search for the quark gluon plasma: The STAR Collaboration’s critical assessment of the evidence from RHIC collisions”, Nucl. Phys. A 757 (2005) 102–183, arXiv:nucl-ex/0501009.
  * [5] PHENIX Collaboration, K. Adcox et al., “Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration”, Nucl. Phys. A 757 (2005) 184–283, arXiv:nucl-ex/0410003.
  * [6] PHOBOS Collaboration, B. B. Back et al., “The PHOBOS perspective on discoveries at RHIC”, Nucl. Phys. A 757 (2005) 28–101, arXiv:nucl-ex/0410022.
  * [7] B. Muller, J. Schukraft, and B. Wyslouch, “First Results from Pb+Pb collisions at the LHC”, Ann. Rev. Nucl. Part. Sci. 62 (2012) 361–386, arXiv:1202.3233 [hep-ex].
  * [8] J.-Y. Ollitrault, “Anisotropy as a signature of transverse collective flow”, Phys. Rev. D46 (1992) 229–245.
  * [9] S. A. Voloshin, A. M. Poskanzer, and R. Snellings, “Collective phenomena in non-central nuclear collisions”, Landolt-Bornstein 23 (2010) 293–333, arXiv:0809.2949 [nucl-ex].
  * [10] U. Heinz and R. Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions”, Ann. Rev. Nucl. Part. Sci. 63 (2013) 123–151, arXiv:1301.2826 [nucl-th].
  * [11] S. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by Fourier expansion of Azimuthal particle distributions”, Z. Phys. C70 (1996) 665–672, arXiv:hep-ph/9407282 [hep-ph].
  * [12] B. Alver and G. Roland, “Collision geometry fluctuations and triangular flow in heavy-ion collisions”, Phys. Rev. C81 (2010) 054905, arXiv:1003.0194 [nucl-th]. [Erratum: Phys. Rev.C82,039903(2010)].
  * [13] B. H. Alver, C. Gombeaud, M. Luzum, and J.-Y. Ollitrault, “Triangular flow in hydrodynamics and transport theory”, Phys. Rev. C82 (2010) 034913, arXiv:1007.5469 [nucl-th].
  * [14] D. Teaney and L. Yan, “Triangularity and Dipole Asymmetry in Heavy Ion Collisions”, Phys. Rev. C83 (2011) 064904, arXiv:1010.1876 [nucl-th].
  * [15] M. Luzum, “Collective flow and long-range correlations in relativistic heavy ion collisions”, Phys. Lett. B696 (2011) 499–504, arXiv:1011.5773 [nucl-th].
  * [16] STAR Collaboration, K. H. Ackermann et al., “Elliptic flow in Au + Au collisions at (S(NN))**(1/2) = 130 GeV”, Phys. Rev. Lett. 86 (2001) 402–407, arXiv:nucl-ex/0009011.
  * [17] PHENIX Collaboration, S. S. Adler et al., “Elliptic flow of identified hadrons in Au+Au collisions at s(NN)**(1/2) = 200-GeV”, Phys. Rev. Lett. 91 (2003) 182301, arXiv:nucl-ex/0305013.
  * [18] STAR Collaboration, L. Adamczyk et al., “Third Harmonic Flow of Charged Particles in Au+Au Collisions at sqrtsNN = 200 GeV”, Phys. Rev. C 88 no. 1, (2013) 014904, arXiv:1301.2187 [nucl-ex].
  * [19] PHENIX Collaboration, A. Adare et al., “Measurements of elliptic and triangular flow in high-multiplicity 3He$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV”, Phys. Rev. Lett. 115 no. 14, (2015) 142301, arXiv:1507.06273 [nucl-ex].
  * [20] ALICE Collaboration, K. Aamodt et al., “Elliptic flow of charged particles in Pb-Pb collisions at 2.76 TeV”, Phys. Rev. Lett. 105 (2010) 252302, arXiv:1011.3914 [nucl-ex].
  * [21] ALICE Collaboration, K. Aamodt et al., “Higher harmonic anisotropic flow measurements of charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV”, Phys. Rev. Lett. 107 (2011) 032301, arXiv:1105.3865 [nucl-ex].
  * [22] ALICE Collaboration, B. Abelev et al., “Elliptic flow of identified hadrons in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV”, JHEP 06 (2015) 190, arXiv:1405.4632 [nucl-ex].
  * [23] ALICE Collaboration, J. Adam et al., “Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Rev. Lett. 116 no. 13, (2016) 132302, arXiv:1602.01119 [nucl-ex].
  * [24] ALICE Collaboration, S. Acharya et al., “Linear and non-linear flow modes in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=$ 2.76 TeV”, Phys. Lett. B773 (2017) 68–80, arXiv:1705.04377 [nucl-ex].
  * [25] ATLAS Collaboration, G. Aad et al., “Measurement of the azimuthal anisotropy for charged particle production in $\sqrt{s_{NN}}=2.76$ TeV lead-lead collisions with the ATLAS detector”, Phys. Rev. C86 (2012) 014907, arXiv:1203.3087 [hep-ex].
  * [26] ATLAS Collaboration, G. Aad et al., “Measurement of the pseudorapidity and transverse momentum dependence of the elliptic flow of charged particles in lead-lead collisions at $\sqrt{s_{NN}}=2.76$ TeV with the ATLAS detector”, Phys. Lett. B707 (2012) 330–348, arXiv:1108.6018 [hep-ex].
  * [27] ATLAS Collaboration, G. Aad et al., “Measurement of the distributions of event-by-event flow harmonics in lead-lead collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV with the ATLAS detector at the LHC”, JHEP 11 (2013) 183, arXiv:1305.2942 [hep-ex].
  * [28] CMS Collaboration, S. Chatrchyan et al., “Centrality dependence of dihadron correlations and azimuthal anisotropy harmonics in PbPb collisions at $\sqrt{s_{NN}}=2.76$ TeV”, Eur. Phys. J. C72 (2012) 2012, arXiv:1201.3158 [nucl-ex].
  * [29] CMS Collaboration, S. Chatrchyan et al., “Measurement of the elliptic anisotropy of charged particles produced in PbPb collisions at $\sqrt{s}_{NN}$=2.76 TeV”, Phys. Rev. C87 no. 1, (2013) 014902, arXiv:1204.1409 [nucl-ex].
  * [30] CMS Collaboration, S. Chatrchyan et al., “Azimuthal anisotropy of charged particles at high transverse momenta in PbPb collisions at $\sqrt{s_{NN}}=2.76$ TeV”, Phys. Rev. Lett. 109 (2012) 022301, arXiv:1204.1850 [nucl-ex].
  * [31] M. Luzum and H. Petersen, “Initial State Fluctuations and Final State Correlations in Relativistic Heavy-Ion Collisions”, J. Phys. G41 (2014) 063102, arXiv:1312.5503 [nucl-th].
  * [32] E. Shuryak, “Strongly coupled quark-gluon plasma in heavy ion collisions”, Rev. Mod. Phys. 89 (2017) 035001, arXiv:1412.8393 [hep-ph].
  * [33] H. Song, Y. Zhou, and K. Gajdosova, “Collective flow and hydrodynamics in large and small systems at the LHC”, Nucl. Sci. Tech. 28 no. 7, (2017) 99, arXiv:1703.00670 [nucl-th].
  * [34] PHENIX Collaboration, A. Adare et al., “Azimuthal anisotropy of neutral pion production in Au+Au collisions at $\sqrt{(}s_{\mathrm{NN}})$ = 200 GeV: Path-length dependence of jet quenching and the role of initial geometry”, Phys. Rev. Lett. 105 (2010) 142301, arXiv:1006.3740 [nucl-ex].
  * [35] ALICE Collaboration, B. Abelev et al., “Anisotropic flow of charged hadrons, pions and (anti-)protons measured at high transverse momentum in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV”, Phys. Lett. B719 (2013) 18–28, arXiv:1205.5761 [nucl-ex].
  * [36] CMS Collaboration, S. Chatrchyan et al., “Azimuthal anisotropy of charged particles at high transverse momenta in PbPb collisions at $\sqrt{s_{NN}}=2.76$ TeV”, Phys. Rev. Lett. 109 (2012) 022301, arXiv:1204.1850 [nucl-ex].
  * [37] ATLAS Collaboration, G. Aad et al., “Measurement of flow harmonics with multi-particle cumulants in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV with the ATLAS detector”, Eur. Phys. J. C74 no. 11, (2014) 3157, arXiv:1408.4342 [hep-ex].
  * [38] ALICE Collaboration, J. Adam et al., “Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Rev. Lett. 116 no. 13, (2016) 132302, arXiv:1602.01119 [nucl-ex].
  * [39] CMS Collaboration, A. M. Sirunyan et al., “Azimuthal anisotropy of charged particles with transverse momentum up to 100 GeV/ c in PbPb collisions at $\sqrt{s}_{{NN}}$=5.02 TeV”, Phys. Lett. B776 (2018) 195–216, arXiv:1702.00630 [hep-ex].
  * [40] PHENIX Collaboration, A. Adare et al., “Heavy Quark Production in $p+p$ and Energy Loss and Flow of Heavy Quarks in Au+Au Collisions at $\sqrt{s_{NN}}=200$ GeV”, Phys. Rev. C84 (2011) 044905, arXiv:1005.1627 [nucl-ex].
  * [41] STAR Collaboration, L. Adamczyk et al., “Elliptic flow of electrons from heavy-flavor hadron decays in Au + Au collisions at $\sqrt{s_{\rm NN}}=$ 200, 62.4, and 39 GeV”, Phys. Rev. C95 no. 3, (2017) 034907, arXiv:1405.6348 [hep-ex].
  * [42] STAR Collaboration, L. Adamczyk et al., “Measurement of $D^{0}$ Azimuthal Anisotropy at Midrapidity in Au+Au Collisions at $\sqrt{s_{NN}}$=200 GeV”, Phys. Rev. Lett. 118 no. 21, (2017) 212301, arXiv:1701.06060 [nucl-ex].
  * [43] ALICE Collaboration, B. Abelev et al., “D meson elliptic flow in non-central Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76TeV”, Phys. Rev. Lett. 111 (2013) 102301, arXiv:1305.2707 [nucl-ex].
  * [44] ALICE Collaboration, E. Abbas et al., “J/$\psi$ Elliptic Flow in Pb–Pb Collisions at $\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV”, Phys. Rev. Lett. 111 (2013) 162301, arXiv:1303.5880 [nucl-ex].
  * [45] ALICE Collaboration, B. Abelev et al., “Azimuthal anisotropy of D meson production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV”, Phys. Rev. C 90 no. 3, (2014) 034904, arXiv:1405.2001 [nucl-ex].
  * [46] ALICE Collaboration, J. Adam et al., “Elliptic flow of electrons from heavy-flavour hadron decays at mid-rapidity in Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=2.76$ TeV”, JHEP 09 (2016) 028, arXiv:1606.00321 [nucl-ex].
  * [47] ALICE Collaboration, J. Adam et al., “Elliptic flow of muons from heavy-flavour hadron decays at forward rapidity in Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV”, Phys. Lett. B753 (2016) 41–56, arXiv:1507.03134 [nucl-ex].
  * [48] ALICE Collaboration, S. Acharya et al., “$D$-meson azimuthal anisotropy in midcentral Pb-Pb collisions at $\mathbf{\sqrt{s_{\rm NN}}=5.02}$ TeV”, Phys. Rev. Lett. 120 no. 10, (2018) 102301, arXiv:1707.01005 [nucl-ex].
  * [49] CMS Collaboration, A. M. Sirunyan et al., “Measurement of prompt D0 meson azimuthal anisotropy in Pb-Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}$ = 5.02 TeV”, Phys. Rev. Lett. 120 no. 20, (2018) 202301, arXiv:1708.03497 [nucl-ex].
  * [50] ALICE Collaboration, S. Acharya et al., “Transverse-momentum and event-shape dependence of D-meson flow harmonics in Pb–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV”, Phys. Lett. B813 (2021) 136054, arXiv:2005.11131 [nucl-ex].
  * [51] ALICE Collaboration, S. Acharya et al., “Prompt D0, D+, and D∗+ production in Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV”, JHEP 01 (2022) 174, arXiv:2110.09420 [nucl-ex].
  * [52] ALICE Collaboration, S. Acharya et al., “Measurement of prompt $D_{s}^{+}$-meson production and azimuthal anisotropy in Pb–Pb collisions at $\sqrt{s_{NN}}$=5.02TeV”, Phys. Lett. B827 (2022) 136986, arXiv:2110.10006 [nucl-ex].
  * [53] M. Gyulassy, I. Vitev, and X. N. Wang, “High p(T) azimuthal asymmetry in noncentral A+A at RHIC”, Phys. Rev. Lett. 86 (2001) 2537–2540, arXiv:nucl-th/0012092.
  * [54] E. V. Shuryak, “The Azimuthal asymmetry at large p(t) seem to be too large for a ‘jet quenching”’, Phys. Rev. C66 (2002) 027902, arXiv:nucl-th/0112042.
  * [55] U. Heinz, Z. Qiu, and C. Shen, “Fluctuating flow angles and anisotropic flow measurements”, Phys. Rev. C87 no. 3, (2013) 034913, arXiv:1302.3535 [nucl-th].
  * [56] F. G. Gardim, F. Grassi, M. Luzum, and J.-Y. Ollitrault, “Breaking of factorization of two-particle correlations in hydrodynamics”, Phys. Rev. C87 no. 3, (2013) 031901, arXiv:1211.0989 [nucl-th].
  * [57] ALICE Collaboration, S. Acharya et al., “Searches for transverse momentum dependent flow vector fluctuations in Pb-Pb and p-Pb collisions at the LHC”, JHEP 09 (2017) 032, arXiv:1707.05690 [nucl-ex].
  * [58] CMS Collaboration, S. Chatrchyan et al., “Studies of azimuthal dihadron correlations in ultra-central PbPb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV”, JHEP 02 (2014) 088, arXiv:1312.1845 [nucl-ex].
  * [59] CMS Collaboration, V. Khachatryan et al., “Evidence for transverse momentum and pseudorapidity dependent event plane fluctuations in PbPb and pPb collisions”, Phys. Rev. C92 no. 3, (2015) 034911, arXiv:1503.01692 [nucl-ex].
  * [60] ALICE Collaboration, K. Aamodt et al., “Harmonic decomposition of two-particle angular correlations in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 2.76 TeV”, Phys. Lett. B708 (2012) 249–264, arXiv:1109.2501 [nucl-ex].
  * [61] R. S. Bhalerao, J.-Y. Ollitrault, S. Pal, and D. Teaney, “Principal component analysis of event-by-event fluctuations”, Phys. Rev. Lett. 114 no. 15, (2015) 152301, arXiv:1410.7739 [nucl-th].
  * [62] Z. Liu, W. Zhao, and H. Song, “Principal Component Analysis of collective flow in Relativistic Heavy-Ion Collisions”, Eur. Phys. J. C79 no. 10, (2019) 870, arXiv:1903.09833 [nucl-th].
  * [63] Z. Liu, A. Behera, H. Song, and J. Jia, “Robustness of principal component analysis of harmonic flow in heavy ion collisions”, Phys. Rev. C102 no. 2, (2020) 024911, arXiv:2002.06061 [nucl-ex].
  * [64] M. Hippert, D. Dobrigkeit Chinellato, M. Luzum, J. Noronha, T. Nunes da Silva, and J. Takahashi, “Measuring Momentum-Dependent Flow Fluctuations in Heavy-Ion Collisions”, Phys. Rev. C101 no. 3, (2020) 034903, arXiv:1906.08915 [nucl-th].
  * [65] P. Bozek and W. Broniowski, “Longitudinal decorrelation measures of flow magnitude and event-plane angles in ultrarelativistic nuclear collisions”, Phys. Rev. C 97 no. 3, (2018) 034913, arXiv:1711.03325 [nucl-th].
  * [66] A. Bilandzic, C. H. Christensen, K. Gulbrandsen, A. Hansen, and Y. Zhou, “Generic framework for anisotropic flow analyses with multiparticle azimuthal correlations”, Phys. Rev. C89 no. 6, (2014) 064904, arXiv:1312.3572 [nucl-ex].
  * [67] Z. Moravcova, K. Gulbrandsen, and Y. Zhou, “Generic algorithm for multiparticle cumulants of azimuthal correlations in high energy nucleus collisions”, Phys. Rev. C103 no. 2, (2021) 024913, arXiv:2005.07974 [nucl-th].
  * [68] ALICE Collaboration, B. Abelev et al., “Performance of the ALICE Experiment at the CERN LHC”, Int. J. Mod. Phys. A29 (2014) 1430044, arXiv:1402.4476 [nucl-ex].
  * [69] ALICE Collaboration, B. Abelev et al., “Centrality determination of Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV with ALICE”, Phys. Rev. C88 no. 4, (2013) 044909, arXiv:1301.4361 [nucl-ex].
  * [70] ALICE Collaboration, K. Aamodt et al., “Alignment of the ALICE Inner Tracking System with cosmic-ray tracks”, JINST 5 (2010) P03003, arXiv:1001.0502 [physics.ins-det].
  * [71] J. Alme et al., “The ALICE TPC, a large 3-dimensional tracking device with fast readout for ultra-high multiplicity events”, Nucl. Instrum. Meth. A622 (2010) 316–367, arXiv:1001.1950 [physics.ins-det].
  * [72] ALICE Collaboration, “The ALICE definition of primary particles”, ALICE-PUBLIC-2017-005 (Jun, 2017) . https://cds.cern.ch/record/2270008.
  * [73] X.-N. Wang and M. Gyulassy, “hijing: A monte carlo model for multiple jet production in $\mathrm{pp}$, $\mathrm{pA}$, and $\mathrm{AA}$ collisions”, Phys. Rev. D44 (Dec, 1991) 3501–3516.
  * [74] R. Barlow, “Systematic errors: Facts and fictions”, in Conference on Advanced Statistical Techniques in Particle Physics, pp. 134–144.  7, 2002.  arXiv:hep-ex/0207026.
  * [75] W. Ke, Y. Xu, and S. A. Bass, “Linearized Boltzmann-Langevin model for heavy quark transport in hot and dense QCD matter”, Phys. Rev. C98 no. 6, (2018) 064901, arXiv:1806.08848 [nucl-th].
  * [76] W. Zhao, H.-j. Xu, and H. Song, “Collective flow in 2.76 A TeV and 5.02 A TeV Pb+Pb collisions”, Eur. Phys. J. C77 no. 9, (2017) 645, arXiv:1703.10792 [nucl-th].
  * [77] W. Zhao, W. Ke, W. Chen, T. Luo, and X.-N. Wang, “From Hydrodynamics to Jet Quenching, Coalescence, and Hadron Cascade: A Coupled Approach to Solving the RAA$\otimes$v2 Puzzle”, Phys. Rev. Lett. 128 no. 2, (2022) 022302, arXiv:2103.14657 [hep-ph].
  * [78] G.-L. Ma and Z.-W. Lin, “Predictions for $\sqrt{s_{NN}}=5.02$ TeV Pb+Pb Collisions from a Multi-Phase Transport Model”, Phys. Rev. C93 no. 5, (2016) 054911, arXiv:1601.08160 [nucl-th].
  * [79] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions”, Comput. Phys. Commun. 199 (2016) 61–85, arXiv:1409.8164 [nucl-th].
  * [80] B. Schenke, S. Jeon, and C. Gale, “Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics”, Phys. Rev. Lett. 106 (2011) 042301, arXiv:1009.3244 [hep-ph].
  * [81] B. Schenke, S. Jeon, and C. Gale, “(3+1)D hydrodynamic simulation of relativistic heavy-ion collisions”, Phys. Rev. C82 (2010) 014903, arXiv:1004.1408 [hep-ph].
  * [82] Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, “A Multi-phase transport model for relativistic heavy ion collisions”, Phys. Rev. C72 (2005) 064901, arXiv:nucl-th/0411110 [nucl-th].
  * [83] S. A. Bass et al., “Microscopic models for ultrarelativistic heavy ion collisions”, Prog. Part. Nucl. Phys. 41 (1998) 255–369, arXiv:nucl-th/9803035 [nucl-th].
  * [84] M. Li, Y. Zhou, W. Zhao, B. Fu, Y. Mou, and H. Song, “Investigations on mixed harmonic cumulants in heavy-ion collisions at energies available at the CERN Large Hadron Collider”, Phys. Rev. C104 no. 2, (2021) 024903, arXiv:2104.10422 [nucl-th].
  * [85] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, “Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium”, Phys. Rev. C94 no. 2, (2016) 024907, arXiv:1605.03954 [nucl-th].
  * [86] P. Bożek, W. Broniowski, M. Rybczynski, and G. Stefanek, “GLISSANDO 3: GLauber Initial-State Simulation AND mOre…, ver. 3”, Comput. Phys. Commun. 245 (2019) 106850, arXiv:1901.04484 [nucl-th].
  * [87] J. S. Moreland, J. E. Bernhard, and S. A. Bass, “Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions”, Phys. Rev. C92 no. 1, (2015) 011901, arXiv:1412.4708 [nucl-th].
  * [88] P. Bozek and R. Samanta, “Factorization breaking for higher moments of harmonic flow”, arXiv:2109.07781 [nucl-th].
  * [89] B.-A. Li and C. M. Ko, “Formation of superdense hadronic matter in high-energy heavy ion collisions”, Phys. Rev. C52 (1995) 2037–2063, arXiv:nucl-th/9505016 [nucl-th].
  * [90] J. Xu and C. M. Ko, “Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV in a multiphase transport model”, Phys. Rev. C 83 (2011) 034904, arXiv:1101.2231 [nucl-th].