Principal Component Analysis of collective flow in Relativistic Heavy-Ion Collisions

TL;DR

This study applies Principal Component Analysis to hydrodynamic simulation data from heavy-ion collisions to identify flow patterns without human bias, revealing similarities and differences with traditional Fourier-based methods.

Contribution

It demonstrates that PCA can effectively extract flow harmonics from collision data, offering a data-driven alternative to traditional Fourier analysis.

Findings

PCA eigenvectors resemble Fourier bases but are not identical.

Flow harmonics $v_n^ ext{'}$ are similar to $v_n$ for n=2,3, but differ for n≥4.

Mode-coupling effects are reduced in PCA-defined flow harmonics.

Abstract

In this paper, we implement Principal Component Analysis (PCA) to study the single particle distributions generated from thousands of {\tt VISH2+1} hydrodynamic simulations with an aim to explore if a machine could directly discover flow from the huge amount of data without explicit instructions from human-beings. We found that the obtained PCA eigenvectors are similar to but not identical with the traditional Fourier bases. Correspondingly, the PCA defined flow harmonics $v_n^\prime$ are also similar to the traditional $v_n$ for $n=2$ and 3, but largely deviated from the Fourier ones for $n\geq 4$. A further study on the symmetric cumulants and the Pearson coefficients indicates that mode-coupling effects are reduced for these flow harmonics defined by PCA.

Tables

Table 1. Table 1 : Event averaged flow harmonics v n ′ superscript subscript 𝑣 𝑛 ′ v_{n}^{\prime} from PCA and v n subscript 𝑣 𝑛 v_{n} from the Fourier expansion, for VISH2+1 simulated Pb+Pb collisions at 10-20% centrality.

$n$	$\overline{v_n^{\prime }}$(PCA)	$\overline{v_n^{\prime }}\times {10}^2$	$\overline{v_n}\times {10}^2$
2	$\sqrt{\frac{m}{2}}\sqrt{{\sigma }_1^2+{\sigma }_2^2}$	6.03	6.08
3	$\sqrt{\frac{m}{2}}\sqrt{{\sigma }_3^2+{\sigma }_4^2}$	2.57	2.53
4	$\sqrt{\frac{m}{2}}\sqrt{{\sigma }_5^2+{\sigma }_6^2}$	1.21	1.25
5	$\sqrt{\frac{m}{2}}\sqrt{{\sigma }_9^2+{\sigma }_{10}^2}$	0.57	0.66
6	$\sqrt{\frac{m}{2}}\sqrt{{\sigma }_{11}^2+{\sigma }_{12}^2}$	0.26	0.37

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}\phi }$$

$$=\frac{1}{2\pi }(1+2\sum _{n=1}^{\infty }{v}_n{e}^{-in(\phi -{\Psi }_n)})$$

$$\frac{dN}{dy{p}_{T}d{p}_{T}d\phi }={\int }_{\Sigma }\frac{g}{(2\pi {)}^3}{p}^{\mu }{d}^3{\sigma }_{\mu }f(x,p)$$

$${\mathbf{M}}_{\mathbf{f}}=\mathbf{X}\mathbf{\Sigma }\mathbf{Z}=\mathbf{V}\mathbf{Z}$$

$$dN/d{\phi }^{(i)}$$

$$S{C}^v{}^{\prime }(m,n)$$

$$r(v_n^{\prime },{\epsilon }_m)=\frac{⟨v_n^{\prime }{\epsilon }_m⟩-⟨v_n^{\prime }⟩⟨{\epsilon }_m⟩}{\sqrt{(v_n^{\prime }-⟨v_n^{\prime }⟩{)}^2({\epsilon }_m-⟨{\epsilon }_m⟩{)}^2}}$$

$${\epsilon }_n{e}^{in{\Phi }_n}=-\frac{\int rdrd\phi r^n{e}^{in\phi }s(r,\phi )}{\int rdrd\phi r^{n}s(r,\phi )},$$

$$\frac{dS}{d\phi }=\int r^{2}drs(r,\phi )$$

$$(\frac{dS}{d\phi }{)}_{smooth}={\int }_{-\pi }^{\pi }K({\phi }^{\prime },\phi )\frac{dS}{d{\phi }^{\prime }}d{\phi }^{\prime }$$

$$dS/d{\phi }^{(i)}$$

Full text

††journal: Physics Letters B

1 Introduction

Collective flow is one of the most important observables in relativistic heavy-ion collisions, which provides valuable information on the initial state fluctuations, final state correlations and the QGP properties. In the past decades, various flow observables have been extensively measured in experiments and studied in theory [1, 2, 3, 4, 5, 6]. In general, these flow observables are defined based on the Fourier decomposition. For example, the integrated flow harmonics are defined as:

[TABLE]

where $\vec{V}_{n}=v_{n}e^{in\Psi_{n}}$ is the $n$-th order flow-vector, $v_{n}$ is the $n$-th order flow harmonics and $\Psi_{n}$ is the corresponding event plane angle. In general, the first coefficient, $v_{1}$, is called the directed flow, the second coefficient, $v_{2}$, is called the elliptic flow and the third coefficient $v_{3}$, is called the triangular flow. For $n\geq 3$, $v_{n}$ is also referred as the higher order flow harmonics.

In spite of the success of the flow measurements and the hydrodynamic descriptions, one essential question is why the Fourier expansion is a natural way to analyze the flow data. In this paper, we will address these questions with one of the machine learning techniques, called the Principal Component Analysis (PCA). In more details, we will investigate if a machine could directly discover flow from the huge amount of data of the relativistic fluid systems without explicit instructions from human beings.

PCA is one of the unsupervised algorithms of machine learning [7] based on the Singular Value Decomposition (SVD) that diagonalize a random matrix with two orthogonal matrices. Compared with other deep learning algorithms, the advantage of PCA lies in its simple and elegant mathematical formulation, which is understandable and traceable to human beings, and is able to reveal the main structure of data in a quite transparent way.

Due to its strong power in data mining, PCA has been implemented to various research area of physics [8, 9, 10, 11, 12, 13]. In molecular dynamics, PCA has been utilized to distinguish break junction trajectories of single molecules [8], which is time efficient and can transfer to a wide range of multivariate data sets. In the field of quantum mechanics, the quantum version of PCA was applied to study quantum coherence among different copies of the system [9], which are exponentially faster than any existing algorithm. In condensed matter physics, PCA has been implemented to study the phase transition in Ising model [11], which found that eigenvectors of PCA can aid in the definition of the order parameter, as well as provide reasonable predictions for the critical temperature without any prior knowledge. Besides, PCA is a widely used tool in engineering for model reduction to make computations more efficient [14].

In relativistic heavy-ion collisions, PCA has been implemented to study the event-by-event flow fluctuations, using the 2-particle correlations with the Fourier expansion [13, 15, 16, 17, 18]. Compared with the traditional method, PCA explores all the information contained in the 2-particle correlations, which reveals the substructures in flow fluctuations [13, 15, 16]. It was found that the leading components of PCA correspond to the traditional flow harmonics and the sub-leading components evaluate the breakdown of the flow factorization at different $p_{t}$ or $\eta$ bins. Besides, PCA has also been used to study the non-linear mode coupling between different flow harmonics[17], which helps to discover some hidden mode-mixing patterns. Recently, the CMS Collaboration further implemented PCA to analyze 2-particle correlation in Pb-Pb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV and p-Pb collisions at $\sqrt{s_{NN}}=$ 5.02 TeV [18], showing the potential of largely implementing such machine learning technique to realistic data in relativistic heavy ion collisions.

These early PCA investigations on flow [13, 15, 16, 17, 18] are all based on the preprocessed data with the Fourier expansion, which still belong to the category of traditional flow analysis. In this paper, we will directly apply PCA to study the single particle distributions from hydrodynamic simulations without any priori Fourier transformation. We aim to explore if PCA could discover flow with its own bases.

This paper is organized as follows. Sec. II introduces relativistic hydrodynamics, principal component analysis (PCA) and the corresponding flow analysis. Sec. III shows and discusses the flow results from PCA and compares them with the ones from traditional Fourier expansion. Sec. IV summarizes and concludes the paper.

2 Model and method

2.1 VISH2+1 hydrodynamics

In this paper, we implement VISH2+1 [19, 20, 21, 22] to generate the final particle distributions for the PCA analysis. VISH2+1 [19, 20, 21, 22] is a 2+1-dimensional viscous hydrodynamic code to simulate the expansion of the QGP fireballs, which solves the transport equations for the energy-momentum tenor $T^{\mu\nu}$ and the second order Israel-Stewart equations for the shear stress tensor $\pi^{\mu\nu}$ and bulk pressure $\Pi$ with an equation of state s95-PCE[23, 24] as an input. The initial profiles for VISH2+1 are provided by TRENTo, a parameterized initial condition model that generates event-by-event fluctuating entropy profiles with several tunable parameters [25, 24]. These parameters, together with the temperature dependent specific shear viscosity and bulk viscosity, hydrodynamic starting time ($\tau_{0}=0.6\ \mathrm{fm/c}$) and decoupling /switching temperature ($T_{sw}=148\ \mathrm{MeV}$) have been fixed through fitting all charged and identified particle yields, the mean transverse momenta and the integrated flow harmonics in 2.76 A TeV Pb+Pb collisions using the Bayesian statistics[24], which also nicely described various flow data at the LHC [26]. In practice, the transition from the hydrodynamic fluid to the emitted hadrons on the freeze-out surface is realized by a Monte-Carlo event generator iss based on the Cooper-Fryer formula[27]:

[TABLE]

where $f(x,p)$ is the distribution function of particles, $g$ is the degeneracy factor, and $d^{3}\sigma_{\mu}$ is the volume element on the freeze-out hypersurface.

For the following PCA analysis, as well as for the traditional flow analysis in comparison, we run the event-by-event VISH2+1 simulations with 12000 fluctuating initial conditions generated from TRENTo for 2.76 A TeV Pb-Pb collisions at 0%-10%,10%-20%, 20%-30%, 30%-40%, 40%-50% and 50%-60% centrality bins. The default iss sampling for each VISH2+1 simulation is 1000 events, which corresponds to the main results presented in Sec. III. In the appendix of this paper, we also investigate the ability of PCA to distinguish signal and noise. We thus implement 25, 100 and 500 iss samplings for each VISH2+1 simulation for such investigation. Note that the default 1000 iss sampling used in this paper has already dramatically suppressed the statistical fluctuations from noises for the final hadron distributions.

With the final particle distributions obtained from hydrodynamic simulations, various flow observables can be calculated based on the traditional flow harmonics defined by the Fourier decomposition in Eq.(1). In Sec.III, the traditional flow results will be served as the comparison to the PCA results.

2.2 Principal Component Analysis (PCA)

Principal Component Analysis (PCA) is a statistical method to analyze complicated data, which aims to transform a set of correlated variables into various independent variables via orthogonal transformations. These obtained main eigenvectors, associated with large or unnegligible singular values, are also called the principal components, which reveal the most representative characteristics of the data. In practice, PCA implements the Singular Value Decomposition (SVD) to a real matrix, which obtains a diagonal matrix with the diagonal elements arranged in a descending order. Therefore, one needs to first construct a related matrix before the following PCA and SVD analysis. Since this paper focuses on investigating the integrated flow with PCA, such final state matrix $\mathbf{M_{f}}$ is constructed from the angular distribution of all charged hadrons $dN/d\varphi$ $(|y|<1.0)$ (obtained from Eq.(2)) of $N=2000$ independent events in each centrality bin, using VISH2+1 simulations with TRENTo initial conditions. In more details, we divide the azimuthal angle $[-\pi,\pi]$ into $m=50$ bins and count the number of particles in each bin. For the $j_{th}$ bin in event ($i$), the number of particles is denoted as $dN/d\phi^{(i)}_{j}$, which is also the $i_{th}$ row and $j_{th}$ column of the matrix $\mathbf{M_{f}}$ 111In practice, we normalize the event vector in $\mathbf{M_{f}}$ to get rid of the multiplicity fluctuations..

Then, we apply SVD to the final state matrix $\mathbf{M_{f}}$ with the size $N\times m$ (Here, $N=2000$ and $m=50$), which gives

[TABLE]

where $\mathbf{{X}}$ and $\mathbf{{Z}}$ are two orthogonal matrices with the size of $N\times N$ and $m\times m$, respectively. $\mathbf{{\Sigma}}$ is a diagonal matrix with diagonal elements (singular values) arranged in the descending order $\sigma_{1}>\sigma_{2}>\sigma_{3}\ \cdots>0$.

With such matrix multiplication, the $i_{th}$ row of matrix $\mathbf{M_{f}}$, denoted as $dN/d\varphi^{(i)}$, can be expressed by the linear combination of the eigenvectors $z_{j}$ (the $j_{th}$ row of matrix $\mathbf{Z}$) with $j=1,2,...,m$:

[TABLE]

where $(i)=1,2,...,N$, represents the index of the event, $m$ is the number of angular bins of the inputting events. $\tilde{v}_{j}^{(i)}$ is the corresponding coefficient of ${z}_{j}$ for the $i_{th}$ event. In the spirit of PCA, one only focuses on the most important components, so there is a cut at the indices ${k}$ in the last approximation of Eq.(4). In Sec. III, we will show that ${k}=12$ is a proper truncation for the integrated flow analysis, and the shape of the bases or eigenvectors ${z}_{j}\ (j=1,...,{k})$ is similar to but not identical with the Fourier transformation bases $\cos(n\varphi)$ and $\sin(n\varphi)$ ($n=1,...,6$) used in the traditional method. Correspondingly, $\tilde{v}_{j}^{(i)}\ (j=1,...,{k})$ is identified as the real or imaginary part of the flow harmonics for event ($i$), and the singular values ${\sigma}_{j}$ are associated with the corresponding event averaged flow harmonics at different orders. For more details, please also refer to Sec. III.

3 Results

In this section, we implement PCA to analyze the single particle distributions $dN/d\varphi$ from hydrodynamics simulations in Pb+Pb collisions at $\sqrt{s_{NN}}=$ 2.76 A TeV. Firstly, we focus on the singular values, eigenvectors as well as the associated coefficients of PCA and explore if such unsupervised learning could discover flow with its own bases.

In practice, we run 2000 event-by-event VISH2+1 hydrodynamic simulations with TRENTo initial conditions to generate the $dN/d\varphi$ distributions for 10%-20% Pb+Pb collisions at $\sqrt{s_{NN}}=$ 2.76 A TeV. With these $dN/d\varphi$ distributions, we construct the final state matrix $\mathbf{M_{f}}$ and then implement SVD to $\mathbf{M_{f}}$ as described in Sec. II. Fig. 1 shows these obtained first 12 eigenvectors ${z}_{j}\ (j=1,2,...,12)$ and the first 20 singular values ${\sigma}_{j}\ (j=1,2,...,20)$ of PCA, arranged by the descending order of magnitude 222Each eigenvector is automatically normalized with $||z_{j}||_{2}^{2}=\sum_{i=1}^{m}(z_{j})_{i}^{2}=1$ ($m=50$), due to the orthogonality of the eigenvector matrix $\mathbf{Z}$.. As introduced in Sec. II, these eigenvectors contain the most representative information on correlations among final particles. Fig. 1 shows that the $1{st}$ and $2{nd}$ eigenvectors from PCA are similar to the Fourier decomposition bases $\mathrm{sin}(2\varphi)$ and $\mathrm{cos}(2\varphi)$, and the $3{rd}$ and $4{th}$ components are similar to $\mathrm{sin}(3\varphi)$ and $\mathrm{cos}(3\varphi)$, etc. Meanwhile, Fig. 1 (b) shows that singular values ${\sigma}_{j}\ (j=1,2,...,12)$ are arranged in pairs. These results indicate that each pair of the singular values may associate with the real and imaginary parts of the event averaged flow vectors at different orders. Therefore, we define the event averaged flow harmonics of PCA with these paired singular values, as outlined in the the second column of Table 1. The values of these PCA flow at different order are compared with the traditional flow harmonics from the Fourier expansion in Table 1, which are close, but not exactly the same values for $n\leq 6$.

As explained in Sec II, one could also read the event-by-event flow harmonics from the results of PCA. In more details, such PCA flow harmonics for event ($i$) is associated with these coefficients $\tilde{v}_{j}^{(i)},j=1...k$ in Eq. (4). Therefore, we define the event-by-event flow harmonics $v_{n}^{\prime}$ with magnitudes projected onto PCA bases, similar to the event averaged ones defined in Table 1. For example, $v_{2}^{\prime}=\sqrt{\frac{m}{2}}\sqrt{\tilde{v}_{1}^{2}+\tilde{v}_{2}^{2}}$ and $v_{3}^{\prime}=\sqrt{\frac{m}{2}}\sqrt{\tilde{v}_{3}^{2}+\tilde{v}_{4}^{2}}$ ($m=50$), etc. Fig. 2 compares $v_{n}^{\prime}$ from PCA and $v_{n}$ from the traditional Fourier expansion at different orders. For the event-by-event elliptic flow $v_{2}$ and $v_{2}^{\prime}$ and triangular flow $v_{3}$ and $v_{3}^{\prime}$, the definitions from PCA and that from Fourier expansion are highly agree with each other, which mostly fall on the diagonal lines. For these higher order flow harmonics with $n\geq 4$, these PCA results are largely deviated from the traditional Fourier ones. We also noticed that the first two PCA eigenvector ${z}_{1}$ and ${z}_{2}$ for $v_{2}^{\prime}$ are similar to but not identical with the Fourier bases $\mathrm{sin}(2\varphi)$ and $\mathrm{cos}(2\varphi)$ with $n=2$, which contain the contributions from $\mathrm{sin}(4\varphi)$ and $\mathrm{cos}(4\varphi)$. Similarly, the PCA eigenvectors ${z}_{3}$ and ${z}_{4}$ also contain the contributions from other Fourier bases. Such mode mixing in the PCA eigenvectors leads to the large deviations between $v_{4}$ and $v_{4}^{\prime}$, as well as between $v_{5}$ and $v_{5}^{\prime}$, etc.

To evaluate the correlations between different PCA flow harmonics $v_{m}^{\prime}$ and $v_{n}^{\prime}$, we calculate the symmetric cumulants as once defined for traditional flow harmonics [28, 29, 30]:

[TABLE]

Correspondingly, the traditional symmetric cumulants $SC^{v}{(m,n)}$ just replace $v^{\prime}_{m}$ and $v^{\prime}_{n}$ with $v_{m}$ and $v_{n}$ from the Fourier expansion.

FIG. 3 compares the symmetric cumulants $SC^{v}{{}^{\prime}(m,n)}$ from PCA and $SC^{v}{{}^{\prime}(m,n)}$ from Fourier expansion, for the event-by-event VISH2+1 simulations in 2.76 A TeV Pb+Pb collisions at various centrality bins. One finds that, except for $SC^{v}(2,3)$, almost all PCA symmetric cumulants $SC^{v}{{}^{\prime}(m,n)}$ reduce significantly compared to the traditional ones. Although $v^{\prime}_{4}$ from PCA largely deviated from the traditional $v_{4}$ from the Fourier expansion, the obtained $SC^{v}{{}^{\prime}(2,4)}$ shows a significant suppression, which contradicts to the long believed idea that the nonlinear hydrodynamics evolution strongly couples $v_{2}^{2}$ to $v_{4}$, leading to an obvious positive correlations between $v_{2}$ and $v_{4}$ obtained from Fourier expansion. Similarly, the non-linear mode coupling between $v^{\prime}_{2}$ and $v^{\prime}_{5}$, $v^{\prime}_{3}$ and $v^{\prime}_{5}$ and $v^{\prime}_{3}$ and $v^{\prime}_{4}$ for these PCA defined flow harmonics also decrease, which results in the reduced symmetric cumulants $SC^{v}{{}^{\prime}(2,5)}$, $SC^{v}{{}^{\prime}(3,5)}$ and $SC^{v}{{}^{\prime}(3,4)}$ correspondingly.

To evaluate the correlations between the initial and final state fluctuations, we use the Pearson coefficients $r(v^{\prime}_{n},\varepsilon_{m})$ and $r(v_{n},\varepsilon_{m})$ to characterize the linearity between the PCA flow harmonics $v^{\prime}_{n}$ and the initial eccentricities $\varepsilon_{m}$, as defined as the following:

[TABLE]

Here, $\varepsilon_{m}$ is the traditional eccentricities defined by Eq.(7). In Appendix A, we will demonstrate that, with a properly chosen smoothing procedure, the event-by-event eccentricities $\varepsilon^{\prime}_{m}$ from PCA is highly similar to $\varepsilon_{m}$ from the traditional method. We thus use $\varepsilon_{m}$ in the Pearson coefficient definition $r(v^{\prime}_{n},\varepsilon_{m})$ for PCA. Meanwhile, we can also calculate the Pearson coefficient $r(v_{n},\varepsilon_{m})$ for the traditional flow with Fourier expansion, which just replaces the flow harmonics $v^{\prime}_{n}$ in Eq. (6) by $v_{n}$. According to the definition, the Pearson coefficient falls in the range $[-1,1]$, with $r>0$ implying a positive correlation, and $r<0$ implying a negative correlation.

Fig. 4 plots the Pearson coefficients $r(v^{\prime}_{n},\varepsilon_{m})$ from PCA and $r(v_{n},\varepsilon_{m})$ from the Fourier expansion, for VISH2+1 simulated Pb+Pb collisions at various centralities. With these Pearson coefficients, we focus on evaluating if the PCA defined flow harmonics reduce or increase the correlations with the corresponding initial eccentricities. As shown in Fig. 3, the event-by-event flow harmonics $v^{\prime}_{2}$ or $v^{\prime}_{3}$ from PCA are approximately equal to the Fourier ones $v_{2}$ or $v_{3}$. As a result, these Pearson coefficients involved with these two flow harmonics $r(v^{\prime}_{2},\varepsilon_{m})$ and $r(v^{\prime}_{3},\varepsilon_{m})$ are almost overlap with the Fourier ones $r(v_{2},\varepsilon_{m})$ and $r(v_{3},\varepsilon_{m})$ as shown by these upper panels in the first two rows. Meanwhile, these diagonal Pearson coefficients $r(v^{\prime}_{2},\varepsilon_{2})$ or $r(v_{2},\varepsilon_{2})$ and $r(v^{\prime}_{3},\varepsilon_{3})$ or $r(v_{3},\varepsilon_{3})$ are much larger than other ones, which confirms the early conclusion that the elliptic flow and triangular flow are mainly influenced by the initial eccentricity $\varepsilon_{2}$ and $\varepsilon_{3}$ with the approximate linear relationship $v_{2}\thicksim\varepsilon_{2}$ ($v^{\prime}_{2}\thicksim\varepsilon_{2}$) and $v_{3}\thicksim\varepsilon_{3}$ ($v^{\prime}_{3}\thicksim\varepsilon_{3}$) [31, 32].

Although $v^{\prime}_{4}$ from PCA is largely deviated from the traditional $v_{4}$ in Fig. 3, such PCA definition largely enhances correlations between $\varepsilon_{4}$, and also largely reduces the correlations between $\varepsilon_{2}$. For example, at 20-30% centrality, the Pearson coefficients $r(v_{4},\varepsilon_{4})$ is only 70% of the $r(v_{4}^{\prime},\varepsilon_{4})$, while $r(v_{4},\varepsilon_{2})$ is 200% larger than $r(v^{\prime}_{4},\varepsilon_{2})$. Traditionally, it is generally believed that $v_{4}$ is largely influenced by $\varepsilon_{2}^{2}$ through the non-linear evolution of hydrodynamics. Our PCA analysis showed that such mode mixing could be deduced through a redefined PCA bases. Meanwhile, such PCA defined bases also significantly reduce the mode mixing for other higher order flow harmonics such as between $v^{\prime}_{5}$ and $\varepsilon_{2}$, $v^{\prime}_{5}$ and $\varepsilon_{3}$, etc.

4 Conclusions

In this paper, we implemented Principal Components Analysis (PCA) to study the single particle distributions of thousands of events generated from VISH2+1 hydrodynamic simulations. Compared with the early PCA investigations on flow that imposed the Fourier transformation in the input data [13, 15, 16, 17, 18], we focused on analyzing the raw data of hydrodynamics and exploring if a machine could directly discover flow from the huge amount of data without explicit instructions from human-beings. We found that the PCA eigenvectors are similar to but not identical with the traditional Fourier basis. Correspondingly, the obtained flow harmonics $v_{n}^{\prime}$ from PCA are also similar to the traditional $v_{n}$ for $n=2$ and 3, but largely deviate from the Fourier ones for $n\geq 4$. With these PCA flow harmonics, we found that, except for $SC^{v}{{}^{\prime}(2,3)}$, almost all other symmetric cumulants $SC^{v}{{}^{\prime}(m,n)}$ from PCA decrease significantly compared to the traditional $SC^{v}{(m,n)}$. Meanwhile, some certain Pearson coefficients $r(v^{\prime}_{n},\varepsilon_{m})$ that evaluate the linearity between the PCA flow harmonics and the initial eccentricities are obviously enhanced (especially for $n\geq 4$), together with an corresponding reduction of the off-diagonal elements.

These results indicate that PCA has the ability to discover flow with its own basis, which also reduce the related mode coupling effects, when compared with traditional flow analysis based on the Fourier expansion. We emphasis that these eigenvectors from PCA are modeled to be orthogonal and uncorrelated to each other. As a result, most of the symmetric cumulants $SC^{v}{{}^{\prime}(m,n)}$ from PCA that evaluate the correlations between different flow harmonics are naturally reduced compared with the traditional ones. Besides, the PCA flow harmonics $v^{\prime}_{n}$ presents an enhanced linear relationship to the corresponding eccentricities $\varepsilon_{n}$, especially for $n=4$. These results seem contradictory to the long believed idea that hydrodynamics evolution are highly non-linear, which leads to strong mode-coupling between different flow harmonics. Our PCA investigation has shown that such mode coupling effects could be reduced with new-defined bases for the flow analysis. With such finding, the non-linearity of the relativistic hydrodynamic systems created in heavy ion collisions should be re-evaluated, which we would like to further explore it with such PCA method in the near future.

Acknowledgements

We would like to thank the fruitful discussions with J. Jia, R. Lacey, D. Teaney and M. Zhou . This work is supported by the NSFC and the MOST under grant Nos. 11675004, 11435001 and 2015CB856900. We also gratefully acknowledge the extensive computing resources provided by the Super-computing Center of Chinese Academy of Science (SCCAS), Tianhe-1A from the National Supercomputing Center in Tianjin, China and the High-performance Computing Platform of Peking University.

Appendix A PCA for initial profiles with smoothing procedure

In this appendix, we focus on analyzing the initial state fluctuations using the PCA method. Traditionally, the initial state fluctuations are evaluated by the eccentricity coefficients $\varepsilon_{n}$ at different order, which are defined as [31]:

[TABLE]

where $\Phi_{n}$ is the participant plan angle, $s(r,\varphi)$ is the initial entropy density and $\varphi$ is the azimuthal angle in the transverse plane [31].

For the PCA analysis, we first construct the initial state matrix $\mathbf{M_{i}}$, using the azimuthal angle distribution of the initial entropy $dS/d\varphi$ which is defined by

[TABLE]

obtained from 2000 event-by-event TRENTo initial conditions. A direct PCA analysis shows that more than 100 eigenvectors are needed to capture the rich structures of the initial state fluctuations. In contrast, 12 PCA eigenvectors are enough to describe the final state ones since the hydrodynamic evolution tends to smear out inhomogeneity of the evolving systems. In order to connect and compare these PCA singular values from the initial and final states, we implement a smoothing procedure for the initial profiles before the PCA analysis.

In more details, we apply a circular convolution with to the initial density profile $dS/d\varphi$, which is written as:

[TABLE]

Here, $K(\varphi^{\prime},\varphi)$ is the convolution kernel, which is taken a gaussian form $K(\varphi^{\prime},\varphi)=\frac{1}{\sqrt{2\pi}a}e^{-\frac{(\varphi^{\prime}-\varphi)^{2}}{2a^{2}}}$. Here, we fine tune the radius $a$ to ensure the the same decaying rate for the PCA singular values from the initial profiles and final profiles. The obtained optimized value for $a$ is 0.251 rad.

With such smoothing procedure, we reconstruct the initial state matrix $\mathbf{M_{i}}$ with 2000 event-by-event $(\frac{dS}{d\varphi})_{smooth}$ distributions from TRENTo for each selected centralities. As the case for the flow analysis in Sec. II B and Sec.III, the implementation of SVD and PCA to the initial state matrix, $\mathbf{M_{i}}=\mathbf{\hat{Y}\hat{\Sigma}\hat{Z}}=\mathbf{\hat{E}}\mathbf{\hat{Z}}$, gives the singular value $\hat{\sigma_{j}}$, eigenvectors $\hat{z_{j}}$ and the corresponding eccentricity coefficients $\hat{\varepsilon}_{j}^{(i)}$, $(j=1,...\hat{k})$ such that

[TABLE]

We find that the PCA eigenvectors $\hat{z_{j}}$ of the initial states are highly similar to traditional Fourier bases $\cos(2\varphi)$, $\sin(2\varphi)$, $\cos(3\varphi)$, $\sin(3\varphi)$, etc. Meanwhile, we could associate the singular value $\hat{\sigma}_{j}$ to the event averaged initial eccentricities of PCA, $\bar{\varepsilon}^{\prime}_{n}$, at different orders and connect the coefficients $\hat{\varepsilon}_{j}^{(i)},(j=1,...\hat{k})$ to the real or imaginary part of the PCA event-by-event initial eccentricities $\varepsilon^{\prime}_{n}(n=1,...\hat{k}/2)$ as the case for flow 333Here, $\bar{\varepsilon}^{\prime}_{2}$=$\sqrt{\frac{m}{2}}\sqrt{\hat{\sigma}_{3}^{2}+\hat{\sigma}_{4}^{2}}$, $\bar{\varepsilon}^{\prime}_{3}$=$\sqrt{\frac{m}{2}}\sqrt{\hat{\sigma}_{5}^{2}+\hat{\sigma}_{6}^{2}}$, $\bar{\varepsilon}^{\prime}_{4}$=$\sqrt{\frac{m}{2}}\sqrt{\hat{\sigma}_{7}^{2}+\hat{\sigma}_{8}^{2}}$, etc. with $m=50$ the number of bins. For event-by-event definition of ${\varepsilon}^{\prime}_{n}(n=1,...\hat{k}/2)$, we could simply replace $\hat{\sigma}_{j}$ with $\hat{\varepsilon}_{j}$ $(j=1,...\hat{k})$, correspondingly..

Fig. 5 compares the event-by-event eccentricites $\varepsilon^{\prime}_{n}$ from PCA and $\varepsilon_{n}$ from the traditional definition (A.1) for the TRENTo initial conditions at 10-20% centrality. It shows, with a properly chosen smoothing procedure of the initial conditions, $\varepsilon^{\prime}_{n}$ and $\varepsilon_{n}$ agree with each other well till ${n=6}$ . Meanwhile, the event averaged eccentricites $\overline{\varepsilon^{\prime}_{n}}$ from PCA and $\overline{\varepsilon_{n}}$ from (7) also fit each other very well, which is much better than ones for flow shown in Table 1 and Fig. 2. Therefore, for the investigation of initial state and final state correlations, we only use the traditional $\varepsilon_{m}$ to define the Pearson coefficient $r(v^{\prime}_{n},\varepsilon_{m})$, and $r(v_{n},\varepsilon_{m})$ for both PCA and traditional flow in Sec.III.

Appendix B Signal and noise distinguishment from PCA

In the event-by-event VISH2+1 simulations, both initial state fluctuations and statistical fluctuations from the iss particle sampling during the Cooper-Fryer freeze-out influence the emissions and distributions of final particles. It is generally believed that the hydrodynamic evolution translate the initial state fluctuations into final state correlations, which directly relate to flow signals. Meanwhile, the statistical fluctuations during Cooper-Fryer freeze-out with a finite number of particle emission introduce statistical noise for the flow definition in each event. As a result, flow harmonics from traditional Fourier expansion are generally analyzed with an event average of millions of events. For the event-by-event flow analysis, one implements the standard Bayesian unfolding procedure to suppress effects from the finite multiplicites and non-flow [33].

In this appendix, we further explore the ability of PCA to distinguish the signal and noise. With such purpose, we implement 25, 100 and 500 iss samplings for each VISH2+1 simulation to generate the $dN/d\varphi$ distributions of final particles and the related final state matrixes $\mathbf{M_{f}}$ with different weighted signal and noise. Then, we implement PCA to analyze these matrixes. As shown in Fig. 6, the distribution of the PCA singular values is changed with the number of iss samplings. For these systems with large statistical fluctuations, for example with 25 iss samplings, the singular values $\sigma_{j}$ at large $j$ tend to have a long and high tail. For these systems with reduced statistical fluctuations with more iss samplings, the height of the tail is largely decreased. Meanwhile, we noticed that these eigenvectors with an index $j$ smaller than a certain “magic number”(12 in this case) is signal-like which has a basis similar to the Fourier one, while these eigenvectors with larger $j$ behave so randomly and chaotically, that we associate these eigenvectors with the noise patterns of the systems.Besides, we check the height of these PCA tails and found the ratios among these heights for different iss samplings approximately satisfy $\frac{1}{\sqrt{25}}:\frac{1}{\sqrt{100}}:\frac{1}{\sqrt{500}}$, such relation is known as the Law of Large Numbers for statistical noise. With more number of samplings, the height of the tail would further decrease. In the main part of this paper, we thus set the iss samplings for each VISH2+1 simulation to 1000, which largely suppresses the noise effects from the statistical fluctuations and makes PCA analysis focus on studying flow signal itself.