A Fourier-cumulant analysis for multiharmonic flow fluctuation: by employing a multidimensional generating function approach.

The Fourier analysis of the final particle distribution followed by cumulant study of the Fourier coefficient event-by-event fluctuation is one of the main approaches for testing the collective evolution in the heavy-ion collision. Using a multidimensional generating function, we propose a method to extract any possible cumulant of multiharmonic flow fluctuations and classify them in terms of the order of cumulants and harmonics involved in them. In particular, we show that there are 33 distinct cumulants with order 2, 3, 4, 5 and harmonics 2, 3, 4, 5. We compute the normalized version of these cumulants from hydrodynamic simulation for Pb-Pb collisions based on T R ENTo+VISH2+1+UrQMD. We compare the simulation with those normalized cumulants that the LHC has measured and predict the unmeasured ones. Comparing the initial and final state fluctuation normalized cumulants, we compute the linear and nonlinear hydrodynamic response couplings. We finally introduce the genuine three-particle correlation function containing information of all third-order cumulants.

Introduction

It is widely accepted in the heavy-ion community that QCD matter in the deconfined phase, quark-gluon plasma, is produced in ultrarelativistic heavy-ion collisions. The collective models equipped with QCD as an underlying physics can successfully explain the majority of experimental observations. However, measuring the accurate transport coefficients such as shear and bulk viscosity over entropy density and finding QCD critical point are still one of the main challenges in the heavy-ion community [1-3]. Moreover, observing similar collectivity signals in small system collisions such as p–p or p–Pb [4-10] raise debates about relation between the observations and the collective processes. These indicate that the ultimate goal has not been achieved, and we need to deepen our understanding of heavy-ion collision physics. To this end, introducing new observables to probe various aspects of the full image is crucial. The present manuscript tries to introduce a systematic procedure to generate observables related to the anisotropic flow and classify them. Some of these observables are studied before, and some others are introduced for the first time.

The anisotropic particle emission in the azimuthal direction [11-13] is one of the most important evidences of collectivity in ultrarelativistic heavy-ion collisions. The final particle momentum distribution in the azimuthal direction can be expanded in a Fourier series,with Fourier coefficients , called flow harmonics. The imprints of different stages of a heavy-ion collision, pre-equilibrium, initial state, collective evolution, and freezeout are cumulatively encoded in the flow harmonics. On account of the quantum mechanical nature of partons inside the nucleons, the initial geometry has a complicated structure and more importantly fluctuates from one event to the other (event-by-event fluctuation) leading to nonvanishing and fluctuating and for any n [4, 14–23]. Practically, we measure the statistical properties of the fluctuation in the experiment. Several attempts focus on the moments [24-33] or on cumulants of fluctuations to achieve information from flow harmonic fluctuations. Single harmonic cumulants [34, 35], symmetric cumulants [36-39], generalized (or higher-order) symmetric cumulants [40, 41], asymmetric cumulant [24] are some examples of these attempts.

It is noteworthy that each cumulant of a probability density function (p.d.f.) contains independent information about the underlying fluctuation. Therefore, a complete study of all possible cumulants up to a given order of cumulant expansion is crucial to gain a comprehensive insight into the fluctuations. In the present work, we start with the standard definition of cumulants based on generating functions for an arbitrary number of variables. We extract and classify all possible cumulants depending on their order and the involved number of harmonics. Considering the space limitation, we focus on cumulants of harmonics up to fifth order, which contains 33 distinct cumulants. Among them, there are already known cumulants together with new cumulants containing symmetry plane correlations. To ease the future application, we have prepared a Mathematica [42] package accessible as an ancillary file for the present manuscript or in the GitHub repository [43]. The package returns the cumulants in three different forms: in terms of and moment, final particle azimuthal angle correlations, and Q-vectors (see Eq. (29)) notations. For cases that the Mathematica software is not available, a list of few first multiharmonic cumulants in terms of and symbolic moments are tabulated in the Appendix A.

Regarding multiharmonic cumulants, other studies have been done previously [44, 45] with some similarities and differences compared to the present work. In Ref. [44], the underlying p.d.f. of the fluctuations are classified into three categories, flow-amplitude (only fluctuations), event-plane-correlation (only fluctuations) and mixed-correlations (both and fluctuations). We find that, instead of classifying the p.d.f. into three types, considering one general p.d.f. is more suitable way to extract all possible cumulants up to a given degree. The author believes that the cumulants extracted in the present paper are achievable with the method introduced earlier in Ref. [45]. In this reference, however, the explicit calculations are limited to the simple few known observables, and the dependence of the cumulant expressions to the harmonics is not presented. Our explicit calculations show that the cumulants’ form is related to the involved harmonics in the cumulant (see also Refs. [44, 46]). For instance, we see that the first cumulant involving simultaneous harmonics appears in the third order, and the next cumulants appear in the fifth order in three different forms. For the case involving simultaneous harmonics , however, cumulants starts from fifth order with two different forms (see Table 1). Besides, our study is equipped with Monte Carlo hydrodynamic simulation as well (TENTo [47]+VISH2+1 [48, 49]+UrQMD [50, 51] with parameter calibrated by a Global Bayesian analysis [52]) to see which unmeasured cumulants would have larger signals in future measurements. Comparing this simulation with future measurements can also validate the heavy-ion collision parameters tuned with Global Bayesian analysis. We compare the simulation with normalized cumulants that are measured at the LHC, which reveals a rather good agreement.

Table 1

List of all cumulants with order and harmonics

	Cumulant	Order	Cumulant expression
1	c_2{2}	2	⟨v22⟩
2	c_3{2}	2	⟨v32⟩
3	c_4{2}	2	⟨v42⟩
4	c_5{2}	2	⟨v52⟩
5	c2,4{4}{2,1}	3	⟨v22v_4cos4ψ_2-ψ_4⟩
6	c2,3,5{-3,5}{1,1,1}	3	⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
7	c_2{4}	4	⟨v24⟩-2⟨v22⟩^2
8	c_3{4}	4	⟨v34⟩-2⟨v32⟩^2
9	c_4{4}	4	⟨v44⟩-2⟨v42⟩^2
10	c_5{4}	4	⟨v54⟩-2⟨v52⟩^2
11	c2,3{0}{2,2}	4	⟨v22v32⟩-⟨v22⟩⟨v32⟩
12	c2,4{0}{2,2}	4	⟨v22v42⟩-⟨v22⟩⟨v42⟩
13	c2,5{0}{2,2}	4	⟨v22v52⟩-⟨v22⟩⟨v52⟩
14	c3,4{0}{2,2}	4	⟨v32v42⟩-⟨v32⟩⟨v42⟩
15	c3,5{0}{2,2}	4	⟨v32v52⟩-⟨v32⟩⟨v52⟩
16	c4,5{0}{2,2}	4	⟨v42v52⟩-⟨v42⟩⟨v52⟩
17	c2,3,4{6,-4}{1,2,1}	4	⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩
18	c3,4,5{8,-5}{1,2,1}	4	⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩
19	c2,3,4,5{3,4,-5}{1,1,1,1}	4	⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩
20	c2,3{6}{3,2}	5	⟨v23v32cos6ψ_2-ψ_3⟩
21	c2,4{4}{2,3}	5	⟨v22v43cos4ψ_2-ψ_4⟩-2⟨v42⟩⟨v22v_4cos4ψ_2-ψ_4⟩
22	c2,4{4}{4,1}	5	⟨v24v_4cos4ψ_2-ψ_4⟩-3⟨v22⟩⟨v22v_4cos4ψ_2-ψ_4⟩
23	c2,3,4{-6,8}{1,2,2}	5	⟨v42v32v_2cos2ψ_2+3ψ_3-4ψ_4⟩
24	c2,3,4{0,4}{2,2,1}	5	⟨v32v22v_4cos4ψ_2-ψ_4⟩-⟨v32⟩⟨v22v_4cos4ψ_2-ψ_4⟩
25	c2,3,5{-3,5}{1,1,3}	5	⟨v53v_2v_3cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v52⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
26	c2,3,5{-3,5}{1,3,1}	5	⟨v33v_2v_5cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v32⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
27	c2,3,5{-3,5}{3,1,1}	5	⟨v23v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v22⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
28	c2,4,5{-8,10}{1,2,2}	5	⟨v52v42v_2cos2ψ_2+4ψ_4-5ψ_5⟩
29	c2,4,5{4,0}{2,1,2}	5	⟨v52v22v_4cos4ψ_2-ψ_4⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v52⟩
30	c3,4,5{-4,10}{2,1,2}	5	⟨v52v32v_4cos6ψ_3+4ψ_4-10ψ_5⟩
31	c3,4,5{4,5}{3,1,1}	5	⟨v33v_4v_5cos9ψ_3-4ψ_4-5ψ_5⟩
32	c2,3,4,5{-3,0,5}{1,1,2,1}	5	⟨v42v_2v_5v_3cos2ψ_2+3ψ_3-5ψ_5⟩-⟨v42⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
33	c2,3,4,5{3,-4,5}{2,1,1,1}	5	⟨v22v_3v_4v_5cos4ψ_2-3ψ_3+4ψ_4-5ψ_5⟩

We study the initial state fluctuation cumulants to see how the final state inherits the fluctuations from the initial state. The initial energy density anisotropy translation to the final state momentum anisotropy is formulated by a set of equations so-called hydrodynamic response [53]. By comparing the initial and final state fluctuations, we obtain the hydrodynamic response coefficients. Moreover, our systematic cumulant study guide us to extend the two-particle correlation function (2PC) notion [5, 54–56] to multiparticle correlation functions, qPC (see the definition in Sect. 6). We specifically study the flow-induced 3PC for harmonics 2 to 5 via our hydrodynamic simulation. This quantity contains information on all third-order cumulants.

This paper is structured as follows: The Sect. 2 is dedicated to introduce generating function method for multiharmonic cumulants. In Sect. 3, the functions available in the Mathematica package are introduced. A realistic Monte Carlo study of the first Fourier-cumulant expansion terms is presented in Sect. 4. The linear and nonlinear hydrodynamic response couplings are extracted by comparing the initial state, and final state fluctuation normalized cumulants in Sect. 5. Finally, we study genuine three-particle correlation functions in Sect. 6. In the Appendix A, a list of few first cumulants is tabulated. In Appendix B, a technical study about the statistical fluctuation of multiharmonic correlations is presented. Detail of some derivations in the text is presented in the Appendices C and D.

Multiharmonic generating function

Despite the simplicity of Eq. (1), experimental measurement of flow harmonics is a challenging task. The number of final produced particles per event in a given collision is not enough to extract statistically accurate flow harmonics values. Many events average resolves the problem of low statistics while leads to the convolution of event-by-event fluctuation of the initial state into the measurement. Strictly speaking, the measurable quantities are the moments or cumulants of the following gigantic probability density function (p.d.f.)where we have used the Cartesian coordinate notation for the flow harmonics, and . Noting the fact that (Eq. (1)) is assigned to each initial state, the p.d.f. (2) is tightly related to the event-by-event fluctuation of the initial state, which is, accordingly, related to the underlying physics of nucleons/nuclei. Moreover, the map which connects the above p.d.f. to the initial state fluctuation is governed by the collective evolution process. Indeed, there are stochastic processes in the collective evolution that can be convoluted into the above p.d.f. as well. Besides the experiment’s statistical uncertainty issues, the reaction plane angle is not experimentally measurable conveniently, which leads to loss of one more degree of freedom in p.d.f. (2). In the polar coordinate the p.d.f. (2) depends on k flow amplitudes s and k symmetry plane angles (or event-planes) s. After imposing the randomness of the reaction plane angle, from k variable , only independent variables remainHere, is chosen conventionally to compute the angle differences. One could choose any other combination of two symmetry plane angle differences. The final cumulants are independent of this choice.

In the present section, we introduce the method of generating function to extract all possible cumulants of the p.d.f. (3). To this end, we start with a standard definition of generating function of a generic multivariate p.d.f. (see for instance [57]),where is a m-dimensional random vector. The cumulants of the distribution function (shown by ) are the Taylor expansion coefficients of the cumulant generating function ,The above generating function returns the cumulants of p.d.f. (2), while in heavy-ion physics, we are interested in cumulants of p.d.f. (3). To clarify the difference between cumulants of p.d.f. (2) compared to p.d.f. (3), let us start with the most simple example. Focusing on one-harmonics , the p.d.f. of the flow fluctuation in the Cartesian coordinate has the form , and its cumulants can be obtained via the definition presented in Eq. (5) for . For instance,are the width of the distribution in the and the skewness in the direction, respectively. In the experiment, however, we lose one of the degrees of freedom in due to the randomness of the reaction plane angle. After imposing event-by-event random reaction plane angle to , we obtain a rotationally symmetric p.d.f. . The information of the latter p.d.f. is encoded in a radial distribution . Traditionally, we assume the cumulants of are [35] where the two first orders of them are given byThese quantities, however, are specific combinations of two-dimensional cumulants of that survived after imposing the randomness of the reaction plane angle. For instance, it turns out [58, 59],Here, we have assumed that and for simplicity. This condition is satisfied when the flow harmonics are sourced purely from fluctuations, namely ellipticity in central Pb–Pb collisions, or triangularity in round nuclei collisions. For nonvanishing , similar relations with more terms in the right-hand side of Eq. (8) are obtained. In any case, the cumulants are the maximum information we can achieve from the cumulants of the original distribution .

Before proceeding, let us briefly review the generating functional method to find [60]. The generating function in Eq. (4) corresponds to the rotationally symmetric reads aswhere in the above is the Bessel function of the first kind and the averaging is performed with respect to . One can expand the logarithm of the generating function in terms of k,withto find the cumulants .

We extend the above procedure into the cases with an arbitrary number of harmonics. To simplify the notation, we rewrite the p.d.f. (3) as whereandConsidering the rotational symmetry imposed by the randomness of the reaction plane angle, one can first average the characteristic function over the azimuthal angle to find a symmetric characteristic function. Then the nonvanishing combination of cumulants of can be obtained directly from the symmetric characteristic function,whereandFor the derivation refer to the Appendix D. Having found , the cumulants of flow fluctuations can be obtained from the following expansion:where the coefficientsare the cumulants ( are defined in Eq. (13)). In the right-hand side of the above equation, we have represented the cumulant with its highest rank moment together with a subscript c. Here, stand for the involving harmonics, are the power of flow amplitudes and are the coefficients of the symmetry plane angle differences. The order of the cumulant is given byBy definition, s are all distinguished while, without loss of generality, we impose a strictly ascending order convention to them, .1 The coefficient is a numerical factor that does not depend on the moments. We fix the coefficient so that the numerical factor of the highest rank moment in cumulant turns to be equal to unity. As we will see in the next section, some combinations of harmonics , , and lead to vanishing cumulant. A trivial example is where the Eq. (18) reduces to . It is known that this cumulant is nonvanishing only for even m (see Eq. (7)).

Each distinguished cumulant contains a piece of independent information about the p.d.f. (3). To numerically compute these cumulants for a given p.d.f., we need to find the explicit form of them written in terms of the moments. Although the cumulants can be found analytically using Eq. (17), its computation is cumbersome, except for the case . By considering one flow harmonic (), the generating function reduces to that mentioned in Eq. (9). Keeping two flow harmonics, the flow fluctuation distribution contains three degrees of freedom , , and . In this case, the function can be written in terms of generalized Bessel function [61], and by expansion, one obtains the cumulants . However, extracting them is arithmetically more involved. For more general cases, the complexity increases, forcing us to choose a more practical and efficient way for the computation.

One package for all cumulants

Using the cumulant generating function in Eq. (17), one can find the cumulants written in terms of symbolic moments of variables and . These results can be immediately used in theoretical studies where and of every single event are accurately accessible. In the experiment as well as some simulations, however, the outcomes are particles azimuthal angles. As a result, one needs to invest an extra effort to translate the azimuthal angle of final particles into the averages containing and . For that, we use the multiparticle techniques to compute the flow harmonic cumulants. We rewrite the moments in terms of Q-vectors (see Eq. (29)) which can be calculated from final particle azimuthal angles at every single event and employed in computing the cumulant [36, 62]. This section presents a practical way to find cumulants and their statistical uncertainties in terms of symbolic moments, correlation of particle azimuthal angles, and Q-vectors.

Cumulants from generating function

As it is mentioned in the previous section, the analytical computation of a generic cumulant is cumbersome. For that reason, we do it symbolically in Mathematica. We encapsulated this computation into different functions (see below) and implemented them into a Mathematica package. One can load the package in a separate Mathematica notebook and recall the functions. This package is available as an ancillary file of this manuscript, or in the GitHub repository [43].

In the package, we first compute the Taylor expansion of the exponential function in Eq. (15) in terms up to order . Then we perform the integral in Eq. (15). After that, we replace the combinations in the expansion with a symbolic variable as a moment. Computing the logarithm of the result, we read the coefficients of the Taylor series for the variables and the Fourier series for the variables (see Eq. (17)). We extract the cumulants in terms of our symbolic moment variables up to a numerical factor by comparing the result with Eq. (17). Finally, we single out the highest rank moment in the result and fix the coefficient such that the numerical factor of this moment turns to unity.

We listed the available functions in a short manual at the header of the package file. One of these functions isthat returns the associated cumulants written in terms of moments of variables and . For instance,By using cMean instead of c, the angle brackets are replaced by Mean[]. For instance, cMean[{2},{},{2},v,] returns Mean[v[2]].

The above representation of cumulant is useful when flow harmonics are accurately accessible in a single event, typically in hydrodynamic simulations. If the final state is the particle azimuthal angles, one needs to employ the particle correlations,Here, is the azimuthal angle of the ith particle in an event, and M is the multiplicity of the event [34–36, 62]. A generic flow harmonic moment withcan be written in terms of particle correlations as . The functionreturns the cumulant written in terms of correlations. For instance,For external applications, it is easier if the output can be simply copied to the codes outside of Mathematica. This can be obtained viaUsing function () in Mathematica, one immediately convert the output of cCorr to an expression applicable in a () code. For instance,where,From the practical point of view, computing contains several nested loops depending on the value of k which are computationally expensive. A technique is introduced in Refs. [36, 62] where the moments are obtained in terms of Q-vectors,2Using this technique, only one loop over particles in an event is needed to compute the correlations.

For those who want to analyze anisotropic flow inside Mathematica, we implemented the recursive algorithm mentioned in Ref. [36] into the package to find in terms of Q-vectors. Substituting and symbolic averages with moments written in terms of Q-vectors, we find cumulants in terms Q-vectors (so-called Q-cumulants),where M and Q are the multiplicity and Q-vector symbols, respectively. For example,which is nothing butas it is already obtained in Refs. [36, 62]. Having defined M, Q[-2], and Q[2] as a list in Mathematica containing multiplicity, , and of many events, the function in Eq. (31) returns the numerical value of the cumulant . The function works accordingly for all other cumulants discuss in the present work.

The other function we advertise here isthat returns a table of all nonvanishing cumulants for flow harmonics and orders between and . It is worth mentioning that considering more flow harmonics and computing higher-order cumulants makes computations more time-consuming. For future studies and cases that Mathematica software is not available, we tabulate the cumulants of flow harmonics up to four harmonics combinations in the Appendix A. In Tables 3, 4, 5 and 6 in the Appendix A, we present the one- and two-harmonics cumulants up to eight, and three- and four-harmonics up to six orders, concerning the length of the final expressions. The interested reader can extract the cases which are not appeared in the tables by using the Mathematica package.

Table 3

One-harmonic cumulants up to order eight

{m}	c_n{m}
{2}	⟨vn2⟩
{4}	⟨vn4⟩-2⟨vn2⟩^2
{6}	⟨vn6⟩-9⟨vn4⟩⟨vn2⟩+12⟨vn2⟩^3
{8}	⟨vn8⟩-18⟨vn4⟩^2-16⟨vn6⟩⟨vn2⟩+144⟨vn4⟩⟨vn2⟩^2-144⟨vn2⟩^4

Table 4

Two-harmonic cumulants up to order eight

{m_1,m_2}	{α}	c2,3{α}{m_1,m_2}
{2,2}	{0}	⟨v22v32⟩-⟨v22⟩⟨v32⟩
{3,2}	{6}	⟨v23v32cos6ψ_2-ψ_3⟩
{2,4}	{0}	⟨v22v34⟩-4⟨v32⟩⟨v22v32⟩-⟨v22⟩⟨v34⟩+4⟨v22⟩⟨v32⟩^2
{4,2}	{0}	⟨v24v32⟩-⟨v24⟩⟨v32⟩-4⟨v22⟩⟨v22v32⟩+4⟨v22⟩^2⟨v32⟩
{3,4}	{6}	⟨v23v34cos6ψ_2-ψ_3⟩-3⟨v32⟩⟨v23v32cos6ψ_2-ψ_3⟩
{5,2}	{6}	⟨v25v32cos6ψ_2-ψ_3⟩-4⟨v22⟩⟨v23v32cos6ψ_2-ψ_3⟩
{2,6}	{0}	⟨v22v36⟩-9⟨v22v32⟩⟨v34⟩-9⟨v32⟩⟨v22v34⟩-⟨v22⟩⟨v36⟩+36⟨v32⟩^2⟨v22v32⟩+18⟨v22⟩⟨v32⟩⟨v34⟩-36⟨v22⟩⟨v32⟩^3
{4,4}	{0}	⟨v24v34⟩-8⟨v22v32⟩^2-4⟨v32⟩⟨v24v32⟩-⟨v24⟩⟨v34⟩-4⟨v22⟩⟨v22v34⟩+4⟨v24⟩⟨v32⟩^2+32⟨v22⟩⟨v32⟩⟨v22v32⟩+4⟨v22⟩^2⟨v34⟩-24⟨v22⟩^2⟨v32⟩^2
{6,2}	{0}	⟨v26v32⟩-⟨v26⟩⟨v32⟩-9⟨v24⟩⟨v22v32⟩-9⟨v22⟩⟨v24v32⟩+18⟨v22⟩⟨v24⟩⟨v32⟩+36⟨v22⟩^2⟨v22v32⟩-36⟨v22⟩^3⟨v32⟩

Table 5

Three-harmonic cumulants up to order six

{m_1,m_2,m_3}	{α_1,α_2}	c2,3,4{α_1,α_2}{m_1,m_2,m_3}
{1,2,1}	{6,-4}	⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩
{1,2,2}	{-6,8}	⟨v42v32v_2cos2ψ_2+3ψ_3-4ψ_4⟩
{2,2,1}	{0,4}	⟨v32v22v_4cos4ψ_2-ψ_4⟩-⟨v32⟩⟨v22v_4cos4ψ_2-ψ_4⟩
{1,2,3}	{6,-4}	⟨v32v43v_2cos2ψ_2-3ψ_3+2ψ_4⟩-2⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩⟨v42⟩
{1,4,1}	{6,-4}	⟨v34v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩-3⟨v32⟩⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩
{2,2,2}	{0,0}	⟨v22v32v42⟩-⟨v22v32⟩⟨v42⟩-⟨v32⟩⟨v22v42⟩-⟨v22⟩⟨v32v42⟩+2⟨v22⟩⟨v32⟩⟨v42⟩
{3,2,1}	{6,-4}	⟨v23v32v_4cos2ψ_2-3ψ_3+2ψ_4⟩-2⟨v22⟩⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩

Table 6

Four-harmonic cumulants up to order six

{m_1,m_2,m_3,m_4}	{α_1,α_2,α_3}	c2,3,4,5{α_1,α_2,α_3}{m_1,m_2,m_3,m_4}
{1,1,1,1}	{3,4,-5}	⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩
{1,1,2,1}	{-3,0,5}	⟨v42v_2v_5v_3cos2ψ_2+3ψ_3-5ψ_5⟩-⟨v42⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{2,1,1,1}	{3,-4,5}	⟨v22v_3v_4v_5cos4ψ_2-3ψ_3+4ψ_4-5ψ_5⟩
{3,1,1,1}	{-3,4,5}	⟨v23v_3v_4v_5cos6ψ_2+3ψ_3-4ψ_4-5ψ_5⟩-3⟨v22v_4cos4ψ_2-ψ_4⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{1,1,1,3}	{3,4,-5}	⟨v53v_2v_4v_3cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩-2⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩⟨v52⟩
{1,1,3,1}	{3,4,-5}	⟨v43v_2v_5v_3cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩-2⟨v42⟩⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩
{1,2,1,2}	{6,-4,0}	⟨v52v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩-⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩⟨v52⟩
{1,3,1,1}	{3,4,-5}	⟨v33v_2v_5v_4cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩-2⟨v32⟩⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩
{2,1,2,1}	{-3,8,-5}	⟨v42v22v_5v_3cos3ψ_3-8ψ_4+5ψ_5⟩-⟨v22⟩⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩
{3,1,1,1}	{3,4,-5}	⟨v23v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v22⟩⟨v_2v_3v_4v_5cos2ψ_2-3ψ_3-4ψ_4+5ψ_5⟩

Apart from cumulant which is obviously a special case of , several multiharmonic cumulants are known before, such as symmetric cumulants [36],generalized (or higher-order) symmetric cumulants for some combination of harmonics [40] (the actual form of the cumulant in terms of moments can be found in the Appendix A),3, and asymmetric cumulant [21, 22, 24],as special cases of what we present here. There are other cumulants, however, missed in other studies so far. For instance,is an example of a cumulant which have not been studied before.

Statistical uncertainty of cumulants

Having computed the flow harmonic cumulants, we discuss the statistical fluctuation of the cumulants. In case the flow harmonics are accurately accessible in a single event, one can follow a standard procedure (see for instance Ref. [63]) to find the covariance matrix of the moments (see Eq. (23)). The covariance matrix of two generic moments are simply given bywhere N is the number of events and the subscript i stands for a generic collective index . The variance of any function of moments, is obtained by [63],For instance, the variance of is given byIf the azimuthal angle of particles in the final state are available, the flow harmonic should be replaced by particle correlations . Although is not an accurate estimation for , in the ultimate many events average it approaches to an accurate estimation. By replacing with in Eq. (38), we obtain the variance of the function of correlations in terms of correlations. For instance, in this notation, we have , and accordingly, the example in Eq. (39) turns to the following form:where we have used the explicit form of covariance matrix in Eq. (38). One can immediately substitute the quantity in terms of Q-vectors for practical computations by using Refs. [36, 62].

The procedure explained above is implemented into the following functions in the package: The output of the above functions is . The input function (func) can be any function of correlations but the correlations should be always written in the form corr[] as the output of Eq. (24).4 For instance, to find by using Nsigma2, we should call the function as follows:The difference between functions in Eq. (41) is in their output presentations. The output of function Nsigma2 is shown in Eq. (42), the function Nsigma2Mean output is similar to cMean function, and Nsigma2Qvec returns the variance in terms of Q-vectors similar to cQvec function (see Eq. (31)). As an example, the statistical error of in terms of Q-vectors can be found below:Referring to Eqs. (40) and (42) (and to any other Nsigma2 outcomes), we find that the statistical uncertainty has no explicit multiplicity dependence. Let us consider two equal-size sets of events when the multiplicity of events in one set is smaller than in the other set. We expect that the statistical uncertainty of cumulants computed from the set of events with smaller multiplicity is larger than the statistical uncertainty obtained from the other set with higher multiplicity events. This apparent contradiction with intuition is due to the presence of correlations in the statistical uncertainty (for instance in Eq. (40)). These forms of correlations are the consequence of in Eq. (38) and have not appeared in previous studies in Refs. [36, 62]. In fact, the correlations depend on multiplicity due to the presence of autocorrelations remained in them. To be more specific, we focus on the first terms in the output in Eqs. (42),The above relation is a four-particle correlation with some remaining autocorrelations. In the technical Appendix B, a more direct approach is employed to extract the statistical uncertainties. This approach is implemented in the function Nsigma2P[func,M]. The result of this function for is as follows:As seen from the above equation, on the one hand, there is no correlation with the form , and on the other hand, there are explicit M dependences as we expect. One can explicitly show that the terms inside the first bracket in Eq. (45) are exactly equal to . It becomes more apparent when we replace all in Eq. (45) with Q-vectors using the results of Refs. [36, 62] where one finds that Eq. (45) in terms of Q-vectors is identical with that mentioned in Eq. (43). In general, the result of function (41b) is identical with that obtained from function Nsigma2P.

As a final remark in this section, we would like to compare our statistical uncertainty results with previous studies. By ignoring the event-by-event fluctuation, all events would be identical. Therefore, we can ignore the outer angle brackets in Eq. (45). Also, by collecting many events we can increase the accuracy such that each estimates the true value of the flow harmonic (see Eq. (23)). In such an ideal scenario, Eq. (45) turns into the following form:which has been obtained before in Refs. [36, 64] for the case with .

Fourier-cumulant expansion to study flow distribution

An important property of cumulant analysis is that the lowest order cumulants capture the global features of the p.d.f., and by moving toward higher-orders, more detailed features become relevant. The same is true for Fourier expansion of in Eq. (1) where lower harmonics contain more coarse-grained pictures than higher harmonics. As a result, we consider a double expansion for studying p.d.f. (3): first, harmonic expansion, and second cumulant order expansion. According to this approach, contains the most coarse-grained information (we ignore direct flow in the present study). The harmonic is the first dominant flow, and is the only possible second-order cumulant for . It is well-known that contains information about the initial state ellipticity together with elliptic flow fluctuation. The first third-order cumulant is which means, in addition to second harmonic , the harmonic should be involved as well. Keeping more harmonics and higher-order cumulants, more detailed information about the p.d.f. (3) and can be achieved. The number of cumulants increases rapidly by including more harmonics and keeping higher-order cumulants.

Henceforth, we study harmonics and cumulants up to fifth order, which includes 33 different cumulants as they are tabulated in Table 1. Apart from those cumulants that contain only one moment, there are cumulants containing more than one event-plane moment (see lines 21, 24, 25, 26, 27, 29, and 32 in Table 1) which are studied here for the first time. We should point out that three different five-particle cumulants have been measured by ALICE [65, 66] two of them (see [66]) contain more than one event-plane moment. The direct flow (flow harmonics with ) is involved in all these cumulants. The cumulant (see Ref. [44] as well) in line 22 is equivalent with measured by CMS up to a normalization factor [67].

Normalized cumulants from initial and final states

To understand the observed cumulants’ origin, let us recall that the central part of the flow fluctuations comes from the initial state. To study the initial state fluctuation, we need to quantify the shape of the initial state with quantities similar to those done for flow harmonics, namely eccentricities [68],where is the averaging with respect to the energy density . The true initial energy density shape is captured by the cumulants of the distribution not its moments (eccentricities) in Eq. (47) [53]. The cumulants of the energy density have been obtained by Teaney and Yan in the seminal paper [53] and follow-up papers [25, 69]. Employing the convention in Ref. [69], the two-dimensional energy density cumulants in the polar coordinate are given by One notes that the eccentricities for and 3 are cumulants as well, while for there are contributions from lower harmonics eccentricities and radial shape of the energy density (see Eqs. (2.9) and (2.11) in Ref. [69]).

The anisotropic flow in the final state is a collective response to the initial anisotropic geometry. This process is formulated via the following response relation [53],Up to the leading order, the physics of collectivity and freezeout is encoded in the linear coupling . Ignoring the nonlinear terms, the flow harmonic cumulants of the p.d.f (3) should be proportional to the cumulants of the p.d.f. . The cumulants of the latter p.d.f. are obtained precisely similar to those studied so far by replacing with .5 To be more precise, let us recall the important properties of cumulants, homogeneity. Based on this property, by scaling the random variable as , the cumulants defined in Eq. (5) are scaled asOne can examine the above relation in two explicit examples shown in Eq. (6). By ignoring the nonlinear parts in Eq. (49), we can relate the cumulants obtained from eccentricities and flow harmonics asThe above relation indicates that the initial state and final state fluctuation cumulants differ by a numerical factor, similar to that mentioned in Eq. (50).

For comparing the initial and final state cumulants in the linear response approximation, the values of are needed. However, it is still possible to modify our cumulant definition such that it is independent of numerical value. To this end, we define normalized (standardized) cumulants,Referring to Eqs. (51) and (52), one finds that the normalized cumulants of the initial anisotropy and flow harmonics fluctuations must be precisely the same at the linear approximation. Any deviation between two normalized cumulants should be sourced from nonlinear terms.

The definition of normalized cumulants in Eq. (52) is compatible with the suggestion in Refs. [25, 70] made for event-plane correlations. According to the above definition, for , we havecompatible with Ref. [71]. Also for , we obtainwhere is the normalized symmetric cumulant.

There is an alternative way of defining the normalized cumulant,This definition is compatible with the scalar product method mentioned in Ref. [30]. For cumulants with all , two approaches coincide. In the present study, we use the first convention in Eq. (52) when all cumulants are normalized with the first single-harmonics cumulants. In comparison between the LHC data and simulation, we employ the alternative approach only when the first approach is not measured.

Normalized cumulants from simulation and the LHC

We study the first 29 normalized cumulants in the Fourier-cumulant expansion by using a hybrid hydrodynamic model, TENTo [47]+VISH2+1 [48, 49]+UrQMD [50, 51] ( by definition). We simulate events for Pb–Pb collisions at TeV. The global Bayesian analysis calibrates the model’s free parameters (consist of the temperature dependence of shear and bulk viscosity over entropy density) to explain the measurements by ALICE experiment [52]. In Fig. 1, the normalized version of flow harmonic cumulants listed in Table 1 are plotted with red filled circles. We kept a fixed range for the vertical axes in all panels to simplify the magnitude comparison of different normalized cumulants. In the calibration, the cumulants for , and are used. This means that only the information of normalized cumulant (in panel (7)) has been used for calibration and the rest of the simulations are predictions.

Fig. 1

All 2, 3, 4, 5-order normalized cumulants for harmonics . Magenta squares and open blue circles are eccentricity and energy density cumulant fluctuations from TENTo initial state model, red filled circles are flow fluctuations after hydrodynamic evolution (VISH2+1+UrQMD). The black triangles are the initial state with a linear and nonlinear hydrodynamic response. In the lower-left part of each panel, the corresponding cumulant line number in Table 1 are presented

A series of normalized cumulants in Fig. 1 has already been measured at the LHC (see Table 2 for a summary). The normalized symmetric cumulants, , are measured by ALICE [20, 37, 38], ATLAS [21, 22], and CMS [72] experiments. The event-plane correlations via scalar product method for Pb–Pb collisions at TeV are measured by ATLAS experiment [30]. The measurement of the normalized version (with alternative convention in Eq. (55)) of cumulants , , , and can be found in this analysis.6 The alternative normalized version of is also measured by ATLAS experiment for Pb–Pb collisions at TeV [21] and TeV [22]. And finally, the normalized cumulants for are measured at TeV by ATLAS experiment [22].

Table 2

Measured normalized cumulants at the LHC for Pb–Pb collisions. The far left column indicates the line number of the corresponding cumulant in Table 1. In the references shown by [...], the alternative version of normalized cumulants (see Eq. (55)) are presented. The method of the measurements, scalar product (SP), standard multiparticle (ST), and subevent multiparticle (SE) methods, are indicated for each reference

Normalized cumulant		ALICE		ATLAS		CMS
		2.76 TeV	5.02 TeV	2.76 TeV	5.02 TeV	2.76 TeV	5.02 TeV
5	nc2,4{4}{2,1}			[30]^′(SP), [21]^′(SE)	[22]^′(ST,SE)
6	nc2,3,5{-3,5}{1,1,1}			[30](SP)
7	nc_2{4}				[22](ST,SE)
8	nc_3{4}				[22](ST,SE)
9	nc_4{4}				[22](ST,SE)
11	nc2,3{0}{2,2}	[37](ST)	[20](SE)	[21](SE)	[22](ST,SE)		[72](ST)
12	nc2,4{0}{2,2}	[37](ST)	[20](SE)	[21](SE)	[22](ST,SE)		[72](ST)
13	nc2,5{0}{2,2}	[38](ST)
14	nc3,4{0}{2,2}	[38](ST)
15	nc3,5{0}{2,2}	[38](ST)
17	nc2,3,4{6,-4}{1,2,1}			[30]^′(SP)
20	nc2,3{6}{3,2}			[30]^′(SP)

Although the simulation is calibrated with ALICE experiment kinematics ( and ), the normalized versions of cumulants are less sensitive to these kinematics. In particular, the transverse momentum range dependence of the normalized cumulants has been examined in Refs. [22, 38]. As a result, in addition to the ALICE experiment measurements, we compare our simulation with ATLAS measurements (GeV and ) for Pb–Pb collisions at TeV. The results are depicted in Fig. 2. In general, the simulations and experimental measurements are in agreement. As seen in the figure, the simulation for cumulants and can accurately explain the ATLAS measurements. The ALICE measurement for normalized symmetric cumulant is overestimated by simulation at centralities below . The same is true for around centrality . However, the normalized cumulant is underestimated for centralities above . The simulation at centralities below 10% reveals a sign change for that is not observed by ALICE measurement.

Fig. 2

The (alternative) normalized cumulants that measured by ALICE [37, 38] and ATLAS [30] compared with TENTo+VISH2+1+UrQMD

Except which is used for calibration, the rest of the normalized cumulants shown in Fig. 1 and not displayed in Fig. 2 contain an independent piece of information that can be examined in the future experiment. Here, we only focus on the Pb–Pb collisions at fixed center-of-mass energy. The cumulants’ system size and energy dependence would lead to interesting information about the underlying initial state fluctuation and collective evolution in the future. Referring to Fig. 1, we find that the normalized symmetric cumulants are not the most pronounced ones compared to the others. For instance, the fourth-order normalized cumulant (panel (19)), and fifth-order normalized cumulants ,7 and (panels (22), (27), and (30)) are the largest cumulants up to harmonics that have not been measured yet.

Concerning the feasibility of the cumulant measurements, there is a resent analysis by ALICE collaboration measuring higher-order normalized symmetric cumulants (shown by NSC) consisting flow amplitudes , , and , , [41]. To calculate these sixth-order cumulants, the moments as and are needed to be measured. Given that these cumulants are one order higher than those we presented in Table 1, we would expect that all the unmeasured normalized cumulants in Table 1 are experimentally accessible in the near future. The cumulant (see line25 in Table 1) contains moments with combination. Compared to NSC(2, 3, 5), it has one power higher for flow amplitude for and one power lower for . Considering its large signal from the simulation, we expect that an experimental observation for this cumulant as well.

So far, we have studied the normalized cumulants of Pb–Pb collisions and investigate to what extend the Monte Carlo simulation can explain the experimental data. To better understand both initial state and collective evolution, we have also compared the initial and final state fluctuation cumulants. In the next section, we investigate it in more detail. As we have discussed earlier, up to the linear order of hydrodynamic response, the normalized cumulants of flow harmonics and energy density cumulant fluctuations must be the same. As a result, the observed difference in Fig. 1 is due to the nonlinear hydrodynamic response.

Linear and nonlinear hydrodynamic response coefficients

In Sec. 4.1, it has been discussed that by comparing the normalized cumulants of the initial and final state fluctuation, we can obtain information about the collective hydrodynamic evolution. In the present section, we quantitatively study the hydrodynamic response coefficients.

In the present study, we keep only the linear terms for , and the first subleasing nonlinear terms for . In particular, we consider the following explicit form of Eq. (49), There is another convention in studying hydrodynamic response (see Ref. [68]) in which the eccentricities are employed to present the hydrodynamic response, We will discuss the interpretation of these two conventions later.

Different approaches can be employed to extract the couplings , , and (or , , and ). We can start from a Gaussian geometrical initial energy density and deform it with one (or a few numbers) of nonvanishing . Then we change the value of the energy density cumulant to probe the hydrodynamic response from this “single shot” simulation [25, 69]. It is shown that (Marcinkiewicz theorem8) reproducing a distribution (except Gaussian distribution) with a finite number of cumulants leads to negative values at some parts of distribution domain. These negative values should be regulated, which consequently produces spurious cumulants [69]. Until spurious cumulants are small, we can neglect their effect in the coupling estimations. The initial state with larger deformation would have more such a problem.

As another approach, one can generate many events with complicated initial geometries and study the hydrodynamic response. For instance, to extract the response coefficient in Refs. [68, 74], the authors multiply both sides of Eq. (49) (ignore nonlinear terms) with and then average over many simulated events in a given centrality class. This approach leads to the following estimation: . By writing this estimator in the Cartesian coordinate, and , we see that this relation is the Pearson correlation between eccentricities and flow harmonics. This estimator is related to the covariance matrix of . To estimate the nonlinear hydrodynamic couplings similar to what is done in Ref. [68], we need to include more harmonics in p.d.f. . Using this approach, we have ignored the non-Gaussian effects encoded in .

Comparing the event-by-event fluctuation of the initial and final states in a simulation is an alternative method (the one we employ here) to study the hydrodynamic response couplings. For that, we technically compare with . By substituting in the flow harmonic fluctuation cumulant (moments) with the corresponding expression for the hydrodynamic response Eqs. (56), we obtain a function written in terms of initial state fluctuation cumulant (moments). For instance, by employing Eq. (56c), we obtainSimilarly, we work out the following flow harmonic cumulants in terms of ,We can rewrite the right-hand sides of the above equations in terms of cumulants as well.

There are six equations (Eqs. (58) and (59)) and six unknown response coupling constants that we can find numerically. We do a naïve analysis by ignoring the statistical errors in our Monte Carlo simulation to extract the unknown coefficients. At the linear level, the initial state’s pressure gradient enforces that the event-plane angle to be the same as the participant plane . As a result, we expect that all s are real and positive. By demanding that the linear response coupling constants , are positive, we obtain only one set of solutions for couplings at each centrality class. The results are depicted in Fig. 3. The couplings and (which are identical with and in our study) are explicitly reported in Ref. [74]. Considering that a different hydrodynamic model (with different tuning) has been used, and a different method (mentioned in Ref. [68]) is employed to extract the couplings, our results for and are compatible with those computed in Ref. [74]. The couplings and are computed in Ref. [25]. The values of these couplings approach zero and change the sign at mid-central collisions. In our method, there are valid positive solutions for and with no sign change. We have also found larger values for these couplings. However, our computations have rather similar behavior as those reported in Ref. [25] for the nonlinear coefficients. The observed differences between the two studies could be due to the different hydrodynamic models or different approaches of extracting the coefficients (single-shot approach has been employed in Ref. [25]). One notes that our estimations depend on the number of nonlinear response terms that we have considered in our hydrodynamic response estimation. By adding more terms, we need to employ more moments (cumulants) as input. This modification can lead to slightly different values for the couplings.

Fig. 3

Linear and nonlinear hydrodynamic response couplings obtained by comparing initial and final state normalized cumulants

As it has been mentioned in Sect. 4.1, to understand that how much flow harmonic fluctuation is originated from the initial state, one can compute the initial state fluctuation normalized cumulants and compare them with those obtained from flow harmonic. In Fig. 1, fluctuation normalized cumulants (blue empty circles) and fluctuation normalized cumulants (magenta squares) are depicted. The observed difference between and cumulants is interpreted as the presence of nonlinear terms similar to what is mentioned in Eqs. (57). The comparison of and gives us a hint about the nonlinear terms similar to what has been shown in Eqs. (56).

We see in the figure that some of the normalized cumulants that are computed from have a different sign from those calculated from (see for instance panels (5), (6), (22), and (30), in Fig. 1). In fact, the “wrong” sign of fluctuation in participant plane correlation has been observed in the previous studied [25, 30, 69, 75–77]. In particular, the quantities and in [77] are similar to and , up to a normalization factor. A sign difference between initial and final state fluctuations has been observed for these quantities. The authors of Ref. [77] conclude that this sign difference is a signature of hydrodynamic response to the initial state. Up to a normalization factor, the same quantities have been studied in [25]. Again a sign change has been observed between eccentricity and flow harmonic fluctuations while the sign change is resolved by replacing with . However, we cannot conclude that replacing with always leads to a compatible correlation sign with final state fluctuation. For instance, one can find cases in panels (18), (23), and (28) of Fig. 1 that correlations have opposite sign compared to flow harmonic fluctuations. Few examples can be found in [25] as well.9

The cumulants calculated from hydrodynamic nonlinear response estimation are shown with black filled triangles in Fig. 1, where the response couplings are those displayed in Fig. 3. The figure shows that the two first cumulants ( and ) have a perfect agreement with the nonlinear response estimation. These two cumulants are inputs to estimate the response couplings. Since normalized cumulants are plotted in Fig. 1, only two ratios and play a role. As a result, two normalized cumulants are enough to fix these ratios. The rest of the normalized cumulants computed from the nonlinear response are, in fact, the predictions. Some of these predictions perfectly match with hydrodynamic computations. For instance, , and has shown in panels (15), (22), and (30) in Fig. 1. There are cases with a poor agreement between hydrodynamic simulation and the nonlinear response estimation, namely and (panels (9) and (10) in Fig. 1). The nonlinear terms could not cure the initial and final correlation sign differences in panels (18), (23), and (28). Including more nonlinear terms would increase the accuracy of the latter cases.

Concerning the interpretation of the nonlinear hydrodynamic response coefficients, we note that the true deformation of the initial energy density is quantified by cumulants . One expects that the final anisotropy to be proportional to the true deformation [53]. Assuming and have solely collective evolution contributions, the nonlinear couplings and in Eq. (57) have contributions from both initial shape and collective evolution. One can relate these two couplings by substituting cumulants (48) into Eqs. (56) and compare them with Eqs. (57). The linear couplings are identical in both conventions, . However, we obtain the following relation for the nonlinear couplings: The above relations mean that the coefficients and receive contributions from the radial shape of the initial energy density as well. In Fig. 4, the coupling ratios and are plotted with square points (blue curves). The couplings and (shown by circles, black curves) are estimated by solving six equations in Eqs. (58) and (59) where energy density cumulants are replaced by eccentricities. Since the coupling ratio estimations are estimated from many events in a given centrality, we compute the right-hand side of Eqs. (60), where the initial energy density contributions are averaged over events in the given centrality class. The quantities and (shown by diamonds, red curve) are directly computed from TENTo events. The estimated and from Eqs. (60) are shown by triangles (green curves) in Fig. 4. We see a rather perfect match between black curves and green curves, as we expected. From the figure, we see a nontrivial centrality dependence of quantities and . This observation indicates that there are contributions from the radial shape of the initial energy density in the values of and . This contribution should be taken into account in interpreting the nonlinear couplings and as hydrodynamic response couplings.

Fig. 4

The TENTo initial state radial shape (diamond points), and the nonlinear over linear hydrodynamic response couplings from VISH2+1+UrQMD (circle and square points)

Flow-induced genuine three-particle correlation

The observation of long-range correlations between particles at and nonzero is one of the first and most important fluid(-like) signals in large (small) system collisions [5, 54–56, 78]. In particular, function , the two-particle correlation function (2PC),10 quantifies the correlation between two particles in the final state. In the present section, we introduce a generalized version of this function and study its relation with flow harmonic cumulants.

Ignoring the experimental complications, in principle, the function is measured as follows: we choose all distinguished pairs of particles in an event and compute and for each pair and fill a histogram from a collection of pairs in many events. To find the connection between 2PC with cumulants, we focus only on the part of the correlation. In other words, we choose all the distinguished particles and compute irrespective of their position in the direction. One can use -gaps (see Refs. [24, 79]) to decrease the contamination of nonflow effects in the correlation function estimation. For two-particle correlation, we find the well-known relation between 2PC and second-order cumulants ,Here, we have used the notation instead of commonly used notation for future generalization. The advantage of measuring compared to is that it contains a cumulative information of all second order cumulants with .

We can extend the notion of two-particle correlation function into q-particle correlation function (qPC). Specifically, in the following, we focus on correlations of three particles in the final state, because it is a function of two variablesand easy to visualize. For the same reason (and finding a clear connection to the flow harmonic cumulants), we ignore the dependence similar to Eq. (61). The -gap method can also be used in this case to decrease the nonflow effects as well. Employing the systematic study of multiharmonic cumulants presented in this manuscript, we can find an expansion of any qPC in terms q-order cumulants, similar to what has been written in Eq. (61). The technical details of finding the relation between qPC and q order cumulants can be found in Appendix C. Here, we show the final result for 3PC where harmonics are involved,The above relation is the correlation function of three particles when all lower-order correlations are removed (see Appendix C.2). The functions and can be found in Eq. (118). By employing and from TENTo+VISH2+1+UrQMD simulation, we have plotted in Fig. 5. Here, we have subtracted to focus only on the nontrivial correlation.

Fig. 5

Genuine three-particle correlation function, 3PC. The correlation is obtained from the T ENTo+VISH2+1+UrQMD simulation for Pb–Pb collision at  TeV in the 10–20% centrality class

As seen in Fig. 5, the correlation reveals repeating patterns which is a consequence of symmetries. Before explaining these symmetries in three-particle correlations, let us discuss them in a more simple case, two-particle correlation functions. To measure two-particle correlation functions, two (charged) particles in a given event are chosen and compute . The signal distribution is obtained by measuring for many events.11 Let us call two distinguished particles in a given event as and with azimuthal angles and . Choosing all pairs of distinguished particles in an event and call them as particle 1 and 2, the particle is labeled as particle 1 and as particle 2 once and the particle as particle 2 and as particle 1 again. The former labeling leads to and the latter to . As a result, in the signal distribution both and have contribution from a single configuration of particles and . This labeling brings us to the following symmetry, . One notes that the mentioned “double counting” is also considered in (see Eq. (22)). As a result, this symmetry is manifestly true in Eq. (61). We also have an obvious periodic symmetry for any integer n. Using these two symmetries, we find that we have independent information in in the range . The correlation (or ) is mostly reported in the range to see a clear ridge and shoulder structures (see for instance [78]). Due to the symmetries explained above, both ridge and shoulder structures are symmetric with respect to the axes pass through their peaks at and and only half of these structures contain independent information in the range .

When we have three particles in the final state, there are 3! possible ways to label three distinguished particles, leading to six different values for with a single configuration of particles. Let us call three particles as , , and with azimuthal angles , , and . If particles , , and are labeled by 1, 2, and 3 (see Eq. (62)), we obtain , and . The particles , and can be labeled as 2, 1, and 3. With this labeling, we have and . Both values for and refer to the same configuration of particles and exist in the signal distribution. As a result, the final distributions have the following symmetry: . Another case is labeling , , and particles with 1, 3, and 2. This labeling leads to the values and . Consequently, the distribution has symmetry . There are three other permutations, but all of them lead to a combination of symmetries explained above. More than the above symmetries, for any configuration of particles, we expect the same probability for the mirrored configurations, , which leads to the symmetry . Finally, we have rotational symmetry for any integer , which eventually leads to the periodic condition for with periodicity . We summarize all the symmetries of 3PC as follows,The symmetries in Eq. (64) are responsible of repeating patterns observed in Fig. 5. Each dashed line in the figure refers to one (or combination of more than one) symmetry(ies). As seen in the figure, the range is divided into twelve repeating triangular regions. In the experiment and the simulations, we lose the statistics twelve times by considering the full region without gaining any more information. For that reason, we confine ourselves into a “unit cell” of and containing all the nontrivial non-repeating patterns. Here, we conventionally choose the following region,as it is shown by a black triangle in Fig. 5. If any choice of particles leads to outside this triangle, the symmetries mentioned in Eq. (64) can be used to map that point inside it. The edges of the black triangle in Fig. 5 refer to the combinations of the angles, and with unequal lengths. It will have more symmetric visualization if we use the following variables,where and . Using this variables, we can plot the unit cell of in a Dalitz-like plot as it is shown in Fig. 6.

Fig. 6

Dalitz-like plot of the “unit-cell” (black triangle in Fig. 5) of genuine three-particle correlation function

Similar to which has cumulative information from all second-order cumulants, experimental measurement of (or equivalently ) contains a cumulative information of all third-order cumulants at the same time. Measuring these correlations in large and small systems can be used as an independent method of testing the collectivity and event-by-event fluctuation in heavy-ion experiments. We leave more investigation of this measurable quantity to future studies.

Conclusion

A multidimensional generating function method was introduced to extract a large class of cumulants related to the flow harmonic fluctuations. We proposed an ordering based on the Fourier-cumulant expansion for these observables to systematically capture the most dominant features of the flow harmonic fluctuations. Using this method and reproducing the already known cumulants, we have found new cumulants consisting of symmetry plane correlations that have not been studied before. We defined the normalized cumulants to compare the cumulant’s magnitude with each other and with initial state fluctuation. We employed hydrodynamic simulation for Pb–Pb collisions (TENTo+VISH2+1+UrQMD) calibrated by a Global Bayesian analysis to predict the unmeasured normalized cumulants’ value. The observables introduced in this study can be used as inputs for calibration of the Bayesian analysis or to validate the already tuned parameters.

We extract the linear and nonlinear hydrodynamic response by comparing the initial anisotropy and flow harmonic fluctuation cumulants. This method can be extended to obtain higher-order nonlinear terms by comparing more initial and final state fluctuations cumulants.

Based on the Fourier-cumulant expansion, we also introduced a general way to find the genuine q-particle distribution function, qPC. In particular, we studied 3PC containing information of all third-order cumulants for all harmonics. The flow-induced 3PC for harmonics was presented by using TENTo+VISH2+1+UrQMD simulation. This measurable helps to study flow and nonflow effects in large and small systems in the future.

{m_1,m_2}	{α}	c2,4{α}{m_1,m_2}
{2,1}	{4}	⟨v22v_4cos4ψ_2-ψ_4⟩
{2,2}	{0}	⟨v22v42⟩-⟨v22⟩⟨v42⟩
{2,3}	{4}	⟨v22v43cos4ψ_2-ψ_4⟩-2⟨v42⟩⟨v22v_4cos4ψ_2-ψ_4⟩
{4,1}	{4}	⟨v24v_4cos4ψ_2-ψ_4⟩-3⟨v22⟩⟨v22v_4cos4ψ_2-ψ_4⟩
{4,2}	{8}	⟨v24v42cos8ψ_2-ψ_4⟩-6⟨v22v_4cos4ψ_2-ψ_4⟩^2
{2,4}	{0}	⟨v22v44⟩-4⟨v42⟩⟨v22v42⟩-⟨v22⟩⟨v44⟩+4⟨v22⟩⟨v42⟩^2
{4,2}	{0}	⟨v24v42⟩-⟨v22v_4cos4ψ_2-ψ_4⟩^2-⟨v24⟩⟨v42⟩-4⟨v22⟩⟨v22v42⟩+4⟨v22⟩^2⟨v42⟩
{2,5}	{4}	⟨v22v45cos4ψ_2-ψ_4⟩-6⟨v42⟩⟨v22v43cos4ψ_2-ψ_4⟩-3⟨v22v_4cos4ψ_2-ψ_4⟩⟨v44⟩+12⟨v22v_4cos4ψ_2-ψ_4⟩⟨v42⟩
{4,3}	{4}	⟨v24v43cos4ψ_2-ψ_4⟩-3⟨v22⟩⟨v22v43cos4ψ_2-ψ_4⟩-2⟨v24v_4cos4ψ_2-ψ_4⟩⟨v42⟩-6⟨v22v_4cos4ψ_2-ψ_4⟩⟨v22v42⟩+12⟨v22⟩⟨v22v_4cos4ψ_2-ψ_4⟩⟨v42⟩
{6,1}	{4}	⟨v26v_4cos4ψ_2-ψ_4⟩-6⟨v24⟩⟨v22v_4cos4ψ_2-ψ_4⟩-8⟨v22⟩⟨v24v_4cos4ψ_2-ψ_4⟩+24⟨v22⟩^2⟨v22v_4cos4ψ_2-ψ_4⟩
{4,4}	{8}	⟨v44v24cos8ψ_2-ψ_4⟩-18⟨v22v_4cos4ψ_2-ψ_4⟩⟨v22v43cos4ψ_2-ψ_4⟩-3⟨v42⟩⟨v24v42cos8ψ_2-ψ_4⟩+36⟨v22v_4cos4ψ_2-ψ_4⟩^2⟨v42⟩
{6,2}	{8}	⟨v26v42cos8ψ_2-ψ_4⟩-20⟨v22v_4cos4ψ_2-ψ_4⟩⟨v24v_4cos4ψ_2-ψ_4⟩-5⟨v22⟩⟨v24v42cos8ψ_2-ψ_4⟩+60⟨v22⟩⟨v22v_4cos4ψ_2-ψ_4⟩^2
{2,6}	{0}	⟨v22v46⟩-9⟨v22v42⟩⟨v44⟩-9⟨v42⟩⟨v22v44⟩-⟨v22⟩⟨v46⟩+36⟨v42⟩^2⟨v22v42⟩+18⟨v22⟩⟨v42⟩⟨v44⟩-36⟨v22⟩⟨v42⟩^3
{4,4}	{0}	⟨v24v44⟩-8⟨v22v42⟩^2-4⟨v42⟩⟨v24v42⟩-⟨v24⟩⟨v44⟩-4⟨v22⟩⟨v22v44⟩+4⟨v24⟩⟨v42⟩^2+32⟨v22⟩⟨v42⟩⟨v22v42⟩+4⟨v22⟩^2⟨v44⟩-24⟨v22⟩^2⟨v42⟩^2-4⟨v22v_4cos4ψ_2-ψ_4⟩⟨v22v43cos4ψ_2-ψ_4⟩+8⟨v22v_4cos4ψ_2-ψ_4⟩^2⟨v42⟩
{6,2}	{0}	⟨v26v42⟩-⟨v26⟩⟨v42⟩-9⟨v24⟩⟨v22v42⟩-9⟨v22⟩⟨v24v42⟩+18⟨v22⟩⟨v24⟩⟨v42⟩+36⟨v22⟩^2⟨v22v42⟩-36⟨v22⟩^3⟨v42⟩-6⟨v22v_4cos4ψ_2-ψ_4⟩⟨v24v_4cos4ψ_2-ψ_4⟩+18⟨v22⟩⟨v22v_4cos4ψ_2-ψ_4⟩^2

{m_1,m_2}	{α}	c2,5{α}{m_1,m_2}
{2,2}	{0}	⟨v22v52⟩-⟨v22⟩⟨v52⟩
{2,4}	{0}	⟨v22v54⟩-4⟨v52⟩⟨v22v52⟩-⟨v22⟩⟨v54⟩+4⟨v22⟩⟨v52⟩^2
{4,2}	{0}	⟨v24v52⟩-⟨v24⟩⟨v52⟩-4⟨v22⟩⟨v22v52⟩+4⟨v22⟩^2⟨v52⟩
{5,2}	{10}	⟨v25v52cos10ψ_2-ψ_5⟩
{2,6}	{0}	⟨v22v56⟩-9⟨v22v52⟩⟨v54⟩-9⟨v52⟩⟨v22v54⟩-⟨v22⟩⟨v56⟩+36⟨v52⟩^2⟨v22v52⟩+18⟨v22⟩⟨v52⟩⟨v54⟩-36⟨v22⟩⟨v52⟩^3
{4,4}	{0}	⟨v24v54⟩-8⟨v22v52⟩^2-4⟨v52⟩⟨v24v52⟩-⟨v24⟩⟨v54⟩-4⟨v22⟩⟨v22v54⟩+4⟨v24⟩⟨v52⟩^2+32⟨v22⟩⟨v52⟩⟨v22v52⟩+4⟨v22⟩^2⟨v54⟩-24⟨v22⟩^2⟨v52⟩^2
{6,2}	{0}	⟨v26v52⟩-⟨v26⟩⟨v52⟩-9⟨v24⟩⟨v22v52⟩-9⟨v22⟩⟨v24v52⟩+18⟨v22⟩⟨v24⟩⟨v52⟩+36⟨v22⟩^2⟨v22v52⟩-36⟨v22⟩^3⟨v52⟩

{m_1,m_2}	{α}	c2,6{α}{m_1,m_2}
{3,1}	{6}	⟨v23v_6cos6ψ_2-ψ_6⟩
{2,2}	{0}	⟨v22v62⟩-⟨v22⟩⟨v62⟩
{3,3}	{6}	⟨v63v23cos6ψ_2-ψ_6⟩-2⟨v23v_6cos6ψ_2-ψ_6⟩⟨v62⟩
{5,1}	{6}	⟨v25v_6cos6ψ_2-ψ_6⟩-4⟨v22⟩⟨v23v_6cos6ψ_2-ψ_6⟩
{2,4}	{0}	⟨v22v64⟩-4⟨v62⟩⟨v22v62⟩-⟨v22⟩⟨v64⟩+4⟨v22⟩⟨v62⟩^2
{4,2}	{0}	⟨v24v62⟩-⟨v24⟩⟨v62⟩-4⟨v22⟩⟨v22v62⟩+4⟨v22⟩^2⟨v62⟩
{6,2}	{12}	⟨v26v62cos12ψ_2-ψ_6⟩-20⟨v23v_6cos6ψ_2-ψ_6⟩^2
{3,5}	{6}	⟨v23v65cos6ψ_2-ψ_6⟩-6⟨v62⟩⟨v63v23cos6ψ_2-ψ_6⟩-3⟨v23v_6cos6ψ_2-ψ_6⟩⟨v64⟩+12⟨v23v_6cos6ψ_2-ψ_6⟩⟨v62⟩^2
{5,3}	{6}	⟨v25v63cos6ψ_2-ψ_6⟩-4⟨v22⟩⟨v63v23cos6ψ_2-ψ_6⟩-2⟨v25v_6cos6ψ_2-ψ_6⟩⟨v62⟩-8⟨v23v_6cos6ψ_2-ψ_6⟩⟨v22v62⟩+16⟨v22⟩⟨v23v_6cos6ψ_2-ψ_6⟩⟨v62⟩
{2,6}	{0}	⟨v22v66⟩-9⟨v22v62⟩⟨v64⟩-9⟨v62⟩⟨v22v64⟩-⟨v22⟩⟨v66⟩+36⟨v62⟩^2⟨v22v62⟩+18⟨v22⟩⟨v62⟩⟨v64⟩-36⟨v22⟩⟨v62⟩^3
{7,1}	{6}	⟨v27v_6cos6ψ_2-ψ_6⟩-10⟨v24⟩⟨v23v_6cos6ψ_2-ψ_6⟩-10⟨v22⟩⟨v25v_6cos6ψ_2-ψ_6⟩+40⟨v22⟩^2⟨v23v_6cos6ψ_2-ψ_6⟩
{4,4}	{0}	⟨v24v64⟩-8⟨v22v62⟩^2-4⟨v62⟩⟨v24v62⟩-⟨v24⟩⟨v64⟩-4⟨v22⟩⟨v22v64⟩+4⟨v24⟩⟨v62⟩^2+32⟨v22⟩⟨v62⟩⟨v22v62⟩+4⟨v22⟩^2⟨v64⟩-24⟨v22⟩^2⟨v62⟩^2
{6,2}	{0}	⟨v26v62⟩-⟨v23v_6cos6ψ_2-ψ_6⟩^2-⟨v26⟩⟨v62⟩-9⟨v24⟩⟨v22v62⟩-9⟨v22⟩⟨v24v62⟩+18⟨v22⟩⟨v24⟩⟨v62⟩+36⟨v22⟩^2⟨v22v62⟩-36⟨v22⟩^3⟨v62⟩

{m_1,m_2}	{α}	c3,4{α}{m_1,m_2}
{2,2}	{0}	⟨v32v42⟩-⟨v32⟩⟨v42⟩
{2,4}	{0}	⟨v32v44⟩-4⟨v42⟩⟨v32v42⟩-⟨v32⟩⟨v44⟩+4⟨v32⟩⟨v42⟩^2
{4,2}	{0}	⟨v34v42⟩-⟨v34⟩⟨v42⟩-4⟨v32⟩⟨v32v42⟩+4⟨v32⟩^2⟨v42⟩
{4,3}	{12}	⟨v34v43cos12ψ_3-ψ_4⟩
{2,6}	{0}	⟨v32v46⟩-9⟨v32v42⟩⟨v44⟩-9⟨v42⟩⟨v32v44⟩-⟨v32⟩⟨v46⟩+36⟨v42⟩^2⟨v32v42⟩+18⟨v32⟩⟨v42⟩⟨v44⟩-36⟨v32⟩⟨v42⟩^3
{4,4}	{0}	⟨v34v44⟩-8⟨v32v42⟩^2-4⟨v42⟩⟨v34v42⟩-⟨v34⟩⟨v44⟩-4⟨v32⟩⟨v32v44⟩+4⟨v34⟩⟨v42⟩^2+32⟨v32⟩⟨v42⟩⟨v32v42⟩+4⟨v32⟩^2⟨v44⟩-24⟨v32⟩^2⟨v42⟩^2
{6,2}	{0}	⟨v36v42⟩-⟨v36⟩⟨v42⟩-9⟨v34⟩⟨v32v42⟩-9⟨v32⟩⟨v34v42⟩+18⟨v32⟩⟨v34⟩⟨v42⟩+36⟨v32⟩^2⟨v32v42⟩-36⟨v32⟩^3⟨v42⟩

{m_1,m_2}	{α}	c3,5{α}{m_1,m_2}
{2,2}	{0}	⟨v32v52⟩-⟨v32⟩⟨v52⟩
{2,4}	{0}	⟨v32v54⟩-4⟨v52⟩⟨v32v52⟩-⟨v32⟩⟨v54⟩+4⟨v32⟩⟨v52⟩^2
{4,2}	{0}	⟨v34v52⟩-⟨v34⟩⟨v52⟩-4⟨v32⟩⟨v32v52⟩+4⟨v32⟩^2⟨v52⟩
{5,3}	{15}	⟨v35v53cos15ψ_3-ψ_5⟩
{2,6}	{0}	⟨v32v56⟩-9⟨v32v52⟩⟨v54⟩-9⟨v52⟩⟨v32v54⟩-⟨v32⟩⟨v56⟩+36⟨v52⟩^2⟨v32v52⟩+18⟨v32⟩⟨v52⟩⟨v54⟩-36⟨v32⟩⟨v52⟩
{4,4}	{0}	⟨v34v54⟩-8⟨v32v52⟩^2-4⟨v52⟩⟨v34v52⟩-⟨v34⟩⟨v54⟩-4⟨v32⟩⟨v32v54⟩+4⟨v34⟩⟨v52⟩^2+32⟨v32⟩⟨v52⟩⟨v32v52⟩+4⟨v32⟩^2⟨v54⟩-24⟨v32⟩^2⟨v52⟩^2
{6,2}	{0}	⟨v36v52⟩-⟨v36⟩⟨v52⟩-9⟨v34⟩⟨v32v52⟩-9⟨v32⟩⟨v34v52⟩+18⟨v32⟩⟨v34⟩⟨v52⟩+36⟨v32⟩^2⟨v32v52⟩-36⟨v32⟩^3⟨v52⟩

{m_1,m_2}	{α}	c3,6{α}{m_1,m_2}
{2,1}	{6}	⟨v32v_6cos6ψ_3-ψ_6⟩
{2,2}	{0}	⟨v32v62⟩-⟨v32⟩⟨v62⟩
{2,3}	{6}	⟨v32v63cos6ψ_3-ψ_6⟩-2⟨v32v_6cos6ψ_3-ψ_6⟩⟨v62⟩
{4,1}	{6}	⟨v34v_6cos6ψ_3-ψ_6⟩-3⟨v32⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{4,2}	{12}	⟨v34v62cos12ψ_3-ψ_6⟩-6⟨v32v_6cos6ψ_3-ψ_6⟩^2
{2,4}	{0}	⟨v32v64⟩-4⟨v62⟩⟨v32v62⟩-⟨v32⟩⟨v64⟩+4⟨v32⟩⟨v62⟩^2
{4,2}	{0}	⟨v34v62⟩-⟨v32v_6cos6ψ_3-ψ_6⟩^2-⟨v34⟩⟨v62⟩-4⟨v32⟩⟨v32v62⟩+4⟨v32⟩^2⟨v62⟩
{2,5}	{6}	⟨v32v65cos6ψ_3-ψ_6⟩-6⟨v62⟩⟨v32v63cos6ψ_3-ψ_6⟩-3⟨v32v_6cos6ψ_3-ψ_6⟩⟨v64⟩+12⟨v32v_6cos6ψ_3-ψ_6⟩⟨v62⟩^2
{4,3}	{6}	⟨v34v63cos6ψ_3-ψ_6⟩-3⟨v32⟩⟨v32v63cos6ψ_3-ψ_6⟩-2⟨v34v_6cos6ψ_3-ψ_6⟩⟨v62⟩-6⟨v32v_6cos6ψ_3-ψ_6⟩⟨v32v62⟩+12⟨v32⟩⟨v32v_6cos6ψ_3-ψ_6⟩⟨v62⟩
{6,1}	{6}	⟨v36v_6cos6ψ_3-ψ_6⟩-6⟨v34⟩⟨v32v_6cos6ψ_3-ψ_6⟩-8⟨v32⟩⟨v34v_6cos6ψ_3-ψ_6⟩+24⟨v32⟩^2⟨v32v_6cos6ψ_3-ψ_6⟩
{4,4}	{12}	⟨v64v34cos12ψ_3-ψ_6⟩-18⟨v32v_6cos6ψ_3-ψ_6⟩⟨v32v63cos6ψ_3-ψ_6⟩-3⟨v62⟩⟨v34v62cos12ψ_3-ψ_6⟩+36⟨v32v_6cos6ψ_3-ψ_6⟩^2⟨v62⟩
{6,2}	{12}	⟨v36v62cos12ψ_3-ψ_6⟩-20⟨v32v_6cos6ψ_3-ψ_6⟩⟨v34v_6cos6ψ_3-ψ_6⟩-5⟨v32⟩⟨v34v62cos12ψ_3-ψ_6⟩+60⟨v32⟩⟨v32v_6cos6ψ_3-ψ_6⟩^2
{2,6}	{0}	⟨v32v66⟩-9⟨v32v62⟩⟨v64⟩-9⟨v62⟩⟨v32v64⟩-⟨v32⟩⟨v66⟩+36⟨v62⟩^2⟨v32v62⟩+18⟨v32⟩⟨v62⟩⟨v64⟩-36⟨v32⟩⟨v62⟩^3
{4,4}	{0}	⟨v34v64⟩-8⟨v32v62⟩^2-4⟨v62⟩⟨v34v62⟩-⟨v34⟩⟨v64⟩-4⟨v32⟩⟨v32v64⟩+4⟨v34⟩⟨v62⟩^2+32⟨v32⟩⟨v62⟩⟨v32v62⟩+4⟨v32⟩^2⟨v64⟩-24⟨v32⟩^2⟨v62⟩^2-4⟨v32v_6cos6ψ_3-ψ_6⟩⟨v32v63cos6ψ_3-ψ_6⟩+8⟨v32v_6cos6ψ_3-ψ_6⟩^2⟨v62⟩
{6,2}	{0}	⟨v36v62⟩-⟨v36⟩⟨v62⟩-9⟨v34⟩⟨v32v62⟩-9⟨v32⟩⟨v34v62⟩+18⟨v32⟩⟨v34⟩⟨v62⟩+36⟨v32⟩^2⟨v32v62⟩-36⟨v32⟩^3⟨v62⟩-6⟨v32v_6cos6ψ_3-ψ_6⟩⟨v34v_6cos6ψ_3-ψ_6⟩+18⟨v32⟩⟨v32v_6cos6ψ_3-ψ_6⟩^2

{m_1,m_2}	{α}	c4,5{α}{m_1,m_2}
{2,2}	{0}	⟨v42v52⟩-⟨v42⟩⟨v52⟩
{2,4}	{0}	⟨v42v54⟩-4⟨v52⟩⟨v42v52⟩-⟨v42⟩⟨v54⟩+4⟨v42⟩⟨v52⟩^2
{4,2}	{0}	⟨v44v52⟩-⟨v44⟩⟨v52⟩-4⟨v42⟩⟨v42v52⟩+4⟨v42⟩^2⟨v52⟩
{2,6}	{0}	⟨v42v56⟩-9⟨v42v52⟩⟨v54⟩-9⟨v52⟩⟨v42v54⟩-⟨v42⟩⟨v56⟩+36⟨v52⟩^2⟨v42v52⟩+18⟨v42⟩⟨v52⟩⟨v54⟩-36⟨v42⟩⟨v52⟩^3
{4,4}	{0}	⟨v44v54⟩-8⟨v42v52⟩^2-4⟨v52⟩⟨v44v52⟩-⟨v44⟩⟨v54⟩-4⟨v42⟩⟨v42v54⟩+4⟨v44⟩⟨v52⟩^2+32⟨v42⟩⟨v52⟩⟨v42v52⟩+4⟨v42⟩^2⟨v54⟩-24⟨v42⟩^2⟨v52⟩^2
{6,2}	{0}	⟨v46v52⟩-⟨v46⟩⟨v52⟩-9⟨v44⟩⟨v42v52⟩-9⟨v42⟩⟨v44v52⟩+18⟨v42⟩⟨v44⟩⟨v52⟩+36⟨v42⟩^2⟨v42v52⟩-36⟨v42⟩^3⟨v52⟩

{m_1,m_2}	{α}	c4,6{α}{m_1,m_2}
{2,2}	{0}	⟨v42v62⟩-⟨v42⟩⟨v62⟩
{3,2}	{12}	⟨v43v62cos12ψ_4-ψ_6⟩
{2,4}	{0}	⟨v42v64⟩-4⟨v62⟩⟨v42v62⟩-⟨v42⟩⟨v64⟩+4⟨v42⟩⟨v62⟩^2
{4,2}	{0}	⟨v44v62⟩-⟨v44⟩⟨v62⟩-4⟨v42⟩⟨v42v62⟩+4⟨v42⟩^2⟨v62⟩
{3,4}	{12}	⟨v43v64cos12ψ_4-ψ_6⟩-3⟨v62⟩⟨v43v62cos12ψ_4-ψ_6⟩
{5,2}	{12}	⟨v45v62cos12ψ_4-ψ_6⟩-4⟨v42⟩⟨v43v62cos12ψ_4-ψ_6⟩
{2,6}	{0}	⟨v42v66⟩-9⟨v42v62⟩⟨v64⟩-9⟨v62⟩⟨v42v64⟩-⟨v42⟩⟨v66⟩+36⟨v62⟩^2⟨v42v62⟩+18⟨v42⟩⟨v62⟩⟨v64⟩-36⟨v42⟩⟨v62⟩^3
{4,4}	{0}	⟨v44v64⟩-8⟨v42v62⟩^2-4⟨v62⟩⟨v44v62⟩-⟨v44⟩⟨v64⟩-4⟨v42⟩⟨v42v64⟩+4⟨v44⟩⟨v62⟩^2+32⟨v42⟩⟨v62⟩⟨v42v62⟩+4⟨v42⟩^2⟨v64⟩-24⟨v42⟩^2⟨v62⟩^2
{6,2}	{0}	⟨v46v62⟩-⟨v46⟩⟨v62⟩-9⟨v44⟩⟨v42v62⟩-9⟨v42⟩⟨v44v62⟩+18⟨v42⟩⟨v44⟩⟨v62⟩+36⟨v42⟩^2⟨v42v62⟩-36⟨v42⟩^3⟨v62⟩

{m_1,m_2}	{α}	c5,6{α}{m_1,m_2}
{2,2}	{0}	⟨v52v62⟩-⟨v52⟩⟨v62⟩
{2,4}	{0}	⟨v52v64⟩-4⟨v62⟩⟨v52v62⟩-⟨v52⟩⟨v64⟩+4⟨v52⟩⟨v62⟩^2
{4,2}	{0}	⟨v54v62⟩-⟨v54⟩⟨v62⟩-4⟨v52⟩⟨v52v62⟩+4⟨v52⟩^2⟨v62⟩
{2,6}	{0}	⟨v52v66⟩-9⟨v52v62⟩⟨v64⟩-9⟨v62⟩⟨v52v64⟩-⟨v52⟩⟨v66⟩+36⟨v62⟩^2⟨v52v62⟩+18⟨v52⟩⟨v62⟩⟨v64⟩-36⟨v52⟩⟨v62⟩^3
{4,4}	{0}	⟨v54v64⟩-8⟨v52v62⟩^2-4⟨v62⟩⟨v54v62⟩-⟨v54⟩⟨v64⟩-4⟨v52⟩⟨v52v64⟩+4⟨v54⟩⟨v62⟩^2+32⟨v52⟩⟨v62⟩⟨v52v62⟩+4⟨v52⟩^2⟨v64⟩-24⟨v52⟩^2⟨v62⟩^2
{6,2}	{0}	⟨v56v62⟩-⟨v56⟩⟨v62⟩-9⟨v54⟩⟨v52v62⟩-9⟨v52⟩⟨v54v62⟩+18⟨v52⟩⟨v54⟩⟨v62⟩+36⟨v52⟩^2⟨v52v62⟩-36⟨v52⟩^3⟨v62⟩

{m_1,m_2,m_3}	{α_1,α_2}	c2,3,5{α_1,α_2}{m_1,m_2,m_3}
{1,1,1}	{-3,5}	⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{1,1,3}	{-3,5}	⟨v53v_2v_3cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v52⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{1,3,1}	{-3,5}	⟨v33v_2v_5cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v32⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{3,1,1}	{-3,5}	⟨v23v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩-2⟨v22⟩⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩
{2,2,2}	{-6,10}	⟨v52v32v22cos4ψ_2+6ψ_3-10ψ_5⟩-4⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩^2
{4,1,1}	{3,5}	⟨v24v_3v_5cos8ψ_2-3ψ_3-5ψ_5⟩
{2,3,1}	{9,-5}	⟨v22v33v_5cos4ψ_2-9ψ_3+5ψ_5⟩
{2,2,2}	{0,0}	⟨v22v32v52⟩-⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩^2-⟨v22v32⟩⟨v52⟩-⟨v32⟩⟨v22v52⟩-⟨v22⟩⟨v32v52⟩+2⟨v22⟩⟨v32⟩⟨v52⟩

{m_1,m_2,m_3}	{α_1,α_2}	c2,3,6{α_1,α_2}{m_1,m_2,m_3}
{2,2,1}	{-6,6}	⟨v32v22v_6cos6ψ_3-ψ_6⟩-⟨v22⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{3,2,1}	{0,6}	⟨v23v32v_6cos6ψ_2-ψ_6⟩-⟨v32⟩⟨v23v_6cos6ψ_2-ψ_6⟩
{2,2,2}	{0,0}	⟨v22v32v62⟩-⟨v22v32⟩⟨v62⟩-⟨v32⟩⟨v22v62⟩-⟨v22⟩⟨v32v62⟩+2⟨v22⟩⟨v32⟩⟨v62⟩

{m_1,m_2,m_3}	{α_1,α_2}	c2,4,5{α_1,α_2}{m_1,m_2,m_3}
{1,2,2}	{-8,10}	⟨v52v42v_2cos2ψ_2+4ψ_4-5ψ_5⟩
{2,1,2}	{4,0}	⟨v52v22v_4cos4ψ_2-ψ_4⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v52⟩
{3,1,2}	{-4,10}	⟨v52v23v_4cos6ψ_2+4ψ_4-10ψ_5⟩
{1,3,2}	{12,-10}	⟨v52v43v_2cos2ψ_2-6ψ_4+5ψ_5⟩
{2,2,2}	{0,0}	⟨v22v42v52⟩-⟨v22v42⟩⟨v52⟩-⟨v42⟩⟨v22v52⟩-⟨v22⟩⟨v42v52⟩+2⟨v22⟩⟨v42⟩⟨v52⟩

{m_1,m_2,m_3}	{α_1,α_2}	c2,4,6{α_1,α_2}{m_1,m_2,m_3}
{1,1,1}	{-4,6}	⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{1,2,1}	{8,-6}	⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩
{1,1,3}	{-4,6}	⟨v63v_2v_4cos2ψ_2+2ψ_4-3ψ_6⟩-2⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩⟨v62⟩
{1,3,1}	{-4,6}	⟨v43v_2v_6cos2ψ_2+2ψ_4-3ψ_6⟩-2⟨v42⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{2,1,2}	{4,0}	⟨v62v22v_4cos4ψ_2-ψ_4⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v62⟩
{3,1,1}	{-4,6}	⟨v23v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩-2⟨v22⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{2,2,2}	{-8,12}	⟨v62v42v22cos4ψ_2+2ψ_4-3ψ_6⟩-4⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩^2
{3,2,1}	{0,6}	⟨v23v42v_6cos6ψ_2-ψ_6⟩-3⟨v22v_4cos4ψ_2-ψ_4⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩-⟨v42⟩⟨v23v_6cos6ψ_2-ψ_6⟩
{1,2,3}	{8,-6}	⟨v42v63v_2cos2ψ_2-4ψ_4+3ψ_6⟩-2⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩⟨v62⟩
{1,4,1}	{8,-6}	⟨v44v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩-3⟨v42⟩⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩
{2,2,2}	{0,0}	⟨v22v42v62⟩-⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩^2-⟨v22v42⟩⟨v62⟩-⟨v42⟩⟨v22v62⟩-⟨v22⟩⟨v42v62⟩+2⟨v22⟩⟨v42⟩⟨v62⟩
{3,2,1}	{8,-6}	⟨v23v42v_6cos2ψ_2-4ψ_4+3ψ_6⟩-2⟨v22v_4cos4ψ_2-ψ_4⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩-2⟨v22⟩⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩

{m_1,m_2,m_3}	{α_1,α_2}	c2,5,6{α_1,α_2}{m_1,m_2,m_3}
{1,2,2}	{-10,12}	⟨v62v52v_2cos2ψ_2+5ψ_5-6ψ_6⟩
{2,2,1}	{10,-6}	⟨v52v22v_6cos4ψ_2-10ψ_5+6ψ_6⟩
{3,2,1}	{0,6}	⟨v23v52v_6cos6ψ_2-ψ_6⟩-⟨v52⟩⟨v23v_6cos6ψ_2-ψ_6⟩
{2,2,2}	{0,0}	⟨v22v52v62⟩-⟨v22v52⟩⟨v62⟩-⟨v52⟩⟨v22v62⟩-⟨v22⟩⟨v52v62⟩+2⟨v22⟩⟨v52⟩⟨v62⟩

{m_1,m_2,m_3}	{α_1,α_2}	c3,4,5{α_1,α_2}{m_1,m_2,m_3}
{1,2,1}	{8,-5}	⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩
{2,1,2}	{-4,10}	⟨v52v32v_4cos6ψ_3+4ψ_4-10ψ_5⟩
{3,1,1}	{4,5}	⟨v33v_4v_5cos9ψ_3-4ψ_4-5ψ_5⟩
{1,2,3}	{8,-5}	⟨v42v53v_3cos3ψ_3-8ψ_4+5ψ_5⟩-2⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩⟨v52⟩
{1,4,1}	{8,-5}	⟨v44v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩-3⟨v42⟩⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩
{2,2,2}	{0,0}	⟨v32v42v52⟩-⟨v32v42⟩⟨v52⟩-⟨v42⟩⟨v32v52⟩-⟨v32⟩⟨v42v52⟩+2⟨v32⟩⟨v42⟩⟨v52⟩
{3,2,1}	{8,-5}	⟨v33v42v_5cos3ψ_3-8ψ_4+5ψ_5⟩-2⟨v32⟩⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩

{m_1,m_2,m_3}	{α_1,α_2}	c3,4,6{α_1,α_2}{m_1,m_2,m_3}
{2,2,1}	{0,6}	⟨v42v32v_6cos6ψ_3-ψ_6⟩-⟨v42⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{2,3,1}	{12,-6}	⟨v32v43v_6cos6ψ_3-2ψ_4+ψ_6⟩
{2,2,2}	{0,0}	⟨v32v42v62⟩-⟨v32v42⟩⟨v62⟩-⟨v42⟩⟨v32v62⟩-⟨v32⟩⟨v42v62⟩+2⟨v32⟩⟨v42⟩⟨v62⟩

{m_1,m_2,m_3}	{α_1,α_2}	c3,5,6{α_1,α_2}{m_1,m_2,m_3}
{2,2,1}	{0,6}	⟨v52v32v_6cos6ψ_3-ψ_6⟩-⟨v52⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{1,3,2}	{15,-12}	⟨v62v53v_3cos3ψ_3-5ψ_5+4ψ_6⟩
{2,2,2}	{0,0}	⟨v32v52v62⟩-⟨v32v52⟩⟨v62⟩-⟨v52⟩⟨v32v62⟩-⟨v32⟩⟨v52v62⟩+2⟨v32⟩⟨v52⟩⟨v62⟩

{m_1,m_2,m_3}	{α_1,α_2}	c4,5,6{α_1,α_2}{m_1,m_2,m_3}
{1,2,1}	{10,-6}	⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩
{1,2,3}	{10,-6}	⟨v52v63v_4cos4ψ_4-10ψ_5+6ψ_6⟩-2⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩⟨v62⟩
{1,4,1}	{10,-6}	⟨v54v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩-3⟨v52⟩⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩
{2,2,2}	{0,0}	⟨v42v52v62⟩-⟨v42v52⟩⟨v62⟩-⟨v52⟩⟨v42v62⟩-⟨v42⟩⟨v52v62⟩+2⟨v42⟩⟨v52⟩⟨v62⟩
{3,2,1}	{10,-6}	⟨v43v52v_6cos4ψ_4-10ψ_5+6ψ_6⟩-2⟨v42⟩⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩
{2,2,3}	{-10,18}	⟨v63v52v42cos24ψ_4+5ψ_5-9ψ_6⟩
{3,2,2}	{0,12}	⟨v62v52v43cos12ψ_4-ψ_6⟩-⟨v52⟩⟨v43v62cos12ψ_4-ψ_6⟩
{4,2,1}	{10,6}	⟨v44v52v_6cos28ψ_4-5ψ_5-3ψ_6⟩

{m_1,m_2,m_3,m_4}	{α_1,α_2,α_3}	c2,3,4,6{α_1,α_2,α_3}{m_1,m_2,m_3,m_4}
{1,2,1,1}	{0,-4,6}	⟨v32v_2v_6v_4cos2ψ_2+2ψ_4-3ψ_6⟩-⟨v32⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{1,2,1,2}	{-6,-4,12}	⟨v62v32v_2v_4cos2ψ_2+3ψ_3+2ψ_4-6ψ_6⟩-2⟨v32v_6cos6ψ_3-ψ_6⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{2,2,1,1}	{-6,4,6}	⟨v32v22v_6v_4cos4ψ_2+6ψ_3-4ψ_4-6ψ_6⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{1,2,2,1}	{0,8,-6}	⟨v42v32v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩-⟨v32⟩⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩
{1,2,1,2}	{6,-4,0}	⟨v62v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩-⟨v32v_6cos6ψ_3-ψ_6⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩-⟨v32v_2v_4cos2ψ_2-3ψ_3+2ψ_4⟩⟨v62⟩
{2,2,1,1}	{6,4,-6}	⟨v32v22v_6v_4cos4ψ_2-6ψ_3-4ψ_4+6ψ_6⟩-⟨v22v_4cos4ψ_2-ψ_4⟩⟨v32v_6cos6ψ_3-ψ_6⟩

{m_1,m_2,m_3,m_4}	{α_1,α_2,α_3}	c2,3,5,6{α_1,α_2,α_3}{m_1,m_2,m_3,m_4}
{1,1,1,1}	{3,5,-6}	⟨v_2v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩
{1,1,1,2}	{-3,5,0}	⟨v62v_2v_5v_3cos2ψ_2+3ψ_3-5ψ_5⟩-⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩⟨v62⟩
{2,1,1,1}	{3,-5,6}	⟨v22v_3v_5v_6cos4ψ_2-3ψ_3+5ψ_5-6ψ_6⟩
{2,1,1,2}	{-3,-5,12}	⟨v62v22v_5v_3cos4ψ_2+3ψ_3+5ψ_5-12ψ_6⟩
{1,3,1,1}	{-9,5,6}	⟨v33v_2v_6v_5cos2ψ_2+9ψ_3-5ψ_5-6ψ_6⟩-3⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩⟨v32v_6cos6ψ_3-ψ_6⟩
{1,1,1,3}	{3,5,-6}	⟨v63v_2v_5v_3cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩-2⟨v_2v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩⟨v62⟩
{1,1,3,1}	{3,5,-6}	⟨v53v_2v_6v_3cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩-2⟨v52⟩⟨v_2v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩
{1,2,2,1}	{6,-10,6}	⟨v52v32v_2v_6cos2ψ_2-3ψ_3+5ψ_5-3ψ_6⟩
{1,3,1,1}	{3,5,-6}	⟨v33v_2v_6v_5cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩-⟨v_2v_3v_5cos2ψ_2+3ψ_3-5ψ_5⟩⟨v32v_6cos6ψ_3-ψ_6⟩-2⟨v32⟩⟨v_2v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩
{3,1,1,1}	{3,5,-6}	⟨v23v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩-2⟨v22⟩⟨v_2v_3v_5v_6cos2ψ_2-3ψ_3-5ψ_5+6ψ_6⟩

{m_1,m_2,m_3,m_4}	{α_1,α_2,α_3}	c2,4,5,6{α_1,α_2,α_3}{m_1,m_2,m_3,m_4}
{1,1,2,1}	{-4,0,6}	⟨v52v_2v_6v_4cos2ψ_2+2ψ_4-3ψ_6⟩-⟨v52⟩⟨v_2v_4v_6cos2ψ_2+2ψ_4-3ψ_6⟩
{1,1,2,2}	{4,10,-12}	⟨v62v52v_2v_4cos2ψ_2-2ψ_4-5ψ_5+6ψ_6⟩
{1,2,2,1}	{8,0,-6}	⟨v52v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩-⟨v52⟩⟨v42v_2v_6cos2ψ_2-4ψ_4+3ψ_6⟩
{2,1,2,1}	{-4,10,-6}	⟨v52v22v_6v_4cos4ψ_4-10ψ_5+6ψ_6⟩-⟨v22⟩⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩

{m_1,m_2,m_3,m_4}	{α_1,α_2,α_3}	c3,4,5,6{α_1,α_2,α_3}{m_1,m_2,m_3,m_4}
{1,1,1,1}	{4,5,-6}	⟨v_3v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩
{1,1,1,2}	{-4,-5,12}	⟨v62v_3v_5v_4cos3ψ_3+4ψ_4+5ψ_5-12ψ_6⟩
{1,2,1,1}	{-8,5,6}	⟨v42v_3v_6v_5cos3ψ_3+8ψ_4-5ψ_5-6ψ_6⟩
{1,1,1,3}	{4,5,-6}	⟨v63v_3v_5v_4cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩-2⟨v_3v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩⟨v62⟩
{1,1,3,1}	{4,5,-6}	⟨v53v_3v_6v_4cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩-2⟨v52⟩⟨v_3v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩
{1,2,1,2}	{8,-5,0}	⟨v62v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩-⟨v42v_3v_5cos3ψ_3-8ψ_4+5ψ_5⟩⟨v62⟩
{1,3,1,1}	{4,5,-6}	⟨v43v_3v_6v_5cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩-2⟨v42⟩⟨v_3v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩
{2,1,2,1}	{-4,10,-6}	⟨v52v32v_6v_4cos4ψ_4-10ψ_5+6ψ_6⟩-⟨v32⟩⟨v52v_4v_6cos4ψ_4-10ψ_5+6ψ_6⟩
{3,1,1,1}	{4,5,-6}	⟨v33v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩-2⟨v32⟩⟨v_3v_4v_5v_6cos3ψ_3-4ψ_4-5ψ_5+6ψ_6⟩