Hydrodynamical noise and Gubser flow

Hydrodynamical noise is introduced on top of Gubser's analytical solution to viscous hydrodynamics. With respect to the ultra-central collision events of Pb-Pb, p-Pb and p-p at the LHC energies, we solve the evolution of noisy fluid systems and calculate the radial flow velocity correlations. We show that the absolute amplitude of the hydrodynamical noise is determined by the multiplicity of the collision event. The evolution of azimuthal anisotropies, which is related to the generation of harmonic flow, receives finite enhancements from hydrodynamical noise. Although it is strongest in the p-p systems, the effect of hydrodynamical noise on flow harmonics is found to be negligible, especially in the ultra-central Pb-Pb collisions. For the short-range correlations, hydrodynamical noise contributes to the formation of a near-side peak on top of the correlation structure originated from initial state fluctuations. The shape of the peak is affected by the strength of hydrodynamical noise, whose height and width grow from the Pb-Pb system to the p-Pb and p-p systems.


I. INTRODUCTION
One of the recent focuses on relativistic heavy-ion collisions carried out at the Large Hadron Collider (LHC) and Relativistic Heavy-Ion Collider (RHIC) lies in small colliding systems, including p-p and p-Pb at the LHC [1][2][3][4][5][6], d-Au and 3 He-Au at RHIC [7,8]. Especially in the collision events with sufficiently high multiplicity, it has been noticed that observables related to multi-particle correlations are consistent with a medium collective expansion scenario [6,9,10]. Relativistic hydrodynamics is a natural candidate to simulate and investigate the collective expansion of a QCD medium, even in small colliding systems [11][12][13][14][15][16][17]. However, applying viscous hydrodynamics to small colliding systems is challenged by a couple of factors. First, the applicability of viscous hydrodynamics is constrained by the convergence of gradient expansion, which is normally quantified by Knudsen number Kn. In smaller systems, the smaller system size leads to a larger value of Kn and the applicability of hydrodynamics is deteriorated accordingly [18]. Second, although hydrodynamical noise was introduced long ago by Landau and Lifshitz [19,20], and was recently generalized to relativistic systems by Kapusta, Müller and Stephanov [21], it is neglected in most of the present hydrodynamic simulations. Since hydrodynamical noise is generically associated with dissipations and expected to be more pronounced in small systems, its influence in the small colliding systems needs to be clarified. The purpose of this work is to quantitatively estimate the effects of thermal noise in Pb-Pb, p-Pb and p-p collisions at the LHC energies, in order to test the applicability of noisy and viscous hydrodynamics in these systems.
Instead of more explicit numerical simulations of hydrodynamics with thermal noise [22,23], this work resorts to an analytical solution of viscous hydrodynamics given by Gubser and Yarom [24,25], known as the Gubser flow. By doing so, the inclusion of thermal noise in hydrodynamics is simplified dramatically both theoretically and numerically. Although Gubser flow is not as realistic as being required for heavy-ion collisions, it mimics the expansion of a hot conformal medium after its thermalization. Therefore, as one preliminary work on hydrodynamical noise, we will restrict ourselves to the initial stage of heavy-ion collisions. We will not address hadronization and freeze-out. This paper is organized as follows. The theoretical framework of viscous hydrodynamics and hydrodynamical noise is discussed on a general ground in Section II. In Section III we briefly review Gubser flow, then thermal noise is introduced into Gubser's solution in Section IV in parallel to the case of 1+1D Bjorken flow [21]. In order to estimate the effects of hydrodynamical noise with respect to the Pb-Pb, p-Pb and p-p systems, we solve the noisy Gubser flow in Section V, with discussions in terms of formal solutions presented in Section V A and numerical simulations in Section V B. Summary and conclusions are given in Section VI.

II. HYDRODYNAMICS AND HYDRODYNAMICAL FLUCTUATIONS
In this section and throughout this work, we work with a metric signature that is mostly positive (−, +, +, +), thereby the flow velocity u µ is normalized as u 2 = −1, and the projection operator ∆ µν is defined as ∆ µν = u µ u ν +g µν . Tensor indices within angular brackets are transverse, traceless and symmetric, while tensor indices inside parentheses are symmetric. Except being specified as in Section V A, we also use angular brackets around a quantity to denote ensemble average, which is defined as the average over simulation events with respect to the same initial condition.
The evolution of quark-gluon plasma in heavy-ion collisions can be well described by relativistic hydrodynamics, together with an equation of state (EOS) originated from lattice QCD. Neglecting baryon number conservation, hydrodynamics regarding heavy-ion collisions is formulated as the conservation of energy-momentum, where we have taken d µ to indicate the covariant derivative and T µν is the energy-momentum tensor. In the Landau frame, with dissipative effects characterized by the stress tensor Π µν . Viscous hydrodynamics generally takes a form of gradient expansion, where σ µν = 2∇ µ u ν and ∇ µ ≡ ∆ µν d ν . Π µν has a determined form up to first order in the expansion, which is known as the Navier-Stokes hydrodynamics.  [19,20] for the first order dissipative hydrodynamics, and extended to the framework of relativistic hydrodynamics by Kapusta, Müller and Stephanov [21]. If one associates for each thermodynamical quantity with a fluctuation term which characterizes thermal noise, e.g., temperature, energy density and pressure are expressed as (suffix '0' indicates ensemble-averaged quantity.) and so for each dynamical quantity in the fluid system, Note that the hydrodynamical noise term S µν is introduced with respect to Π µν . One thus arrives at a new expression for the energy-momentum tensor, T µν 0 in Eq. (2.6) indicates contributions from ensemble-averaged hydro quantities, while δT µν is determined by thermal fluctuations. Apart from cases, such as those where phase transition plays an significant role in the system evolution (cf. [26]), which is beyond the scope of this work, thermal fluctuations are relatively small variables. Accordingly, one can treat thermal fluctuations perturbatively with respect to the background evolution d µ T µν 0 = 0. To the first order in fluctuations, hydro equations of motion for fluctuations d µ δT µν = 0 are linearized and lead to where w = +P is the enthalpy density. δΠ µν in Eqs. (2.7) is a term induced by δT and δu µ . Without the hydrodynamical noise term S µν , Eqs. (2.7) also describe the hydro evolution of initial state fluctuations to linear order. In Eqs. (2.7) and in the following, we ignore the suffix '0' for the ensemble-averaged quantities for simplicity, when confusions do not arise.
When the system is in thermal equilibrium, two-point autocorrelations of these fluctuations in Eqs. (2.4) and δu µ are known to be local in space and time, with correlation strength constrained by thermodynamical variables in equilibrium [19,20]. Autocorrelations of thermal quantities are related with each other through the equation of state. When the system is out-of-equilibrium and driven by hydrodynamics, especially when dissipation is taken into account, the autocorrelation of these fluctuations must be determined with respect to that of hydrodynamical noise S µν , according to the fluctuation-dissipation theorem, where the Onsager coefficient γ µναβ relates the symmetric tensor ∇ (µ u ν) /T to Π µν , In addition to the symmetry between pairs of indices (µν) and (αβ), which is guaranteed by Onsager's relation [20], γ µναβ is symmetric, traceless (with conformal EOS) and transverse in µ and ν, as well as in α and β. It can be shown that for Navier-Stokes hydrodynamics with the tensor structure ∆ µναβ defined by projection operators, The two-point auto-correlation determined by Eqs. (2.8) and (2.9) are characteristic only for the first order dissipative hydrodynamics, which represents white noise with correlation strength constrained by first order transport coefficients. When higher order viscous corrections are taken into account with respect to the causality, space-time dependence must be altered from a Dirac delta function and one accordingly obtains colored noise which depends also on higher order transport coefficients [27]. Although Navier-Stokes hydrodynamics suffers from causality problem, which contains superluminal modes corresponding to sufficiently large momentum. In this work, we shall only solve Navier-Stokes hydrodynamics with finite number of hydrodynamic modes, in a way that the issue of causality becomes less significant. where ∆ µναβ ζ = 1 2 ∆ µν ∆ αβ is factorized automatically in terms of the projection operator.

III. GUBSER FLOW
In this section we briefly review some of the essential ingredients of Gubser flow, which are relevant to our present work. More details of Gubser's solution to relativistic hydrodynamics and discussions can be found in [24,25].
Analytical solutions to viscous hydrodynamics can be achieved with respect to certain symmetry constraints, which allows one to recover the solution in a more general coordinate system by coordinate transformations. For instance, the well-known Bjorken flow relies on a Bjorken boost invariance regarding the space-time geometry in heavy-ion collisions, which is explicit if one writes the metric tensor in Milne coordinates (beam axis along z), and one immediately finds the solution of flow velocity profile with u ξ = 0, which then gives rise to v z = z/t in the original space-time. Recently, Gubser and Yarom developed a new type of analytical solution by imposing rotational symmetry with respect to the beam axis [24,25], in addition to Bjorken boost invariance. Starting from the Milne space-time Eq. (3.1), the manifold R 3,1 is firstly transformed into dS 3 × R via a Weyl rescaling. The following mapping, further allows one to rewrite the metric tensor as dŝ 2 = −dρ 2 + dξ 2 + cosh 2 ρ dθ 2 + sin 2 θdφ , We follow the same notation in [25] that a 'hat' over a quantity indicates that the quantity is defined in the coordinate system Eq. (3.4). Although flow velocity profile in Eq. (3.5) describes a fluid cell at rest, there exists non-zero expansion rate due to the geometry, which results in a diagonal shear tensor σ µν , whereP µν is a traceless projection operator which is identical to∆ µν except for the ξξ component:P ξξ = −2∆ ξξ . Simplification of hydro equations of motion is expected accordingly. And indeed the only non-trivial equation left from Eq. (2.1) is an equation of continuity, which still requires solution with respect to the gradient expansion order by order. The analytical solution to Navier-Stokes hydrodynamics can be found for a conformal equation of state:ˆ = 3P. For instance, the energy density in the 'hat' system is [25] (ρ) = (cosh ρ) − 8 where F (ρ) is a function defined by the following integral and H 0 is a constant proportional to η/s which parameterizes shear viscosity in the 'hat' coordinate system.T 0 is a dimensionless parameter constrained by the total multiplicity of the collision event, whose details will be given later in Section V regarding numerical simulations.
Once hydrodynamics is fully solved in the 'hat' coordinate system, quantities in the (τ, x, y, ξ) system can be recovered via the following transformations, It has been noticed that the gradient expansion breaks down in Gubser's solution at early de Sitter times for Navier-Stokes hydrodynamics [25], and also at late de Sitter times if one investigates second order viscous hydrodynamics. Such constraints indicate initialization times at which the thermal system can be approximated by viscous hydrodynamics. Accordingly, regarding realistic systems in heavy-ion collisions, we shall avoid very early or late de Sitter times in our analysis.

IV. NOISY GUBSER FLOW
To a large extent, effects of hydrodynamical noise can be investigated analytically with respect to solvable hydro models, such as Bjorken flow in 1+1D and Gubser flow. The discussion on hydrodynamical noise with respect to Bjorken flow in 1+1D was given previously in [21], from which we generalize to the case of Gubser flow.

A. Bjorken flow
Following the discussions in [21], for a 1+1D Bjorken's solution with respect to the Navier-Stokes hydrodynamics, the essential simplification is a factorization of the tensor structure of the two-point autocorrelation function in Eq. (2.11) and one thus expects where P µν is the same tensor defined under Eq. (3.7), but now in the Milne space-time.
Notice that we write the factorization in Eq. (4.2) in terms of the traceless projector P µν instead of ∆ µν [21], which differs by a factor of two in this particular case of Bjorken flow in 1+1D. 2 Generally, for a conformal fluid system, P µν is preferred since γ µναβ is traceless in µ and ν. With dimension being saturated by the enthalpy density w(τ ) in Eq. (4.2) and tensor structure being represented by P µν , hydrodynamical noise associated with shear viscous tensor is reduced to one unknown scalar function f (τ, ξ), which satisfies A ⊥ characterizes transverse size of the system. Noisy hydro equations of motion can be simplified with respect to the factorization. For later convenience, we expand the ξ dependence of a hydrodynamical variable into its conjugate k ξ , through a Fourier transformation, e.g., temperature is One finds that for each k ξ -mode, Eqs. (2.7) can now be recast into a Langevin equation, where prime indicates derivatives with respect to 3 ρ = ln(τ /τ 0 ) and We follow the same notations as in [21], such thatñ = dξe ik ξ ξ δs/s stands for the relative fluctuations of entropy density, andα = dξe ik ξ ξ τ u ξ the fluctuation of flow velocity along ξ. Γ is a 2 × 2 matrix which is entirely determined by Bjorken's solution of hydrodynamics [21]. K incorporates k ξ -mode of the random scalar function f (τ, ξ), 2 Since the bulk part is automatically factorized, one can simply restore the result in [21] for the ζ = 0 case with an extra term proportional to ∆ αβ 3 We use the same notation ρ, which should be distinguished from the de Sitter time ρ used in the 'hat' coordinate system for Gubser flow.
It is interesting to note that the effect of thermal noise onα vanishes for the k ξ = 0 mode.
Without transverse expansion, Bjorken flow in 1+1D is an oversimplified model regarding heavy-ion collisions. Nonetheless, it is worth analyzing, at least to a qualitative level, the effect of hydrodynamical noise. We start by rewriting Eq. (4.3) for each k ξ -mode as, Despite the constant ν, which is proportional to η/s, the amplitude of the autocorrelation is totally determined by the factor A ⊥ wτ in the denominator. Apparently one would expect a stronger effect of hydrodynamical noise in a system with smaller size, due to the appearance of A ⊥ . However, with respect to heavy-ion collisions, parametrically one can write the factor as a whole as which is interpreted as the transverse energy deposited per rapidity, and is equivalent to the multiplicity production of one collision event. Eq. (4.8) thus demonstrates the fact that the strength of hydrodynamical noise is essentially controlled by the multiplicity, instead of transverse size of a colliding system.

B. Gubser flow
Following the same strategy for Bjorken flow in the previous section, we generalize the discussion of hydrodynamical noise to the case of Gubser flow. Note that in the (ρ, θ, φ, ξ) coordinate system, flow velocity profile is essentially the same as a Bjorken's solution in 1+1D, except the fact that transverse expansion is now taken into account as well. Therefore, it is not surprising that the symmetry considered in Gubser flow leads to a similar factorization of the tensor structure in Eq. (2.11). Especially for Navier-Stokes hydrodynamicsΠ µν ∝P µν , one has∆ µναβ ∝P µνP αβ (4.10) and thus the autocorrelation of hydrodynamical noise becomes (4.11) With the above factorization, we can accordingly writeŜ µν in terms ofP µν , so that the correlation function of hydrodynamical noise reduces to correlation of a dimensionless scalar function, Note that enthalpy densityŵ(ρ) is used to saturate conformal dimension, even though quantities in the 'hat' system are dimensionless by construction. Accounting for the SO(3) × SO(1, 1) × Z 2 -symmetry of Gubser flow, the mode decomposition can be done with respect to spherical harmonics. In particular for the scalar function f (ρ, θ, φ, ξ), one has 4f and the two-point autocorrelation of each mode can be correspondingly found as A couple of comments are in order with respect to Eq. (4.15). First, it should be emphasized that hydrodynamical noise contains only scalar modes, as a direct consequence of the factorization Eq. (4.12). However, it is not generally true since the tensor structure of the hydrodynamical noise term can be more involved and leads to fluctuations in both vector and tensor modes. Second, one can parametrically estimate the effects of hydrodynamical noise in heavy-ion collisions, by examining the magnitude of correlation in Eq. (4.15), as we did previously in Section IV A for the Bjorken flow. The essential quantity that determines the strength of hydrodynamical noise isŵ in the denominator, in addition to the constant parameter ν which is proportional to η/s. In Gubser's solution to hydrodynamics,ŵ relies solely on the parameterT 0 , which is determined by the total multiplicity, but not system size. Therefore, in accordance with what was noticed in Section IV A, although one usually expects the system size to play a significant role in the estimate of hydrodynamical noise, we conclude that the absolute effect of hydrodynamical noise in a expanding medium in heavy-ion collisions is dominated by multiplicity.
When decomposing fluctuations of hydro variables, scalar modes and vector modes must be considered with respect to spherical harmonics Y lm and vector spherical harmonics Φ i lm respectively. We expand the following independent fluctuation variables, with extra dependence on the transverse dimension captured by indices l and m. The prime in Eq. (4.17) indicates derivative with respect to ρ, andṼ(ρ) is, (4.18) We purposely assign the vector mode to be the fourth element inṼ lm .Γ(l, m, ρ, k ξ ) in Eq. (4.17) is a 4 × 4 matrix with a rather complicated from, with its components given in [25].Γ is block-diagonalized, since vector modes and scalar modes are decoupled due to parity.K lm (ρ, k ξ ) depends on the scalar modes of hydrodynamical noise h lm , The last element inK vanishes as it should, since hydrodynamical noise does not contribute to vector modes. Also, we notice that for the k ξ = 0 mode, the evolution of v ξ lm (ρ, 0) is insensitive to the hydrodynamical noise as well. Eq. (4.19) introduces two sources of instabilities, corresponding to the zeros ofT andT + H 0 tanh ρ. As has been discussed throughly in [25], both sorts of instabilities have been noticed in the structure of the matrix Γ already in the case without hydrodynamical noise.

V. SOLVING GUBSER FLOW WITH HYDRODYNAMICAL NOISE
Hydro equations of motion are coupled with the equation of state, which in general resorts to numerical solutions. When hydrodynamical noise is taken in account, in addition to the ensemble-averaged background flow, one also needs numerical simulations of stochastic equations [22,28] regarding the random unknown tensor variables S µν . However, with respect to a background Gubser flow, the tensor structure of hydrodynamical noise is well determined and has been discussed in the previous section, so that numerical simulations are largely simplified.
A general procedure of solving noisy hydrodynamics on top of Gubser flow comprises the following steps. First, one solves Eqs. (4.17) mode-by-mode, with initial conditions and parameters specified with respect to desired collision systems. Second, hydrodynamical variables, such as flow velocity and temperature, are obtained through mode summation as in Eqs. (4.16). Similar strategy has also been applied in analyses of perturbations on top of solvable hydro models [29][30][31]. To recover quantities in the original (τ, r, φ, ξ) coordinate system, one is additionally required to do a coordinate transformation according to Eqs. (3.2).
For each mode, although Eqs. (4.17) represent four-coupled equations, simplifications can be made as follows. First, as we emphasized before, the hydrodynamical evolution does not couple vector modes to scalar modes. Besides, hydrodynamical fluctuations do not contribute to vector modes. As a consequence, we shall not consider vector modes in our analysis. Second, we shall restrict ourselves to the case of k ξ = 0. By doing so, the matrix Γ is further block-diagonalized in a way that equation of motion for v ξ is decoupled. The third component ofK also vanishes when k ξ = 0. Therefore, the only non-trivial equations of motion which receive extra contributions from hydrodynamical noise are the two coupled equations for the scalar modes δ lm and v s lm , with Similar equations have been investigated in [31] without hydrodynamical noise. Note that all m-modes of the same index l evolve identically since there is no dependence on index m in theΓ matrix andK.

A. Formal solution
Before numerically solving Eqs. (4.17), we investigate qualitatively the behavior of hydrodynamical noise and its evolution. The following discussion is made ad hoc with respect to a background Gubser flow, but it can be applied to hydrodynamical noise on top of Bjorken flow as well. One can write a formal solution of Eq. (4.17) (and similarly Eq. (4.5)) in terms of a Green function, where K is an abbreviated notation for the conjugate variables to specify modes, i.e., K = (l, m, k ξ ) for Gubser flow and K = k ξ for Bjorken flow. Green functionG is determined by the following equation (note that we take ρ = ln(τ 0 /τ ) for Bjorken flow) with the initial conditionG (0, K) = 1 .
where T indicates a time ordering with respect to ρ. Except in some extreme limits, there is no simple analytical expression for the Green functionG, regarding a specified form ofΓ from hydrodynamics. However, some of the behaviors of mode evolution are known qualitatively. For instance, it has been shown that viscosity damps mode evolution [25]. In particular, in the large l and small viscosity limit, one finds thatG(∆ρ) ∼ exp[−l 2 H 0 ∆ρ].
From Eq. (5.2), we notice that at any time, each mode comprises contributions from fluctuations in the initial state (the second term) and hydrodynamical noise (the first term). Especially, due to the fact that the ensemble average of one-point function of thermal fluctuations vanishes, K = 0 , the mode evolution of the one-point function is governed by hydrodynamic response to initial state fluctuations. In heavy-ion collisions, initial state fluctuations are quantum fluctuations associated with the probability distribution of nucleons inside the colliding nucleus. Therefore, one would expect initial state fluctuations to be independent of ensemble average, and the one-point function evolution has a form of linear response The effect of hydrodynamical noise can be investigated in terms of the evolution of twopoint correlation function. The equal-time two-point correlation can be found according to Eq. (5.2) More explicitly, collision events are distinguished by their initial conditions, while ensemble events stick with the same initial condition but evolve with random hydrodynamical noise on an event-by-event basis. Therefore, we are allowed to drop the mixing between hydrodynamical noiseK at any time and initial state fluctuationsṼ(ρ 0 , K) in the average, namely K (ρ, K)Ṽ(ρ 0 , K ) = 0, which factorizes into the ensemble average of one-point function of hydrodynamical noise, and the two-point correlation in Eq. (5.7) is written again as a term from initial state fluctuations plus a term from hydrodynamical noise.
The structure of two-point autocorrelations of hydrodynamical noise is known locally in space-time, with the amplitude Λ th fixed with respect to the form ofK given in Eq. (4.19) and the fluctuation-dissipation theorem in Eq. (4.15). Note that vector modes do not contribute. Two-point correlations of initial state fluctuations are determined by averaging over collision events, instead of ensemble average. We consider two extreme scenarios in this work. One is inspired by independent sources [32], from which two-point correlations of initial state fluctuations are expected to be local in the transverse plane, and all K-modes are initialized correspondingly with specified values. In the other scenario, we initialize with selected modes accounting for the SO(2) rotational symmetry in the transverse plane. Especially, we are allowed to deform initial energy density profile with a desired eccentricity. More details of these two types of initialization will be given later in the next section.

B. Numerical simulations
We solve the two coupled equations of motion numerically on top of Gubser flow for the scalar modes δ l and v s l , with parametersT 0 and q specified with respect to ultra-central collision events of Pb-Pb, p-Pb and p-p at the LHC energies.T 0 is determined according to the event multiplicity. Following [24] and [31], we adopt the relation, where the constant f * = /T 4 = 11 is the effective degree of freedom extracted from Lattice QCD calculations, and dS dξ = 7.5 dN ch dy . (5.10) T 0 is found to be 7.3 for the 0∼5% central PbPb collisions with √ s N N = 2.76 TeV [31], corresponding to multiplicity production per rapidity dN ch /dy ∼ 1600 [33]. For the p-Pb collisions with √ s N N = 5.02 TeV [3], we takeT 0 = 3.1 which leads to an estimate of dN ch /dy ∼ 150. Very recently, long-range correlations were measured in the ultra-central proton-proton collision events with multiplicity higher than 100 in the rapidity gap |y| < 2.4 [1], for which we takeT 0 = 2.0. The parameter q constrains the finite transverse size of the fluid system. For a Pb-Pb system, we take q = (4.3fm) −1 [24], while for both p-Pb and p-p systems we take q = (1.1fm) −1 . In all the simulations in this work, the viscosity parameter is taken to be H 0 = 0.33 [25], which corresponds to η/s = 0.134. Fig. 1 displays the solved temperature distribution of one random Pb-Pb, p-Pb and p-p event at τ = 2.5 fm, without initial state fluctuations. Specifying τ = 2.5 fm for the analysis of all the three systems breaks conformal symmetry, which is however of phenomenological interest since τ = 2.5 fm is a typical time scale that all the three systems experience in the early stages of evolution. The effect of hydrodynamical noise is qualitatively captured by the bumpiness of the temperature profile, which presents a clear trend of becoming more pronounced from the Pb-Pb system ( Fig. 1(a)) to the p-p system (Fig. 1(c)). We have checked that when taking ensemble average, the bumpiness disappears in the temperature evolution, which corresponds to δ l (ρ) = 0.
In order to quantify the effect of hydrodynamical noise, one must investigate two-point auto-correlations of hydrodynamical variables. For the sake of numerical simplicity, in this work we focus on the two-point auto-correlation of radial flow velocity u r of the same τ (equal time) and r (equal radius), which is determined by the equal-time two-point autocorrelations of modes v s l (ρ)v s l (ρ) . We have checked that our conclusions are not changed from the analysis of other types of two-point auto-correlations, such as the auto-correlations of temperature fluctuations which depends on δ l (ρ)δ l (ρ) . We thereby define the following correlation function that describes the two-point auto-correlation of radial flow velocity of the medium, C urur (τ, ∆φ, r, φ) = u r (τ, r, φ)u r (τ, r, φ + ∆φ) − u rb (τ, r) 2 . (5.11) Note that trivial contributions from the azimuthally symmetric background flow are subtracted in Eq. (5.11). One can further focus on the correlation structure with respect to the relative azimuthal angle ∆φ by integrating over r and φ, The angular structure in Eq. (5.12) depends not only on the hydrodynamical noise, but also contains a fraction induced from initial state fluctuations, which we denote as C T urur and C I urur respectively. The significance of hydrodynamical noise in heavy-ion collisions is then captured by the relative contributions from C T urur and C I urur . Therefore, in our numerical results, we shall concentrate on the ratio between these two contributions, C T /C I .
It should be emphasized that the separation of two-point radial flow velocity autocorrelation into one term from initial state fluctuations and the other induced by hydrodynamical noise in Eq. (5.12) is one particular feature in our analysis, as has been demonstrated in the formal solution Eq. (5.7). In practical analysis, one is thus allowed to simulate the system evolution independently in cases with only initial state fluctuations, and in cases with only hydrodynamical fluctuations. However, in more involved studies where contributions to two-point correlations from initial state fluctuations and hydrodynamical fluctuations are not simply separable, one must carry out simulations with both sources of fluctuations considered simultaneously. In this work, we follow the conventional procedure of solving noisy hydro equations of motion, by initializing the system with initial state fluctuations in two extreme scenarios: One with a specified initial azimuthal anisotropy and the other with a Dirac delta function. By doing so, we claim that effects of hydrodynamical noise can be estimated with respect to the experimentally measured long-range and short-range correlation structures respectively.

Effects of hydrodynamical noise on long-range correlations
Long-range correlations of the observed spectra in heavy-ion collisions are associated with harmonic flow, which in turn depends on the evolution of azimuthal anisotropies of density  2. (Color online) Structure of the radial flow velocity auto-correlation determined by hydrodynamical noise, dipicted as the ratio between C T urur (τ, ∆φ) and the magnitude of C I urur (τ, ∆φ) (i.e, the absolute value of C I urur (τ, ∆φ) at ∆φ = 0) at τ = 2.5 fm. Results are obtained from numerical simulations of 10000 events with respect to initial condition Eq. (5.13) with an (a) ε 2 , (b) ε 3 , (c) ε 4 and (d) ε 5 , for p-p (red solid lines), p-Pb (orange dashed lines) and Pb-Pb (green dash-dotted lines) systems.
profile. For each harmonic order m, anisotropy is characterized by the so-called eccentricity ε m which represents an invariant deformation of the density profile under azimuthal rotation φ → φ + 2π/m. Therefore, we perturb the azimuthally symmetric initial profile of Gubser flow by certain azimuthal modes. Following the discussions originated in [25], for each harmonic order m, when initial state fluctuations are characterized by the following nonzero modes 5 δT (θ, φ, ρ 0 , ξ) one accordingly realizes a density profile with a non-zero ε m . It should be emphasized that the eccentricity generated from Eq. (5.13) is closely related to the eccentricity defined in terms of cumulants with a r m weighting [34], yet to an approximate level. More detailed analysis of Eq. (5.13) with respect to phenomenology can be found in [25]. The sign convention in Eq. (5.13) is taken so that the gradient of the deformed density is maximal along x-axis. Λ ini in Eq. (5.13) reduces to a constant parameter to be fixed by the values of ε m in the colliding systems. In our simulations with respect to the ultra-central Pb-Pb and p-Pb colliding systems, we take ε 2 (Pb-Pb) ∼ 0.05 and ε 2 (p-Pb) ∼ 0.15 at τ = 0.6 fm, which are typical values considered in phenomenological studies of heavy-ion collisions (cf. [35,36]). For ultra-central p-p collisions, due to the smaller multiplicity production we take a larger value of initial anisotropy ε 2 (p-p) ∼ 0.2. We have checked that the same values of Λ ini result in eccentricities of higher harmonics (up to ε 5 ) of the similar order of magnitude in the corresponding systems.
For each of the particular deformations introduced in the initial state, we solve for the noisy Gubser flow and calculate the correlation function defined in Eq. (5.12). Since we are only interested in the effect of hydrodynamical noise on the long-range correlations, i.e., the evolution of azimuthal anisotropies, in the mode summation we ignore contributions from modes other than those of relevance to the corresponding initial anisotropies. For instance, if one calculates the evolution of ellipticity, to quantify Eq. (5.12) the mode summation only involves (v s 2,2 (ρ)) 2 and (v s 2,−2 (ρ)) 2 . In Fig. 2, two-point auto-correlation C T urur (τ, ∆φ) are plotted as a function of ∆φ at τ = 2.5 fm. Although the structure of C I urur (∆φ) is not shown, it is worth mentioning that C I urur and C T urur share the same structures as a function of ∆φ, but are different in magnitudes. The periodic correlation structures shown in Fig. 2 are rooted in the azimuthal symmetries considered for each of these cases. For an initial ε 2 , ε 3 , ε 4 and ε 5 , correlations of radial flow velocity exhibit periodicity in π, 2π/3, π/2 and 2π/5 respectively. We take the ratio between C T urur and the magnitude of the correlation function C I urur , i.e., C T urur (τ, ∆φ)/C I urur (τ, 0), so that one is allowed to read off directly the relative change of anisotropies due to hydrodynamical noise in Fig. 2. As can be seen in Fig. 2, hydrodynamical noise results in extra contributions to the development of azimuthal anisotropies, which are getting stronger from the ultra-central Pb-Pb collision systems, to ultra-central p-Pb and p-p, and also from lower order harmonics to higher order harmonics. The increasing contribution from hydrodynamical noise according to harmonic order can be understood as follows: On one hand, hydrodynamcial noise is insensitive to the harmonic order, which is subject to the independence of index m in the evolution equation (cf. Eqs. (5.1)). On the other hand, however, evolution of anisotropies of higher order harmonics suffers stronger viscous suppression. Accounting for both effects, one would expect that relatively hydrodynamical noise becomes more important for higher order harmonics. Nonetheless, the over-all magnitude of enhancement is not significant, in particular for the Pb-Pb systems in which it is less than 3%.

Effects of hydrodynamical noise on short-range correlations
For the analysis of short-range correlations, initial condition is chosen with temperature fluctuating in the transverse plane with respect to a Dirac delta function, (Color online) Two-point auto-correlations of radial flow velocity defined by only an integration over φ in Eq. (5.11), C urur (τ, ∆φ, r), as a function of r and ∆φ. Results are obtained with respect to initial condition Eq. (5.14) with θ 0 = 0.2 and φ 0 = π, at τ = 2.5 fm, without hydrodynamical noise (upper row) and with hydrodynamical noise (lower row) from numerical simulations of 5000 events. To have all the three systems plotted on the same scale, magnitudes of the correlations are rescaled in a way that the height of the bump at |∆φ| = π is equal to 1.
One may check that on an event-by-event basis, Eq. (5.14) is consistent with a Dirac delta in the transverse plane in the original Milne space-time. 6 The same values of Λ ini fixed by eccentricities in the previous section are kept regarding the colliding systems under consideration. Without losing generality, we take θ 0 = 0.2 and φ 0 = π throughout this work, which corresponds to a peak near the origin on the x-axis. Since higher l-modes receive stronger viscous suppression, our simulations are limited to the modes l < 30. Fig. 3(a), (b) and (c) depict the corresponding results of the correlation function C urur (τ, r, ∆φ) (without r integration) in the ultra-central Pb-Pb, p-Pb and p-p colliding systems respectively without the inclusion of hydrodynamical noise, at τ = 2.5 fm. Regarding an initial density profile which is perturbed by a delta function, hydro evolution results in sound-wave propagation along with the medium expansion, which is exactly seen as a bump in the over-all structure of C urur . To have the correlation functions plotted on the same scale, we have rescaled the height of the bump so that it is unity at |∆φ| = π. This rescaling reveals a generic feature of Gubser's solution that the medium expands much faster in smaller systems p-Pb and p-p than in Pb-Pb, which however does not affect our analysis of the effect of hydrodynamical noise. The position of the sound horizon reflects hydro response to the position of the initial delta function. The shape of the bump is altered from a delta function by diffusion and dissipation. As one would expect, hydrodynamical noise affects the fine structure of the fluid system, which is then reflected as an excess around (Color online) Details of the near-side peak in radial flow velocity correlations due to hydrodynamical noise, depicted in terms of the correlation function C T urur (τ, ∆φ) divided by the magnitude of correlation due to initial fluctuation, i.e., value of C I urur (τ, ∆φ) taken at |∆φ| = π, with τ = 2.5 fm. ∆φ = 0. As seen in Fig. 3(d), (e) and (f), when hydrodynamical noise is considered in our simulations, an excess in the auto-correlation structure indeed appears at ∆φ = 0 and persists from around origin (r = 0) to large radii. Especially, on top of the bump structure the hydrodynamical noise leads to a near-side peak which is marginal in the Pb-Pb system as shown in Fig. 3(d), but becomes significant in p-p as shown in Fig. 3(f).
We next concentrate on the details of the near-side peak from the hydrodynamical noise in Fig. 3. To avoid overestimation of the effect of hydrodynamical noise, we limit the integration over r and φ around the bump. Rescaled again by the magnitude of C I urur (now with ∆φ = π), we obtain consequently the near-side peak in C T urur due to hydrodynamical noise relative to the correlation structure induced by initial state fluctuations. Results corresponding to the ultra-central Pb-Pb, p-Pb and p-p systems are shown in Fig. 4, which exhibit a clear trend in the three colliding systems. The height, as well as the width, of the peak increase from Pb-Pb, to p-Pb and p-p.

VI. SUMMARY AND CONCLUSIONS
In this work we solved the noisy Gubser flow, by implementing hydrodynamical noise on top of Gubser's solution to Navier-Stokes hydrodynamics. Hydrodynamical noise is formulated in a standard way according to the fluctuation-dissipation relations, from which we noticed that the absolute amplitude of hydrodynamical noise in ultra-central heavy-ion collisions is essentially determined by the total multiplicity of the collision event, instead of the system size as one might have anticipated. Regarding the ultra-central Pb-Pb, p-Pb and p-p collisions carried out at the LHC energies, we quantitatively analyzed the effects of hydrodynamical noise with emphasis on long-range (evolution of azimuthal anisotropies) and short-range correlations. A clear trend of enhancement of the hydrodynamical noise was confirmed in both cases from the Pb-Pb system, to p-Pb and p-p systems, which as we have claimed, is mostly due to the decrease of multiplicity, but not system size.
Azimuthal anisotropies receive extra contributions from the hydrodynamical noise during medium evolution, which implies an increase of harmonic flow v n . In addition, higher order harmonics are found more sensitive to the hydrodynamical noise, which confirms the results from more sophisticated hydrodynamical simulations [23]. However, at least from our simulations for the conformal and azimuthally symmetric systems, the increase of anisotropies due to hydrodynamical noise is not significant. Especially, we expect that for the ultracentral Pb-Pb collision systems, the effect of hydrodynamical noise on harmonic flow v n is negligibly small. On the contrary, short-range structure of the fluid system is more affected by the effect of hydrodynamical noise, which is demonstrated as the appearance of a near-side peak. It is understandable in the sense that higher order hydro modes, which correspond to the finer structure of a fluid system, are dominated by contributions from hydrodynamical noise, rather than hydro response to initial state fluctuations. We investigated in details the structure of of the peak, relative to the short-range two-point flow velocity correlations induced by initial state fluctuations. We noticed that both the height and the width are enhanced when the effect of hydrodynamical noise becomes larger.
Bearing in mind that our analysis of the noisy Gubser flow is less realistic in several aspects, especially the caveat that we have no freeze-out process and there exists a gap between the two-point flow velocity correlation and the experimentally measured two-particle correlations, it is implied from our results that simulations of viscous hydrodynamics without hydrodynamical noise are reliable for the investigation of harmonic flow v n , even in the smaller colliding systems like p-Pb. However, hydrodynamical noise must be taken into account if one studies the fine structure of the near-side two-particle correlations, which also contains informations related to the dissipative properties of the QCD medium. However, ambiguity arises when one takesρ i = ρ + i∆ρ/n orρ i = ρ + (i + 1/2)∆ρ/n in the sum, which corresponds to the Itô integral and the Stratonovich integral respectively. For multiplicative noise (whenK depends onṼ) these two integrals lead to different values, while for noise that is not multiplicative (whenK is independent ofṼ) Itô and Stratonovich integrals coincide. Since hydrodynamical noise in our work is not multiplicative, recipes derived from either the Itô integral or the Stratonovich integral can be applied equivalently. Indeed, we have checked that the Euler-Maruyama method and the Heun's method discussed in Ref. [28] result in compatible solutions.