{'fraction_non_alphanumeric': 0.037883311672345356, 'fraction_numerical': 0.010660718676455091, 'mean_word_length': 4.302570311044887, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Differential HBT Method to Analyze Rotation

Two particle correlations are studied in the reaction plane of peripheral relativistic heavy ion reactions where the initial state has substantial angular momentum. The earlier predicted rotation effect and Kelvin Helmholtz Instability, leads to space-time momentum correlations among the emitted particles. A specific combination of two particle correlation measurements is proposed, which can sensitively detect the rotation of the emitting system. Here the method is presented in simple few source models where the symmetries and the possibilities of the detection can be demonstrated in a transparent way.


I. INTRODUCTION
Collective flow is one of the most dominant observable features in heavy ion reactions up to the highest available energies, and its global symmetries as well as its fluctuations are extensively studied. Especially at the highest energies for peripheral reaction the angular momentum of the initial state is substantial, which leads to observable rotation according to fluid dynamical estimates [1]. Furthermore the low viscosity quark-gluon fluid may lead to to initial turbulent instabilities, like the Kelvin Helmholtz Instability (KHI), according to numerical fluid dynamical estimates [2], which is also confirmed in a simplified analytic model [3]. These turbulent phenomena further increase the rotation of the system, which also leads to a large vorticity and circulation of the participant zone. [4]. In ref. [2] it is estimated that the increased rotation can be observable via the increased v 1flow, but the v 1 signal at high energies is weak, so other observables of the rotation are also needed.
The two particle correlation method is used to determine the space-time size of the system emitting the observed particles, thus providing valuable information on the exploding and expanding system at the freeze out stage of a heavy ion collision. This method is based on the Hanbury Brown and Twiss (HBT) method, originally used for the determination of the size of distant stars [5]. In heavy ion collisions the HBT method was used first for the same purpose, the determination of the system size [6], but later also the ellipsoidal shape of the system and its tilt [7,8]. It was also observed relatively early that the expansion of the system modifies the size estimates due to the collective radial flow velocity of the emitting system [9], while the effect of flow on two particle correlations was also analysed in great detail [10]. Transport model studies have indicated that the HBT radius shows a minimum at the phase transition threshold [11].
A detailed two particle correlation study of flow rotation was not performed up to now [12]. This became actual as at higher beam energies, where the initial angular momentum of the participant system is increasing in peripheral reactions, the system may rotate, causing a significant and detectable effect. Here we do not want to use the most advanced state of the art developments of the HBT method, rather we want to show in the most simple way, how rotation can be detected by the method.
Different theoretical approaches were worked out up to now to evaluate the two particle correlation functions in different reaction models. We here picked one method, which is not used very frequently up to now, but it can be generalized and used to fluid dynamical models with arbitrary flow patterns well. With this method hereby we study simplified, idealized fluid dynamical systems with different symmetry structures. These studies show how we can detect rotation via two particle correlation functions and what effects may cause difficulties in identifying rotation.
Based on these studies we also present a Differential Hanbury Brown and Twiss (DHBT) method, which can sensitively determine the strength (vorticity or circulation) of the rotating flow and the direction of this rotation.

A. The Emission Function
Following [13] the definition of a particle 4-current is: where f (x, p) is the invariant scalar phase space density distribution of the emitted particles. The total flow of N nucleons across a space-time (ST) hypersurface, with the surface element, dσ µ is: If we do not perform the integration over the momentum, p, then we get the Cooper-Frye formula [14] for invariant momentum distribution: 2 where E = p 0 . This assumes that there is a 3-dimensional Freeze Out (FO) hypersurface, which can also be generalized to become a layer, as we will show. Using d 3 p = dp || dp ⊥ = p 0 dydp ⊥ ) we can write this in another form, and we can extend it to a 4-volume integral of a 4-dimensional Source Function, S, as [15,16]: where we assumed that the emission appears on a 3dimensional hypersurface with the outward pointing normal, dσ µ . The emission function gives the distribution of the ST positions of momenta of emitted particles. The emission function gives the number of the particles, ∆N , emitted in the phase-space element ∆ 3 x∆ 3 p per unit time, ∆t. We get a Lorentz invariant scalar if we also multiply it by the energy, p 0 , of the emitted particles: Along the lines of this introduction we can describe the emission over a hypersurface also as a 4-volume integral. We still assume a ST hypersurface with an outward pointing normal vector, dσ µ , but then the ST integral is interpreted as a 4-volume integral with a delta function for the given surface, δ(x − x) as: where the emission is constrained to a 3-dimensional hypersurface in the ST and it is directed in a given direction characterized by the unit vector in the source described byσ µ (x ) δ (4) (x ). Even if we assume that the source is in a ST layer (which is not too thick, e.g. 2-3 fm or 2-3 fm/c), in this layer we can have a maximum of the emission. This FO layer in case of pions is narrow if we consider the rapid and simultaneous hadronization and FO from the plasma. This could even be idealized as a 3dim hypersurface in the 4-dim ST if the thickness of the layer is neglected. Thus,σ µ , is the unit normal vector of this surface or layer can be both timelike or spacelike: σ µσ µ = +1 for timelike hypersurface or layer , σ µσ µ = −1 for spacelike hypersurface or layer . (7) The idealization of FO in a 3D hypersurface is not necessary, however it makes the presentation more transparent. The more realistic emission distribution must happen in a 4D ST surface layer, which can still can have an effective space-like or time-like normal vector. See refs. [10,[17][18][19]. Although, one might naively believe that in case of a normal vector,σ µ = (1, 0, 0, 0) the emission is uniform in all spatial directions, this is not true, as the local flow velocity also influences the emission probability [10,18]. When the flow velocity points in the direction of the detector (or the FO normal) the probability for the emission into theσ µ -direction is bigger, so that the emission probability should be proportional to (u µ ·σ µ ). This will be important later on when the observability of rotation is discussed.
At ultrarelativistic energies most of the FO happens in time like directions, frequently idealized as t =constant or τ =constant hypersurfaces or layers. In addition even in case of time-like FO the deeper (or earlier) points of the FO layer have a smaller emission probability because of the opacity of QGP, indicated also by the strong jet quenching. This effect causes additional asymmetries in the emission.

B. Hydrodynamical Parameterization
Let us assume that at the points of the source the matter is still in local equilibrium. Then we can describe the phase-space distribution of the particles, f (x, p) by a Jüttner distributions [13] or a relativistic Bose-Einstein or Fermi-Dirac distribution. Furthermore, instead of the delta function we may assume the emission distributed in a ST layer, which still has a preferred direction of emissionσ µ .
For our purposes the most suitable parametrization of the emission function is introduced in the special "Buda-Lund" model, see section 8 of ref. [20] . Here the emission function is parametrized as where the emission probability and its dependence of the FO direction is already included in the p µ d 4 Σ µ (x) term. The 4-volume integral is directed and yields a maximum for k µ which is closest to d 4 Σ µ (x). In the Buda-Lund model the FO direction points into the flow 4-velocity, d 4 Σ µ (x) ∝ u µ . This is also a frequent approximation in other fluid dynamical models, although, in the general case it is not a valid approximation. Thus, this approach is identical to the dynamical volume FO in a layer [17][18][19], discussed above. We also assume that the FO layer depends parametrically on the proper time, τ , (or distance s) in the FO factor d 4 Σ µ (x) ∝σ µ d 4 x. So, that the emission probability is proportional to G(x) H(τ ): where and that the widths of the emitting sources do not change significantly during the course of emission from a given source. (The emission probability will be discussed in more detail when the source asymmetry is analysed in section III.) The emission function characterized with a locally thermalized volume-emitting source is then: where in place of f (x, p) we inserted the relativistic Bose-Einstein distribution. The factor g stands for the degeneracy, u µ (x) is the 4-velocity field, T (x) is the temperature field, µ(x) is the chemical potential and d 4 x = dτ dx dy dz. G(x) is the ST emission density across the layer of the particles (e.g. pions). This can be approximated with a ST hypersurface and then with Θ → 0 we obtain the Cooper-Frye FO description. [14] For the phase space distribution we frequently use the Jüttner (relativistic Boltzmann) distribution: which is normalized to the invariant scalar density of particles where C n = 4πm 2 T K 2 (m/T ) and we use the c = k = 1 convention. Thus in terms of the local invariant scalar particle density the Jüttner distribution is [13] f We can also define dσ µ = dx dy dzσ µ (for a timelike surface or layer) where the norm of dσ µ is the 3-volume of the source element (like a fluid cell), similarly to Refs. [17][18][19] Our source function in this case, with Eq. (13), similarly to Ref. [20] in the frame whereσ µ = (1, 0, 0, 0) is: where N F D (x) is the density of particles arising from fluid-dynamical (FD) evolution (or from other transport models). Then taking inti account the flow velocity in this frame: The H(τ ) (or H(s)) freeze-out probability along the τ (or s) parameter, across the layer can be integrated separately from the remaining three orthogonal coordinates.
Here we assume that the primary direction of the emission isσ µ . Thus, in case of an explosively expanding system it points towards the detector, soσ µ ≈k µ . The emission happens from the 4-volume of ST surface layer with an effective normal directionσ µ , and not from a ST hypersurface.
In ref. [10] the correlation function was analyzed in detail in dependence of the direction of the primary emission direction, k µ . A possibility of longitudinal momentum difference was considered in terms of the rapidity of the two emitted particles, where the difference of the longitudinal momenta, the width parameter, was studied in detail. While this analysis provides a deeper insight into the features of the correlation function, in our studies we concentrate to the rotation of the system and restrict ourself to a simplest presentation of the correlation function, which can be realized experimentally without much additional effort.
We can assume FO from a narrow layer at a proper time hyperbola τ F O = cons. like in the Buda-Lund model, see ref. [20]. This can be practical if the CFD model uses proper time and rapidity coordinates. Recent studies indicate that irrespective of the coordinate system choice, in the major part of FO in high energy collisions the FO happens near to a constant proper time hyperbola, although the origin of this FO-hyperbola is at an earlier point of time than the intersection of the centers of the projectile and target trajectories [21].
For the first test purpose we take an oversimplified model of 4 fluid elements, which may or may not expand or rotate. We assume that these are in the reaction plane, [x − z]-plane, and will characterise parameters of a heavy ion reaction based on CFD results. We assume that the system is stationary so the time emission probability is a Gaussian (like) distribution in time. Later on we intend to study and see that these methods can be applied for realistic full scale FD calculations also.

C. Pion Correlation Functions
The pion correlation function is defined as the inclusive two-particle distribution divided by the product of the inclusive one-particle distributions, such that [15]: where p 1 and p 2 are the 4-momenta of the pions. We will assume pions are created at two points, x 1 and x 2 , which are distributed in space. The particle distribution is given by the reduced phase space source distribution: For two identical pions with momenta p 1 and p 2 the twoparticle distribution is: where the wave equation ψ 12 is given by: We now introduce the center-of-mass momentum 1 and the relative momentum where assuming the mass-shell constraint for the two particles and so we have q · k = (p 1 − p 2 ) · (p 1 + p 2 )/2 = (p 2 1 − p 2 2 )/2 = (m 2 π − m 2 π )/2 = 0, which leads to q 0 = q · k/k 0 With the relative and center-of-mass momentum we can write the wave equation as: and then We can then insert Eq. (23) into Eq. (18), and we obtain: where the last term in the brackets is cos[q(x 1 − x 2 )]. Similarly for the one particle distribution we get: We now use a method for moving sources presented in ref. [22]. Using Eqs. (24,25), together with the definition of the correlation function we have: where Here R(k, q) can be calculated [22] via the function and we obtain the R(k, q) function as This can be verified, by using Eq. (28), forming a double integral over d 4 x 1 d 4 x 2 from J(k, q) J(k, −q), yielding to a term exp[−iq(x 1 − x 2 )]. Then taking the real part of the double integral leads to a term cos[q(x 1 − x 2 )] and this recovers Eq. (27).
Source with local Jütner distribution: Let us take the S(x 1 , p 1 )S(x 2 , p 2 ) term in Eq. (24), and assume that the single particle distributions, f (x, p), in the source functions are Jüttner distributions, which depend on the local velocity, u µ (x) , via the term as shown in Eq. (13). Here the local flow velocity may be different in different locations, x 1 and x 2 , and this influences the correlations of the observed momenta. Thus, the scalar products in terms of k and q become: where we used the notation u 1 = u(x 1 ) = u µ (x 1 ). We assume that for a given detector position the normal direction of the emission is approximately the same, so for the two sources the term p µσ µ (x) is the same and it cancels in the nominator and denominator.
Thus, the expression of the correlation function, Eq. (27) will be modified to and the corresponding J(k, q) function will become In Eq. (32) the term including the momentum component q and the flow velocity u becomes unity if the source has a uniform distribution of u(x)/T (x), and in this case we may be able to use the so called smoothness approximation: S(x, k+q/2)S(y, k−q/2) ≈ S(x, k)S(y, k), and the correlation function, expression (26), takes the form

D. One Fluid Cell as Source
We now assume a source function, which is reduced to one Freeze Out (FO) time moment. Thus the integration over the 4-volume of an emission layer is reduced to the 3-volume of a FO hypersurface. For simplicity, we assume FO along the timelike cordinate, t, where we assume a local Jüttner distribution. Thus, we have the source function as where k µσ µ is an invariant scalar, and for a single cell we use a simple quadratic parametrization for n(x) as: Here n s is the average density of the Gaussian source (or fluid cell) of mean radius R.
Single source at rest: The invariant scalar k µ u µ can be calculated in the frame where the cell is at rest. We have then In this simplest case we also assume that the FO direction isσ µ = (1, 0, 0, 0), so the τ -coordinate coincides with the t-coordinate, and it is orthogonal to the x, y, x− coordinates. Then we can make use of the following integral: We can perform the integral along the t direction of H(t), which gives unity and then the single particle distribution is where T s is the temperature of the source, and E k = k 0 in the rest frame of the fluid cell. Due to the normalization of H(t) the integral over the time t is unity. The contribution to the nominator from Eq. (33) is where we used [23] 3.323/2. In the time integral the present choice ofσ µ would give (q 0 ) 2 , but we wanted to indicate that other choices are also possible and they would yield (σ µ q µ ) 2 . In the J(k, q)J(k, −q) product the terms exp[±q 0 /(2T s )] cancel each other. Inserting these equations into (26) we get If we have a source at a point in the FO layer, which is at a longer distance from the external side of the FO layer than Θ, then the contribution of the time integral from this point is reduced. In a few source model it is more transparent to describe this reduction by assigning a smaller weight factor to the contribution of the deeper lying source. See in section III.
If we tend to an infinitely narrow FO layer, Θ → 0, i.e. to a FO hypersurface, then The k dependence drops out from the correlation function, C(k, q) as the k dependent parts are separable. See Fig. 1. The size of the fluid cells in a high resolution 3+1D fluid dynamical calculation is (0.3fm) 3 . With this resolution the numerical viscosity of the fluid dynamical calculation [24] is the same as the estimated minimal viscosity of the QGP [25] which occurs at the critical point of the phase transition [26]. As Fig. 1 shows the correlation for such a cell size yields to an extended distribution in the relative momentum q. For the study of the rotation of the system the thickness of the FO layer is of secondary importance, especially if we discuss only a few fluid sources. In this case the role of the depth of a source point within the layer is given by its reduced contribution to the particle emission. This can be represented much simpler with assigning emission weights to the small number of sources. Thus, in the following discussion, we do not go into the details of the time structure of the emission.
Single moving source: Let us take a single source which moves in the x-direction with a velocity v x . Then we have, u µ s = γ s (1, v x , 0, 0), and the scalar product k · u s /T s = k µ u µ s /T s provides an additional contribution to the correlation function. However, in the case of a single fluid cell or a single source the velocity and the temperature do not change within the cell, so the modifying term in eq. (32) becomes unity. We use , and the source function becomes where Within the source (or fluid element) the velocity u s and temperature T s are assumed to be the same. The source or fluid element may have a density profile, but this profile should be the same for all cells (although the average density, n s is not the same for all cells. The spatial integrals can be performed in the rest frame of the cell, giving the same integral result as above (39), because the moving cell-size shrinks, but the apparent density increases, so that the total number of particles in a cell remains the same as it is an invariant scalar. Then the integral of the single particle contribution is Then the two particle distribution: When calculating R(k, q), in the J(k, q)J(k, −q) product the terms exp[±q · u s /(2T s )] cancel each other. In the formulae the = 1 convention is used and k and q are considered as the wavenumber vectors. We then insert these equations into equation (26) and we get for one moving Gaussian source Again, this result does not depend on k, just as the previous single source at rest, Eq. (42).

A. Two Steady Fluid Cells
For emission from two steady sources, two particle correlations were studied in ref. [20]. Here we use the present method. We assume that the two source system is symmetric both their positions are placed symmetrically and also their FO normal vectors,σ µ , are the same. If the normal wereσ µ = (1, 0, 0, 0), then the invariant scalar k µσ µ would be k 0 = E k , although we do not need this additional requirement to illustrate the correlation function, which would arise from an idealized symmetric system.
We also assume that the time distributions, H(τ ) for the two sources are identical, so these can be integrated simultaneously and yield unity. If we have two sources then the source function is while the J function in the Jüttner approximation is where x s is the position of the center of the source, and the spatial integrals run separately for each of the identical sources, i.e. we assume fluid cells with identical density profiles, but with different densities, n s and temperatures, T s . In case of steady sources u µ s = (1, 0, 0, 0), and the spatial integral for one source is the same as for a single source. Thus, and In the J(k, q)J(k, −q) product the terms exp[±q 0 /(2T s )] cancel each other. Both J(k, q) and J(k, −q) include a sum [exp(iqx s )+exp(−iqx s )], and their product leads to a factor 2[1+cos(2qx s )]. Here we assumed that the timelike extent of the emission layer is negligible compared to the space-like size. Consequently, if the two sources have the same parameters, just different locations, This result agrees with ref. [20], section 9.1 (p. 41), and in the limit of x s = 0 it returns the single source result, Eq. (42). See Fig. 3. If the distance of the two sources is 2d, i.e. x 1 = d and x 2 = −d, then 2qx s = 2q x d, thus the modification appears in the q x -direction only. In the other directions, q y and q z , the single source result (42) is returned.
The appearance of the zero points is to a large extent an artifact of the used very simplistic two source model. In case of other additional sources these zero points would disappear. Nevertheless, this feature illustrates that the correlation function can be more complex than a set of Gaussians of the momentum difference q in different directions or at different rapidities.

B. Two moving sources
We study the system the same way as before, but now the two sources are moving in opposite directions, so that For now we also assume that FO happens at a t =const. FO hypersurface, so dσ µ = (1, 0, 0, 0) and so t 1 = t 2 .
If we have several sources then the source function in Jüttner approximation is while the J function is where x s is the 4-position of the center of source s, and the spatial integrals run separately for each of the identical sources, i.e. we assume fluid cells with identical density profiles, but with different densities, n s , velocities, u s and temperatures, T s .
The spatial integral for one source is the same as for a single source. Thus, where the factor exp(iq 0 x 0 s ) can be dropped if the FO time distribution is simultaneous for the two sources, because then x 0 s = 0. This returns Eq. (50) if u µ s = (1, 0, 0, 0). Now we can divide the two particle correlation with the square of the single particle distribution Ts Consequently, if the two sources have the same parameters, just opposite locations with respect to the center, and opposite velocities, then the correlation function is If we have two sources placed at x = ±d x , and with the velocity in the ±z-direction, ±v z , then the correlation function is for the different directions becomes: (58) 5. (color online) The correlation functions, C(k, q), for two moving sources where the displacement of the sources is in the x-direction, and the center-of-mass momentum, k, of emitted particles is in the x and y direction. The solid black line is for the momentum difference, qz, the dashed red line is for qy and dotted blue line is for qx. The radius of the sources is R = 1 fm, (same as in Fig. 3), the displacement is d = 1fm, and the source velocity is, γvz/Ts = 1.0 fm. This can be satisfied e.g. by us = 0.6 c and T = 0.12 GeV.
Therefore only the correlation functions in the k z , q xdirection and in the q z -directions are affected by the zdirected velocity of the source. In this direction, k z , unfortunately it is difficult to detect the two particle correlations. For the k x and k y -directions the q x -distribution is affected by the displacement of the two sources by ±d x . The q y -distribution is not effected by either the displacement or the source velocities. The correlation function for different source locations and velocities are similar. The cosine term appears in the same direction as the axis at which the sources are located and the hyperbolic cosine in the direction of the velocity. See Figs. 5 and 6. The zero points discussed for the two static sources at Eq. (51), appear in the distributions C(k x , q x ) and C(k y , q x ). These distributions do depend on the magnitude of the flow velocity, v z , but not on its direction! This arises from the fact that the detectors are assumed to be reached from both sides of the system with opposite velocities with equal probability.
Unfortunately the dominant direction of flow (see Fig.  8) is the beam direction (z−direction), where we have no possibility to place high acceptance detectors. At the same time the strongest effect of the flow appears in this direction.
The rotation in the reaction plane can also be characterized with another configuration of the two moving sources, when the displacement is in the z-direction while the flow velocities are pointing into the x-direction, so that the source at x 1 = d z has a negative velocity, −v x while the source at x 2 = −d z has a positive velocity, v x . (See Fig. 7.) The detailed description of the correlation functions from this configuration can be obtained in a straightforward way similarly to the previous case, see Eq. (61). In this case the flow has the most dominant effect in the k x -direction, which is accessible for detection. The x-directed flow, however, is more sensitively dependent on secondary effects, like the Kelvin-Helmholtz Instability [2]. In this configuration of the sources the magnitude of the flow velocity makes visible change in C(k, q), in the (k x , q x )-direction also, which is detectable by the usual detector configurations. Still the direction of the rotation does not appear in the observables with the approach presented here. This actually arises from the simplifying assumption, that the freeze out is happening instantly at a timelike hypersurface withσ µ = (1, 0, 0, 0), where particles from all sides of the system can reach each detector with the same probability. We will return to this problem after having discussed the more complex source configurations.
For these two-particle correlation measurements it is necessary to identify independently, event by event the global collective reaction plane azimuth, Ψ RP , experimentally and the corresponding event by event center of mass of the system (e.g. with the method [27]). Knowing these we can identify the k x -direction (and the k ydirection also.
In this section we derived a relatively simple formula for two sources with opposite positions and opposite velocities. These kind of systems were analysed earlier for radially expanding systems.
Recently due to the angular momentum in peripheral heavy ion collisions strong rotation [1] and turbulence (Kelvin-Helmholtz Instability) [2] were predicted in fluid (Color online) Typical orientation of the spatial axes in case of an ultra-relativistic heavy ion reaction shortly after the impact. In the configuration space the projectile and target appear to be flat due to the Lorentz contraction.
dynamical models arising from the symmetries, shear and vorticity of the initial state.
In the simple two source example shown in the previous section the two sources may describe a rotation if the sources are at a distance from the center in the xdirection, x 1 = (+d, 0, 0) and x 2 = (−d, 0, 0), while these have opposite velocities pointing into the z-direction, It is important to mention that to detect rotation the accurate identification of the reaction plane and its proper orientation is necessary. In the so called "cumulative" methods the reaction plane is identified but its projectile and target sides are not. This makes it impossible to detect directed flow, and odd components of the global collective flow. (All harmonic components of random fluctuations of course can be detected.) Furthermore, not only the reaction plane with proper direction but also the event by event center of mass (c.m.) should also be identified [27]. This hardly ever done! In both cases the use of zero degree calorimeters are provide an adequate tool as these are sensitive to the spectator residues.
The correlation function depends both on vectors k and q. To detect rotation the choices should be correlated correctly with the beam and the directed reaction plane as illustrated in Fig. 8. The positive x-axis points in the direction of the projectile, which moves in the positive direction along the z-axis.
In Eq. (57), in the above situation, ku s = γk z v z , qu s = γq z v z and qx s = q x d. Thus, the Correlation function, apart of the single cell source size, R, sensitivity, has a specific dependence on k z and q z , as well as on q x . Unfortunately it is difficult to measure the particle momenta in the z-direction as it coincides with the beam. The q x dependence would enable us to estimate the distance of the two sources.

C. Four Fluid Cell Sources
Four sources can be treated as a combination of two moving double source systems. We use the same parameters as under paragraph II B, where s 1 and s 2 will be the two different pairs of sources with different locations and velocities. See Fig. 9.
So we have If s 1 = s 2 then we recover Eq. (57) In the case of a rotating but symmetric system the displacements and velocities are of equal magnitude and are orthogonal to each other in the two pairs: x s1 ⊥ x s2 and u s1 ⊥ u s2 . Thus a simple sign change of the velocity for one of the pairs or both does not change the result, and so the rotation can be identified, but this evaluation does not provide sensitivity to the direction of the rotation. The reason is in the simplified freeze out assumption as we mentioned already at the end of paragraph II B.
If the two pairs are not completely identical, i.e. the magnitude of the characteristic quantities of the two source pairs are not equal then a sensitivity to the direction of the rotation may in principle occur. However, if we change the direction of the velocities of the two source pairs simultaneously (as it happens in changing the direction of rotation) the result still does not change.
Four Sources with Flow Circulation: Recent fluid dynamical studies indicate [1,2], that due to the initial shear and angular momentum the early fluid dynamical development has significant flow vorticity and circulation on the reaction plane. These were recently evaluated [4].
At the present LHC Pb+Pb collision energy in the mentioned fluid dynamical model calculation the maximum value of vorticity, ω, was found exceeding 3 c/fm , and the circulation after 6 fm/c flow development and expansion was still around 4-5 fm·c. This vorticity in the reaction plane was more than an order of magnitude bigger than in the transverse plane estimated from random fluctuations in ref. [28].
In this section we will look at the four source correlation function with similar circulation as in the above mentioned fluid dynamical model estimates in the reaction plane. See Fig. 9. We will simulate a circulation value Γ = 5f m · c. We use Eq. (62) where the center-ofmass momentum, k points in the x − direction.
Since the position and velocity are of the same value and because of symmetry the correlation functions C(k x , q x ) and C(k x , q z ) provide the same values. So we take the correlation function C(k x , q x ) and we have afterwards some simplifications. See Fig. 10. For large values of the center-of-mass momentum kx the correlation functions C(kx, qx) and C(kx, qz) will approach the correlation function C(kx, qy). The larger displacement and smaller rotation velocity leads to stronger deviation from the unaffected correlation function C(kx, qy).
For C(k x , q y ) we have the same result as we had for the two moving sources. Here the flow and displacement have no effect.
Let us look at comparisons for similar circulations and for similar displacements.
By comparing Figs. 11 and 12 we see that an increase in the displacement of the sources gives a increase in the apparent size of the system (narrower q−distribution. We also see that the measured size of the system increases with decreasing velocity. At the same time the shape of correlation functions are becoming less and less Gaussian as the flow velocities increase. At the same time the structure of the correlation function is also very different in different directions, which is not the case for spherical For large values of the center-of-mass momentum kx the correlation functions C(kx, qx) and C(kx, qz) will approach the correlation function C(kx, qy). Now for the dashed green lines with even higher velocity and smaller displacement, the deviation is significant and it is in the positive direction.
or linear expansion. This indicates that the rotating system contributes to essential non-Gaussian modifications, which can be seen directly in the correlation function, but they would become invisible if we would like to fit these data with a set of Gaussians. Earlier works studied the correlation function at different angles or pseudorapidities with Gaussian parametrizations [7,10], however, for rotating systems this is not the most sensitive way of presenting the results. Thus rotation can be detected even in "symmetric" few source systems where the emission is equally probable from all emitting sources. This emission scenario is less applicable to emission from heavy ion reactions where the absorption of particles in QGP is not negligible, and this affects the emission from the interior of a timelike (spacelike) FO layer, where the emission of earlier (deeper) emitted particles are quenched.
The correlation function is symmetric in all these cases as sources from opposite sides of the system contribute equally. Thus the correlation function is not sensitive to the direction of rotation.

III. ASYMMETRIC SOURCES
We have seen in the previous few source model examples that a highly symmetric source may result in highly symmetric correlation functions, however, this re- sults were not sensitive to the direction of the rotation, which seems to be unrealistic. We saw that this result is a consequence of the assumption that both of the members of a symmetric pair contribute equally to the correlation function even if one is at the side of the system facing the detector and the other is on the opposite side. The dense and hot nuclear matter or the Quark-gluon Plasma are strongly interacting, and for the most of the observed particle types the detection of a particle from the side of the system, -which is not facing the detector but points to the opposite direction, -is significantly less probable. The reason is partly in the diverging velocities during the expansion and partly to the lower emission probability from earlier (deeper) layers of the source from the external edge of the timelike (or spacelike) FO layer. This feature is recognized for a long time and discussed in detail by now. This influences the particle emission (or freeze out (FO)) process and modifies the post FO particle distribution. This topic has an extended literature, and this feature destructs the symmetry of emission of from source pairs at the opposite sides of the system [10,14,[17][18][19][29][30][31][32][33][34][35][36][37][38].
For the study of realistic systems where the emission is dominated by the side of the system, which is facing the detector, we cannot use the assumption of the symmetry among pairs or groups of the sources from opposite sides of the system. Even if the FO layer has a time-like normal direction,σ µ the (k µσµ ) factor yields a substantial emission difference between the opposite sides of the system. Now we want to demonstrate this effect on few source examples, and we will demonstrate the consequences of the non-symmetric emission.

A. The Emission Probability
It was first recognized that the freeze out with the Cooper-Fry description [14], may lead to negative contributions for particles, which move towards the center of the system and not in the direction out, towards the detectors. The first proposal to remedy this problem came from Bugaev [30], which led to the introduction of an improved post freeze out distribution in the Cooper-Frye description, first with the Cut-Jüttner distribution [30,32] and then by the Cancelling-Jüttner distribution [36].
Subsequently it was realized that for the realistic treatment of the freeze out process in transport theory one has to modify the Boltzmann transport equation by replacing the local molecular chaos assumption with a non-local one, where the point of origin is also included in the phase space distributions of the colliding particles. This led to the Modified Boltzmann Transport equation (MBT), and also the necessity to introduce an escape probability, P esc was pointed out.
The escape probability was then introduced and analysed in a series of publications [17][18][19]38], in transport theoretical approaches. It was pointed out that even if the pre FO distribution is a locally equilibrated isotropic distribution, the freeze out process and the escape probability will provide a nonisotropic distribution which eliminates the earlier observed problems. This developing anisotropy in the freeze out process occurs for freeze out both in space-like and time-like directions.
The escape probability introduced in the works [17][18][19]38], for a space-time surface layer of the system of thickness L, pointing in the four directionσ µ was given at a point x µ inside the freeze out layer as where p µ is the momentum of the escaping particle, u µ (x) is the local flow velocity and s = x µσ µ is the distance of the emission point from the inside boundary of the layer. The first multiplicative term describes higher emission probability to the particles, which are emitted closer to the outside boundary of the layer, the second multiplicative term describes the higher emission probability for the particles, which move in the normal direction of the surface, because these should cross less material in the layer. The last term secures that only those particles can escape, which move outwards through the layer.
The last two momentum dependent factors are important in transport theoretical models, to determine the shape of the post FO momentum distribution, e.g. [36], which would replace the Jüttner distribution. This shape modification happens to the single and two particle distributions equally, and it acts in all emission directions, k, equally, so this effect is secondary from the point of view of the flow velocity dependence of the correlation function.
In order to describe the complete freeze out process for a reaction the system had to be surrounded with a freeze out layer in the space-time, and the phase space distribution of the escaping, frozen out particles can be obtained by integrating over the whole 4-volume of the freeze out layer the local (usually isotropic) phase space distribution with the escape probability P esc (x). This procedure would then play the role of function G(x) in the source function in Eq. (10) instead of the simplified assumptions, as e.g. in Eq. (14).
The correlation function, C(k, q) is always measured in a given direction of the detector, k. Obviously only those particles can reach the detector, which satisfy k µσ µ > 0. Thus in the calculation of C(k, q) for a givenk-direction we can exclude the parts of the freeze out layer where k µσ µ < 0 (see Eq. (10) of ref. [10] or ref. [30]). For time-like FO a simplest approximation for the emission possibility is P esc (x) ∝ k µ u µ (x) [20].
In a model calculation we therefore have to define the freeze out layer also, this realistically should not include the whole space-time volume of the reaction. In case of calculating C(k, q) for a givenk we should select the relevant part of the freeze out layer, which may contribute to emission in the k direction. This should be a layer of 2-3 m.f.p facing the detector at the directionk. This can eliminate the symmetric pairs of fluid cells in the previous calculations of the correlation function, even if the emission normal is timelike, because the FO particle from an earlier emission point in the ST has to propagate through the plasma for some finite time, with considerable quenching.
Therefore in the following models we should apply the escape probability and we should define ak-dependent freeze out layer also! The most simple approximation is to select an emission layer from the system for a givenkdirection with uniform emission probability from within this layer. The next to most simple approximation is to introduce an emission probability within the layer, increasing towards the outside boundary of the layer. (Here it is important to mention that the spatial emission probability should be sufficiently smooth, so that one fluid cell and its contribution to C(k, q), should not be effected by this emission probability.
When we have up to 4 sources we can always addkdependent emission weights to these sources. This still would qualitatively change the outcome. As we discuss here up to four sources only a detailed formal evaluation of the emission probability would be an exaggerated approach, by defining more parameters than the outcome, so we just define the weights themselves here. In a full 3+1D fluid dynamical model with 100000+ fluid cells of course we have to apply a realistic and general evaluation of emission probability for every point of the ST.

B. Emission probabilities for few sources
Two sources: The previous discussion included two sources (i) in the beam-, z−direction and (ii) in the trans-verse direction in the reaction plane, x−direction. In case (i) the emission could be different from the two sources if the detector is in the z− direction, which is difficult to achieve, so we do not have to discuss this possibility. In configuration (ii) the observation can be in different k-directions. Ifk points into the ±y−direction, then the probabilities must be identical so emission probabilities do not lead to any change.
Ifk points into the ±x−direction, then one of the sources is closer to the detector and may shadow the more distant one. Thus, we can just introduce two positive weight factors so that w c is the weight for the cells closer to the detector and w s is for the cells which are far from the detector measuring the average momentum k. These weights are the same for the calculation of the nominator and denominator of the correlation function, so their normalization does not influence the correlation function.
As not all emitted particles reach a given detector the normalization is also dependent on the direction of the detector. Thus, we evaluate the correlation function this way. This immediately changes the earlier result (58), because it breaks the symmetry between the two sources. We can simply repeat the calculation for two moving sources in section II B, modifying the derivation of Eq. (56) and obtain the general result Note that this result is valid for the case whenk points to the +x direction, because the weights depend on this and w c > w s . See Fig. 13. The fact that the emission from the source, which is closer to the detector is stronger makes the direction of the flow detectable.
If we introduce the notation w c = 1 + and w s = 1− , the deviation from the symmetric result will become apparent If = 0 we have the symmetric situation where both sources have equal contribution, the asymmetric terms vanish, and the result becomes to be symmetric for the change of the direction of the flow velocity. If reaches its maximal value, = 1 the contribution of the far side source is eliminated (w s = 0, w c = 2), and only the single nearby source contributes to the correlation function. In this case the asymmetric term in the nominator vanishes, the remaining terms in the nominator and denominator are equal, and we recover the single static source result.
This result has terms, which change sign if the flow velocity, u s changes sign. The result is valid only if the detector is in thek = (1, 0, 0) direction. For this direction, however, if the flow velocity points in the z-direction, i.e. orthogonal to k the asymmetric term does not provide any contribution, so it will not show up in C(k x , q). To circumvent this problem we should study detector directions, which do not coincide with the primary axes of the given event (where x is the direction of the impact parameter vector, b, pointing to the projectile; y is the other transverse direction; and z is the direction of the projectile beam).
Correlation in Tilted Directions: The form of the correlation function is the same if k is in the same plane, the reaction plane, but it has a z component also, i.e. k = (k x , 0, ±k z ). This is possible for all LHC heavy ion experiments, ATLAS, CMS and even ALICE, where the longitudinal acceptance range of the TPC (∆η < 0.8) is the smallest. See Fig. 14.
Earlier for spherically or longitudinally expanding systems the dependence of the correlation function on the tilt angle or width parameter was analysed in detail in ref. [10]. We do not go into similar fine details, just demonstrate the possibilities for an arbitrary configuration.
Depending on the detector acceptance we should chose a detector direction where |k z | is as big as the detector acceptance allows it. For this configuration the form of the correlation function is the same as (66) with keeping the different weights, w c , w s or so that the forward shifted and backward shifted directions have the same weights. These weights are not specified up to now anyway.
For detection of the correlation function we have to introduce here the usual, k-dependent coordinate system to classify the direction of q. Thus if where a 2 + b 2 = 1, see Fig. 14. Then the difference vector, q, can be measured in the directionŝ q out = (a, 0, ±b), q x = a|q|, q z = ±b|q| q side = (0, 1, 0), q y = |q| q long = (∓b, 0, a), q x = ∓b|q|, q z = a|q|. (69) This leads to the following correlation functions Ts Ts Ts Ts Although, it seems that C(k (±) , q out ) and C(k (±) , q long ) are the same, this is in fact not the case, because the values of the components of the different types of k and q are not the same as described in Eqs. (68,69). In all cases, the out-, side-and long-q = |q|. We will also use the notation k = |k| and γv x = u x , γv y = u y , γv z = u z , so that u s = (u x , u y , u z ). For example for the out component the difference of the forward and backward shifted correlation functions is As Eq. (71) and Fig. 15 show, the Differential Correlation Function (DCF), ∆C(k ± , q out ), is sensitive to the speed and direction of the rotation, and it is also sensitive to the amount of the tilt in the directions of the detection, regulated here by the parameters a and b. ∆C(k ± , q out ) tends to zero both if q → 0 and if q → ∞. The structure of ∆C(k ± , q out ) is determined by the cosh (u z b q/T s ) cos (ad x q) product. If in both arguments the coefficients of q, u z b/T s and ad x are positive, smaller than one, and u z b/T s ≤ ad x , then the DCF is positive. If the coefficient ad x exceeds one the cos function changes sign at high q out values (e.g. above q = 1 − 2fm −1 ), and the DCF becomes negative at high q out values. Note that the ratio of the two coefficients is influenced be the choice of the tilting angle, i.e. by the parameters a and b.
If the parameter ad x remains constant, about 1 fm, and then when u z b/T s becomes larger (than one) the Differential Correlation Function becomes negative at small q values.
If the parameter u z b/T s remains constant, and about 1 fm, then when ad x becomes less (than one) the DCF becomes negative at small q values.
In case if the detector has a narrow pseudorapidity acceptance, then k ± is close to k x , i.e. b a and then the weights are maximal for the source in the x−direction, as indicated in Fig. 16.
If we change the direction of rotation to the opposite the Differential Correlation Function changes sign due to the sinh function in the nominator. In this configuration with the change of the tilt of the detector directions we can adjust the DCF, to the threshold value where the ∆C(k ± , q out ), is still positive, which provides a sensitive estimate for the rotation velocity at Freeze Out. This very sensitive behaviour is rather special and it appears in this special two source model this way. With an increased resolution and with more source elements this strong and specific structure will be smoothed out to some extent. The sinh(2u z bk/T s ) term changes sign in the nominator when u z changes sign the difference of the two correlation functions, ∆C(k ± , q out ) changes sign also because all other terms are symmetric to the sign change of the velocity.
This is an important observation as we can detect the direction and magnitude of the rotation in the reaction plane. This difference is also increasing with the longitudinal shift, b, of the average momentum vector, k, so that detectors with larger pseudorapidity acceptance can detect the rotation better.
In order to perform this measurement, one has to determine the global reaction plane (e.g. from spectator residues in the ZDCs), and determine the projectile side of this plane. Furthermore the event by event center of mass should also be identified (using e.g. the method shown in ref. [27]). This will be the positive x-direction. Then the correlation function can be measured for four different k-directions in the global reaction plane. These four directions are shifted forward and backward from the center of mass symmetrically on the projectile side, and there should be a symmetric pair of detection points in the target side of the reaction plane too.
The k directions opposite to each other across the c.m. point give the same result, while the difference, ∆C(k ± , q out ), between the Forward (F) and Backward (B) shifted contributions will characterize the speed and direction of the rotation. This symmetry can be used to eliminate the contribution from eventual random fluctuations. The observed F/B asymmetry depends on the parameters , v z and d x , these can be estimated by measuring the correlation functions at all possible moments k. Fig. 15 indicates that the differential correlation function has a larger amplitude for smaller k values, and the zero points are sensitively dependent on the rotation velocity.
The zero points come from the term The values of k are: for the solid black lines 0.10, 1.00, 3.00, and 5.00 fm −1 and for the dashed red lines 0.10, 1.00, 3.00, and 5.00 fm −1 . For a detector with narrow pseudorapidity acceptance we could use a value a > b. For this figure we have used a = 4b. We can compare this with Fig. 17, and we see that for q > 1.5 fm −1 they are very similar. The DCF has a smaller amplitude in the positive direction for values around q < 1.0 fm −1 , but a larger amplitude in the negative direction for values around 0.75 fm −1 < q < 1.5 fm −1 . We also see that the DCF is not close to zero for values around 2.0 fm −1 < k < 5.0 fm −1 .
and it is not dependent on k, so for the values used in Fig. 15 and q = x fm −1 we have Since the cosine term must be positive, there are no zero points for π/ √ 2 ≤ x ≤ 3π/ √ 2 or 2.23 ≤ x ≤ 6.66.
For x > 3π/ √ 2 the exp(−R 2 x 2 ) term will be small and there would be no correlation difference. So we will look at the values 0 < x < π/ √ 2 If v z is smaller than 0.57c then there will be no zero point for 0 < x < π/ √ 2. For v z larger than 0.58c there will be 1 zero point for 0 < x < π/ √ 2. The zero point for a given velocity can be found by solving equation (73)

C. Emission from four sources
With four sources we can illustrate the possibilities of differential HBT method studies in different directions. The correlation functions can be calculated in general for four sources and two detector positions. This can then be applied to different detector configurations.
The out component of the four source correlation function with weight factors ω a , ω b , ω c , ω d is given by similar structure in the two source model and the four source model in a resembling configuration. Fig. 20 has the same shape as Fig. 17, but the amplitude is larger.
For a set of large number of sources, forming a system with close to perfect rotational symmetry, a single correlation function would not depend on the (polar) angle of the detection, and the DCF would vanish. Thus, the DHBT method would not be applicable for highly symmetric systems, like for a rotating star observed from within the plane of the rotation. At the same time for a rotating binary star system the DHBT method would work. Also the weighting of the sources should be different: If the observer is in the plane of rotation, the distant star is shadowed by the front one at some periods, just like emission from a highly opaque plasma (evidenced by jet quenching). If the observer is slightly out of the plane of rotation then the two stars are visible all the time and then the (time dependent) correlation function would change between the configurations of Figs. 4 and 7. This also illustrates the role of symmetric and asymmetric weightings.
The rotating and expanding final state of a relativistic heavy ion reaction is of course does not look like a perfect wheel, so the four source model is a more adequate approximation than a wheel would be.

IV. CONCLUSIONS
In this work we attempted to study the possibility of detecting and evaluation the rotation of a source by the specific use of the Hanbury-Brown and Twiss method for rotating systems. Our primary interest was the application for ultra-relativistic heavy ins where in peripheral collisions at ultra-relativistic energies the system can gain large angular momentum. Nevertheless, some of the conclusions can be applied to macroscopic systems also, like for past rotating stars.
We selected one of the several methods to evaluate two particle correlations, which was suitable to study collective fluid systems with significant and well defined internal fluid dynamical motion. The obtained standard correlation functions were showing the consequences of the flow, but for highly symmetric sources the correlation functions gave symmetric results, which were invariant for the change of the direction of rotation.
It turned out that it is important to take into account that the particles reaching the detector cannot reach it with equal probability from the near side and the far side of the emitting object. With this fact considered we could obtain correlation functions, which reflected the properties and also the direction of the flow. These results can be used rather generally.
The obtained results have shown that the correlation function is most sensitive to the rotation if it is measured in the beam direction (or close to it). This, unfortunately, is not possible in most heavy ion accelerator experiments, so we introduced and investigated a Differential Hanbury Brown and Twiss method, which made it possible to trace down the rotation in relativistic heavy ion collisions by measuring the correlation functions in the reaction plane at nearly transverse angles to the beam direction. The method is promising and can be performed in most heavy ion experiments without difficulties, as well as it can be implemented in different reaction models, like fluid dynamical models, microscopic transport models and hybrid models. In full scale theoretical models, the emission probabilities from the FO layer have to be considered. From the general formulas derived in the beginning of the paper apparently these dependencies can be factorized.
To complement these studies we also applied the method to a high resolution, 3+1D, computational fluid dynamics model [39], which was used earlier to predict roation, KHI, flow vorticity, and polarization [1,2,4,40]. The result shows that the method can detect rotation, while the effects of irregular shape, sperical flow, and specific flow patterns, require a more extended analysis to separate all these effects form one another.