An MLE analysis on the relationship between the initial-state granularity and final-state flow factorization

Relativistic hydrodynamics is a widely accepted theoretical framework for modeling the temporal evolution of strongly coupled quark-gluon plasmas produced in relativistic heavy-ion collisions [1, 2, 3, 4, 5, 6, 7]. This macroscopic approach treats the plasma as a continuum, which is crucial for capturing the underlying physics that gives rise to observable phenomena. Key observables in relativistic hydrodynamics include the particle spectrum at intermediate and low transverse momentum, with particular emphasis on collective properties such as flow harmonics and correlations [8, 9, 10, 11, 12, 13, 14]. Experimentally, measurements of azimuthal distributions have been instrumental in demonstrating the concept of a perfect liquid, first observed at RHIC [15]. Consequently, azimuthal anisotropy has emerged as a crucial observable for extracting information about the properties of the underlying physical system [16, 17, 18, 19, 20, 21]. Recently, these observables have been further studied in the context of deformed nuclei [22, 23, 24].

The essence of hydrodynamical evolution can be viewed mainly as the dynamic response to fluctuating initial conditions (IC). Given the inherently nonlinear nature of hydrodynamics, numerous studies have investigated these complex interactions. In particular, significant research has focused on understanding the relationship between initial-state eccentricities and the final-state azimuthal anisotropies [25, 26, 27, 28, 29, 30, 31, 32]. The quantitative decomposition of anisotropic IC was first introduced in Refs. [25, 26]. This approach relies on a cumulant expansion, where each expansion coefficient captures the “connected” eccentricity component at a given order. Consequently, higher-order cumulants are defined by subtracting contributions from the “disconnected” combinations of lower-order terms. Furthermore, flow harmonics represent the hydrodynamic response to IC fluctuations, classified according to their cumulant order, with the lowest-order cumulants generally assumed to have the most significant impact.

In the literature, contributions proportional to cumulants of the same azimuthal order are often identified as linear responses. In contrast, contributions from combinations of lower-order cumulants that result in the same azimuthal order are attributed to non-linearity. In practice, the response strength varies across different cumulant combinations. Specifically, the established mapping between IC and flow harmonics is defined by the strongest correlation between an optimized linear combination of cumulant products and the corresponding flow harmonic [27, 30]. Numerical studies confirm that this mapping is particularly effective for more central collisions, inspiring further research in this area [28, 31, 32]. For example, one might examine whether each individual component of a given azimuthal order results in a linear response by isolating them from the IC [29]. It should be noted that each azimuthal harmonic order corresponds to an infinite set of moments or cumulants. Therefore, it remains unclear whether a specific one-to-one mapping exists or to what extent different terms of the same azimuthal order mix during dynamical evolution. Consequently, this ambiguity leaves room for alternative descriptions of the radial expansion.

In this regard, some authors argued that the cumulant expansion may not be optimal for decomposing the IC. Alternatively, the Bessel-Fourier expansion proposed in Refs. [33, 34] utilizes an orthonormal basis to decompose IC fluctuations. A vital advantage of this approach is that it orders fluctuations based on their wavelength, allowing shorter wavelength radial modes to be effectively suppressed. This enables a more precise separation of different modes within a hydrodynamic framework [35]. Another intuitive method for capturing the granularity of the IC is the peripheral tube model [36, 37, 38, 39], which provides an intuitive interpretation for the generation of triangular flow and ridge structures observed in di-hadron correlations. Unlike Fourier-based eccentricity decompositions, this model aims to capture the localized feature in the event-by-event fluctuating IC. It is motivated by heuristic arguments and numerical simulations, which suggest that high-energy-density peaks are formed near the surface region due to elementary binary collisions in the transverse plane. These localized regions naturally form a tube-like structure along the longitudinal direction and are shown to correlate strongly with the observed ridge-like structures in two-particle correlations. This model emphasizes the impact of localized IC fluctuations, as opposed to global sinusoidal expansions based on the azimuthal moment expansion. For instance, a tube located deep inside the medium, though it has a sizable contribution to the corresponding moments, would have its hydrodynamic effect significantly suppressed by the surrounding matter. In contrast, a tube near the surface can induce notable distortions in the single-particle azimuthal distribution and influence two-particle correlations. This model has been successfully employed to study various features in di-hadron correlations, showing good agreement with experimental data [36, 40, 41, 39].

The anisotropic distribution of final-state particles in momentum space is characterized by flow harmonics $v_{n}$, defined through the one-particle distribution function [9, 42]:

where $\phi$ is the azimuthal angle of the emitted particle and $\Psi_{n}$ represents the event plane for the harmonic order $n$. Elliptic flow ($v_{2}$) and triangular flow ($v_{3}$) are particularly relevant observables, with $v_{2}$ primarily arising from the almond-shaped geometry of the overlap region [8] and $v_{3}$ stemming from event-by-event IC fluctuations [43]. Considerable research has focused on exploring the relationship between initial geometric anisotropies and final-state flow harmonics, primarily to investigate non-linear effects [44, 45, 46, 47, 48, 49, 50], eccentricity and flow fluctuations [51, 52, 53, 54], and multi-particle correlations [55, 56].

Various techniques have been developed to determine flow harmonics $v_{n}$ from experimental data. The traditional event plane method [9, 57] estimates the event plane angles $\Psi_{n}$ to evaluate the harmonics in Eq. (1), reflecting the fact that the reaction plane[43] cannot be directly measured. Other methods, such as particle correlations, use Q-vectors and cumulants [58, 11, 59, 60] to eliminate the need for $\Psi_{n}$. This approach allows cumulants to be compactly expressed using generating functions [11, 61]. It includes techniques like particle cumulants [11, 61, 59], Lee-Yang zeros [62, 63], and symmetric cumulants [64]. Recently, an alternative method based on the maximum likelihood estimator (MLE) was proposed [65], treating the flow coefficients as unknown parameters of a hypothetical probability distribution derived from experimental data. The MLE approach offers several advantages. First, it effectively handles nonflow effects, which refer to correlations that cannot be attributed to collective behavior as per Eq.(1). Momentum conservation, for example, can cause deviations from a purely flow-driven spectrum [66]. Second, MLE is efficient at the limit of large samples and can incorporate additional constraints in the flow distribution function, making it robust in scenarios with nonflow effects. Furthermore, MLE’s asymptotic normality ensures unbiased or nearly unbiased results, making it well-suited for the extensive datasets obtained at RHIC and LHC. Given these attributes, MLE provides a unique tool for analyzing multi-particle correlations and assessing the connection between initial anisotropy and final-state flow harmonics.

Beyond flow harmonics, the discussion can be extended to more complex observables constructed using generic multi-particle correlators. As highlighted in the literature [67, 68], such quantities are often more sensitive to the underlying IC’s specific characteristics than their global averaged properties. Specifically, event-by-event fluctuations in the IC can introduce significant effects in multi-particle correlators, leading to transverse momentum ($p_{\rm T}$) dependent event planes [53]. To capture these variations, the final-state event planes are typically estimated using the azimuthal distribution of particles over a broad $p_{\rm T}$ range on an event-by-event basis. Experimental data [69, 70] and hydrodynamic simulations [43, 53] indicate that the event planes fluctuate significantly across different $p_{\rm T}$ intervals. This observation, which is supported by studies such as Refs. [71, 72, 73], reveals that the correlation matrix in transverse momentum approximately factorizes. Such factorization is considered compelling evidence for the hydrodynamic picture of heavy-ion collisions. Consequently, flow factorization [67, 68, 53] has emerged as a powerful tool for probing the properties of initial-state fluctuations.

The breakdown of flow factorization is primarily attributed to event-by-event fluctuations in the initial energy distribution [67, 68, 53, 71, 72] rather than to variations in the transport properties of the medium. This suggests that a detailed analysis of factorization breakdown can yield valuable insights into the IC of the system. Moreover, distinct differences have been observed between flow Pearson correlations constructed using different moments [74]. Given the ability of the MLE to preserve the structure of these correlations owing to its equivariance properties, it provides a robust and meaningful framework for assessing multi-particle correlators and their sensitivity to initial-state granularity.

The present study is motivated by the above considerations. We aim to explore the relationship between initial-state anisotropies and final-state flow harmonics and their correlations by employing MLE as a statistical estimator for the flows and other higher-order moments. Specifically, we use the peripheral tube model to quantify IC granularity and employ flow factorization to capture nuances in high-order flow harmonics and multi-particle correlators. In particular, MLE has been applied to specific correlators that are otherwise inaccessible for an arbitrary combination.

The remainder of this paper is organized as follows. In the next section, we discuss MLE as a statistical estimator for multi-particle correlators and flow factorization. In Sec. III, we introduce the peripheral tube model and describe its application in modeling ICs with varying granularity. Numerical studies are carried out in Sec. IV, where we investigate the relationship between IC granularity and final-state flow factorization. The results are compared with those obtained using cumulants and the event-plane method. The final section provides further discussions and concluding remarks.

The most prominent methods for extracting flow harmonics are based on multi-particle correlations [57]. These techniques rely on the following definition for $k$-particle correlator [52]:

where $\langle\cdots\rangle$ denotes the average over distinct tuples of particles, assuming independent particle emission as described by Eq. (1) in the limit of infinite multiplicity.

To isolate $v_{n}$, one typically chooses a specific set of indices $(n_{1},\cdots,n_{k})$ such that [52]

which ensures that all contributions from the event planes cancel in the exponential. This condition effectively reduces the expression to a formalism independent of the event-plane angles [55]. The simplest example is the two-particle correlation ($k=2$), where $n_{1}=-n_{2}=n$. This choice yields

which directly relates the two-particle correlation to the square of the flow harmonic $v_{n}$.

In realistic scenarios, however, the number of particles ($M$) in a given event is finite. Thus, the analysis is performed using discrete azimuthal angles $\phi_{1},\phi_{2},\cdots,\phi_{M}$ corresponding to the measured particles. Rather than employing the integration in Eq. (4), which is not viable for a finite number of particles $M$, it is intuitive to use the following summation [64]:

to estimate the flow harmonic $v_{n}$.

As a more sensitive observable, flow factorization is defined as a Pearson correlation in terms of flow vectors evaluated at different transverse momenta, namely,

where $p_{\rm T}^{\rm t}$ and $p_{\rm T}^{\rm a}$ are transverse momenta of the trigger and associated particles, and $V_{n\Delta}(p_{\rm T1},p_{\rm T2})$ is the $n$th harmonic of the underlying di-hadron azimuthal distribution with transverse momenta $p_{\rm T1}$ and $p_{\rm T2}$, namely,

where

is known as the flow vector [75, 76]. Because of its explicit consideration of transverse momentum dependence, Eq. (7) can be viewed as the differential counterpart of the two-particle correlation defined in Eq. (4). Similar to Eq (5), in practice, Eq. (7) is implemented by enumerating all possible combinations of the measured particles, as given in the Appx. A of [74]

From the statistical inference perspective, Eq. (5) serves as an estimator, providing an estimation for $v_{n}^{2}$ rather than $v_{n}$ itself, based on a finite number of observations. The first two moments of this estimator can be readily evaluated [65] and generally do not vanish. Specifically, this estimator has a finite variance that decreases with increasing multiplicity. For higher-order correlators, one can employ generating functions [11, 77, 78] or compute them directly [59] using $Q$-vectors [58] or flow vectors [75, 76]. Notably, quantities constructed using $Q$-vectors or flow vectors can also be interpreted as estimators, serving as generalized extensions of the unweighted sums in Eq. (2). Notably, the variance of these quantities largely remains finite [79], indicating an inherent statistical uncertainty due to the finite number of measured events and particle multiplicity. A finite variance implies that uncertainties in the estimated flow harmonics are inevitably influenced by limited statistics, particularly for a finite number of events. Therefore, such limitations of statistical origin must be considered when comparing results obtained from different flow measurement methods.

Following this line of thought, MLE could serve as an alternative estimator for flow and related observables, a possibility explored in a previous study [65]. Expressly, for a given set of observations $y\equiv(y_{1},y_{2},\cdots,y_{M})$, we assume that they are sampled from a joint probability distribution governed by several unknown parameters $\theta\equiv(\theta_{1},\theta_{2},\cdots,\theta_{m})$. As mentioned in the Introduction, the likelihood function $\mathcal{L}$ for the observed data is given by:

which represents the joint probability density for the given observations evaluated at the parameters $\theta$. The goal of MLE is to determine the parameters for which the observed data attains the highest joint probability, namely:

where $\Theta$ is the domain of the parameters. In particular, for independent and identically distributed (i.i.d.) random variables, $f(y;\theta)$ is given by a product of likelihood functions $f^{\mathrm{uni}}$:

This scheme can be readily applied to collective flow in heavy-ion collisions. Considering an event consisting of $M$ particles, the likelihood function reads:

where the likelihood function $\mathcal{L}$ is governed by the one-particle distribution function Eq. (1). The last equality is based on the assumption that the particles’ azimuthal angles are i.i.d. variables. The parameters of the distribution, $\theta=(v_{1},v_{2},\cdots,\Psi_{1},\Psi_{2},\cdots)$, are the flow harmonics and the event planes.

In practice, one often chooses the objective function to be the log-likelihood function ${\ell}$:

Numerical calculations indicate that Eq. (13) is more favorable than Eq. (12), although as multiplicity $M$ increases, an appropriate implementation should be adopted to avoid the increasing truncation error.

The maximum of $\ell$ occurs at the same value of $\theta$, which maximizes $\mathcal{L}$. For $\ell$ that is differentiable in its domain $\Theta$, the necessary conditions for the occurrence of a maximum are:

As discussed in the Introduction, MLE has asymptotic normality, which attains the Cramér-Rao lower bound as the sample size increases. In other words, no consistent estimator has a lower asymptotic mean squared error than the MLE. In the context of relativistic heavy-ion collisions, all events of a given multiplicity asymptotically form a (multivariate) normal distribution:

where $\theta_{0}$ represents the true value, and $I_{M}(\theta)$ is the Fisher information matrix, defined as:

where the expectation $E_{\theta}$ is taken with respect to the distribution $f(\phi_{1},\cdots,\phi_{M};\theta)$. For i.i.d. data, the Fisher information possesses the form:

where $I_{1}$ is the Fisher information matrix for a single observation. As a result, the standard deviation of MLE is expected to be roughly proportional to $\frac{1}{\sqrt{M}}$. As discussed in [65], these properties can be further quantified using the Wald, likelihood ratio, and score tests.

The peripheral tube model [39] offers a simplified yet effective framework for understanding the generation of triangular flow and the resulting particle correlations. This approach is intrinsically connected to event-by-event fluctuating hydrodynamics, where the IC are represented by a few high-energy tubes near the system’s surface. Each of these high-energy tubes independently influences the local hydrodynamic evolution, and these tubes’ collective contributions combine to produce the observed particle correlations.

The present study makes use of the tube model to quantify the granularity of the fluctuating IC, while the bulk of the hot matter is substituted by an average energy density distribution obtained from multiple events of the same centrality class. We systematically vary the size and number of these tubes to investigate their impact on flow harmonics and particle correlations, focusing on flow factorization. The IC in this model is devised to reflect the underlying event-by-event fluctuating initial energy density distribution, consisting of two main components: a smooth background and several high-energy-density tubes close to the surface. The smooth background captures the bulk properties of the system, while the tubes represent localized fluctuations on an event-by-event basis. The energy density profile in the model is expressed as:

where the average background distribution $\epsilon_{\mathrm{bgd}}$ is given by:

where the radial coordinate is defined as:

where $K,~{}L,~{}M,~{}a,~{}b,c$ are parameters obtained by numerical calibration to the average ICs generated by the NeXuS and EPOS models [80, 81, 82, 83, 84].

The profile of an individual high-energy tube is given by:

where $A_{\mathrm{tube}}$ and $R_{\mathrm{tube}}$ denote the maximum energy density and radius of the tube, respectively. The radial position $r_{\mathrm{tube}}$ is defined as:

with the spatial coordinates $(x_{\mathrm{tube}},y_{\mathrm{tube}})$ given by:

where $r_{0}$, $a$, $b$, and $\theta$ determine the radial location and orientation of the tube. For the present study, the number of tubes will be varied to reflect the different degrees of granularity of the IC. Previous analysis suggests that this parameter has a limited effect on the resulting flow harmonics and two-particle correlations [85]. Nonetheless, as shown below, different granularity has a rather subtle effect on the resulting collect flows. The parameters used in this model are calibrated using typical IC profiles for Au+Au collisions in the $0\%-5\%$ centrality class at $\sqrt{s_{NN}}=200$ GeV. For a randomly generated event, the radii and azimuthal angles of the tubes are drawn from the following uniform distributions $r_{0}\sim\mathrm{U}(0,0.546)$ and $\theta\sim\mathrm{U}(0,2\pi)$. The remaining parameters are summarized in Tab. 1.

Using the devised IC furnished by the peripheral tube mode, the final-state particles are obtained by numerical simulations using the hydrodynamic code NeXSPheRIO [86]. We evaluate the temporal evolution, differential flow harmonics, and multi-particle correlations in terms of flow factorization. The calculations are performed for IC configurations with different numbers of tube $N_{\mathrm{tube}}$. The numerical results are shown in Figs. 1-4.

As shown in Fig. 1, the dynamic evolutions can drastically differ depending on different IC configurations. As demonstrated in the left and middle columns of Fig. 1, when there are only a few tubes, the fluid is observed to be deflected in the vicinity of each tube to both sides. In particular, the evolution around an individual high-energy tube leads to two peaks separated by roughly $120$ degrees in the azimuthal distribution, as indicated by the left-most column of Fig. 1. Subsequently, this gives rise to the desired two-particle distribution where a double peak is formed on the away side, as pointed out in previous studies [40, 38]. However, as the number of tubes increases, the hydrodynamic evolution associated with different tubes becomes significantly interfered. As shown in the right column, the resulting evolutions are rather complex.

However, the rather drastic difference in the initial conditions and hydrodynamic evolutions are not directly reflected in the flow harmonics, which measure the global anisotropies in the momentum space. As indicated by Fig. 2, the resulting differential elliptic and triangular flow harmonics are mainly irrelevant to the difference in the IC. Specifically, the resultant differential flows are rather robust against the variation of the number of tubes. This result can be understood as the background primarily governs the overall elliptic and triangular shape of the IC, while the number of tubes impacts mainly the granularity. Consistent with previous findings [85], this result implies that the resulting two-particle correlation structure remains, by and large, unchanged. Notably, the results shown in Fig. 2 are carried out using two approaches, the MLE and the multi-particle cumulant. On the one hand, the two-particle cumulant estimates flow harmonics through $v_{n}^{2}$ by Eq. (5). By definition, the estimation is skewed owing to the presence of flow fluctuation [87, 79]. On the other hand, the MLE evaluates the flow harmonics by maximizing the likelihood, Eq. (10), which also might lead to a biased estimation before attaining the statistical limit. Nonetheless, it is observed that these two methods give consistent results as good agreements are manifestly reached for the differential flow. Therefore, one must explore more sensitive quantities to scrutinize the observational impact from different granularities of the fluctuating IC.

In this regard, we proceed to evaluate flow factorization [67, 68], which is more susceptible to initial state fluctuations. For flow factorization, we elaborate on three different scenarios, and the results are presented in Figs. 3 and 4. The numerical results indeed indicate that a significant difference is observed in such quantities. Specifically, we focus on two types of deviations. First, we explore the flow factorization’s dependence on the event’s granularity by varying the number of tubes that constitute the IC. Secondly, a sizable difference is also observed across different flow estimators, namely, the MLE and multi-particle cumulants. In other words, even though flow harmonics are found to be broadly consistent across different estimations, observables associated with the high-order moment of the one-particle distribution Eq. (1) entail rather substantial deviations.

Specifically, we consider three types of factorization ratios. The first scenario involves a ratio regarding two identical flow harmonics, as defined in Eq. (6). This quantity is evaluated for elliptic and triangular flows by employing IC with different numbers of tube $N_{\mathrm{tube}}$. As presented in the left column of Fig. 3, the calculated factorization ratio $r_{n}$ is shown as a function of the difference between the transverse momenta of the trigger and associated particles. The calculations are carried out using the MLE method. At the origin, the two transverse momentum intervals of $p_{\rm T}^{\rm a}$ and $p_{\rm T}^{\rm t}$ coincide, and therefore, any substantial deviation from the unit comes solely from the correlation within the small given interval. Our numerical calculations indicate that the deviation from perfect factorization mostly vanishes at the origin of the coordinates, consistent with the experimental data [71, 88]. In theory, this is understood because at the limit, when the size of the interval vanishes, the Pearson correlation falls back to that between two identical quantities, which is guaranteed to have a perfect correlation. Moreover, based on Eq. (8), the flow vector can be further factorized if there were no fluctuations. Even though such a factorization is not exact, it indicates that $r_{n}$ can be roughly viewed as receiving contributions from two factors: the Pearson correlation of flow harmonics and event-plane correlation. These results are shown in the middle and right columns of Fig. 3. By comparing the three columns, it is observed that correlations between flow harmonics and event planes are both crucial in forming the resulting flow factorization. The results obtained for different numbers of tubes indicated that a more significant breakdown of the factorization occurs when the number of tubes is small and gradually recovers when the number of tubes increases. For a large number of tubes, the resulting IC is somewhat reminiscent of that of realistic IC generated by NeXuS [80, 81]. This is particularly true for the elliptic flow, as the resulting factorization ratios approach that of Au+Au collisions, as indicated by black filled squares in Fig. 3. This result demonstrates that flow factorization is a more sensitive quantity to quantify the initial state fluctuations, precisely, the degree of granularity.

The second scenario involves a mix of three different flow harmonics of the form

where the three indices satisfy the relation

Such a quantity has been investigated in the literature regarding flow fluctuations [89, 90], as the term related to event-plane correlation readily vanishes if one has a global constant event plane. As the above quantity depends on three transverse momenta, in practice, one assigns two particles as triggers and one remaining particle as an associated particle, whose transverse momentum will be integrated over a given interval, as shown in the figures. Here, we are also interested to compare the resulting factorization ratios using two different approaches: the MLE and particle cumulant methods. The resulting flow factorization ratios are shown in Fig. 4, where one considers the range $2.0<p_{\mathrm{T}}^{a}<2.5$ GeV.

One observes that the factorization ratios obtained using two different methods possess a similar trend. Nonetheless, the difference between these two methods is rather pronounced. In addition to the analysis performed in Fig. 3, this further reinforces our conclusion that factorization is a more sensitive quantity to the initial state fluctuations. Besides that it distinguishes different degrees of granularity, the quantity is rather sensitive to the employed method, particularly when the quantity in question involves an explicit role of higher order moments. In other words, the consistency between different approaches observed for differential flow in Eq. (2) does not generalize straightforwardly to factorization ratios. Also, for such a scenario, it is noted that the factorization ratio does not traverse the origin of the coordinates. This is understood as, in this case, the underlying correlation does not assume the Pearson correlation of the same quantity as its limit, and therefore the factorization ratio does not fall back to unity even when the momentum intervals coincide.

Lastly, we are interested in a further variation of the factorization ratio, essentially when the relation between indices Eq. (24) is no longer satisfied. For this case, even if the event plane is a global constant, it will not be canceled out by the particle tuple combination in question. As a result, such a case corresponds to a specific scenario where the MLE plays a more distinct role. In the left panel of Fig. 5, we explore the dependence of the results on IC of different granularities furnished by different numbers of tubes. A feature similar to Fig. 3 is observed, as the degree of factorization breaking decreases as the number of tubes increases. Nonetheless, we note that it would still be possible to form particle pairs to assess the harmonics through the higher moments of the underlying one-particle distribution. However, in such a way, we argue that the results might be quite sensitive to the specific construction. This can be shown by explicitly evaluating their MLE counterparts in terms of high moments. The calculations can be carried out straightforwardly owing to the equivariance of MLE [91]. The results are presented in the right panel of Fig. 5, for which the calculations are carried out for the IC generated using $N_{\mathrm{tube}}=100$ randomized tubes, corresponding to the black filled circles shown in the left panel. For both panels, again, the factorization ratio does not traverse the origin of the coordinates. It is observed that the use of different moments affects not only the magnitude of the factorization ratio but also its momentum dependence. Even for the same construction, the results obtained by using the MLE and multi-particle cumulant methods are different. In particular, as the MLE and particle correlators are not mathematically equivalent statistical estimators, such a difference is not unsurprising. We, therefore, argue that the MLE provides a meaningful alternative to assess collective flows besides existing means.

Using the MLE, we investigated the relationship between initial-state fluctuations and final-state flow anisotropies in relativistic heavy-ion collisions. By reflecting on existing results that the mostly linear relation between initial state eccentricities and final state anisotropies [28, 31], the present study proceeds further based on the observation that distinct IC may produce almost identical momentum-space anisotropy in terms of flow harmonics [32]. In this regard, our primary focus was on evaluating how the granularity of IC, modeled using a peripheral tube approach, impacts flow factorization, a more sensitive probe compared to traditional flow harmonics. Specifically, while differential flow harmonics, such as $v_{2}$ and $v_{3}$, showed minimal sensitivity to changes in the number and configuration of tubes, the flow factorization ratio displayed substantial variation, highlighting its potential as an effective observable for uncovering detailed information about the initial state.

By employing MLE, we extracted some specific mix-order correlators that are otherwise challenging to access, and the results showed distinct differences compared to standard techniques. The present study further generalizes our initial proposal [65, 74] to a more specific subject. A primary challenge of the approach was its computational cost. Our findings indicate that such an extension is computationally feasible. Moreover, MLE provides a more nuanced understanding of flow factorization and its dependence on initial-state granularity. As an asymptotically normal estimator, MLE’s robustness and flexibility offer a new approach for analyzing complex multi-particle correlations. In particular, its implementation does not rely on constructing particle pairs or tuples to cancel event planes, substantially leading to broader applications. Also, compared to the standard methods, it is more flexible to deal with scenarios where a template is not known prior to the analysis. This feature indicates a promising aspect in flow analysis that MLE can be applied. We plan to continue exploring related topics in future studies.

We are thankful for enlightened discussions with Sandra Padula, Takeshi Kodama, and Yogiro Hama. We gratefully acknowledge the financial support from Brazilian agencies Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). A part of this work was developed under the project Institutos Nacionais de Ciências e Tecnologia - Física Nuclear e Aplicações (INCT/FNA) Proc. No. 464898/2014-5. This research is also supported by the Center for Scientific Computing (NCC/GridUNESP) of São Paulo State University (UNESP). SFS acknowledges the funding from the Physics Master Teacher Studio of Guangsha University. CY acknowledges the support of the Postgraduate Research & Practice Innovation Program of Jiangsu Province under Grant No. KYCX22-3453. JL acknowledges the support of the National Natural Science Foundation of China under Grant No. 12347101.