Detecting a First-Order Transition in the QCD Phase Diagram with Baryon-Baryon Correlations

We suggest baryon-baryon correlations as an experimentally accessible signature for a first-order phase transition between a baryon-rich phase, like quarkyonic, and a baryon-suppressed hadronic phase in the QCD phase diagram. We examine the consequences of baryon-rich bubble formation in an expanding medium and show how the two-particle correlations vary in the transverse and longitudinal direction depending on the strength of the radial flow, the bubble temperature, and the time when the baryons are emitted.


Introduction and Motivation
Exploring the Quantum-Chromodynamics (QCD) phase diagram has been the focus of significant experimental and theoretical research in the last few decades. In particular, the effect of heating and compressing nuclear matter on confinement and chiral symmetry breaking has been studied, new phases of matter have been suggested, and the existence of a critical endpoint. predicted [1]. Relativistic heavy-ion experiments also aim to explore the temperature-baryon density plane and to provide evidence for a critical endpoint or a first order phase transition. Therefore, it is important to define reliable signatures for identifying a critical point, a new phase, or a first-order transition line.
In attemtping to understand QCD, it has often proven instructional to explore QCD-like theories in which some of the parameters of the Lagrangian (quark mass m q , number of colors N c , number of flavors N f ) are varied [3]. McLerran and Pisarski used such a study to argue that at large N c and sufficiently small temperature to baryon chemical potential ratio T /µ B the phase diagram exhibits a strongly first-order phase transition to what they termed quarkyonic matter [2]. For large N c and all N f , baryon number has been suggested as an order parameter for the transition [4]: In large N c the baryon mass is of the order M B ∼ N c Λ QCD . For moderate temperatures, T ∼ Λ QCD , the expectation value of baryon number is n B ∼ exp (µ B /T − M B /T ) ∼ e −Nc . This is negligibly small (identically zero at large N c ), and it remains that way Email addresses: amocsy@pratt.edu (Ágnes Mócsy), prsorensen@bnl.gov ( Paul Sorensen) as long as µ B ≤ M B . For larger µ B the baryon number density becomes non-zero. In the deconfined quarkgluon plasma phase there are no baryon masses, so that there is no baryon-number suppression [4]. Figure 1 shows the phase diagram with a hadronic (mesonic) phase and a high energy density quarkyonic phase. At small temperatures the transition from the mesonic to quarkyonic phase should be first order [4]. Although large N c is not QCD, it is not excluded that the QCD phase diagram bears some resemblance to the large N c diagram. For realistic N c and N f the boundary of these regions may be crossovers. The question we asked ourselves was, if there is a firstorder transition between the quarkyonic and mesonic phase, what are the observable consequences? The arguments laid down in the rest of the paper are, however, more general, because they only require the existence of a baryonrich and a baryon-poor phase. Phase conversion between these can happen either by bubble nucleation or by spinodal decomposition. Here we explore the phenomenolog-ical consequences of bubble nucleation in QCD matter. Spinodal decomposition has been addressed in [5]. It has been suggested by Voloshin [6] that correlations in coordinatespace can be transferred to momentum-space through radial flow. Baryon number fluctuations were discussed in [7], but the effects of radial expansion and azimuthal correlations were not considered in that work. We focus on the transfer of spatial correlations between baryons confined to baryon-rich regions, into baryon-baryon correlations in momentum space via radial flow.

Baryon-rich Bubbles
In this work we consider a bubble nucleation scenario, similar to the one proposed for strangelet condensation in the early universe [9]. Imagine an expanding system in which the initial temperature is higher than that of the transition between quarkyonic and mesonic matter. Initially there is non-zero baryon number density. As the system cools a small amount below the phase transition temperature, bubbles of stable mesonic phase start to nucleate. It has been shown in a dynamical calculation by Paech and Dumitru [8] that baryon inhomogeneities do form. The mesonic bubbles expand rapidly, reheating the system back to the transition temperature. The system then is in a mixed phase, which at first is mostly quarkyonic with small droplets of hadron gas. Then these mesonic droplets increase in size. At the time when roughly half the matter is in the quarkyonic phase one can think of bubbles of quarkyonic matter evaporating in the hadron gas, as shown in Figure 2. The ratio of baryon densities of the hadronic phase and that of the quarkyonic phase n H /n Qy ∼ e −Nc is tiny (identically zero in large N c ). Because of this exponential suppression of the baryon density in the mesonic phase, all of the baryon number remains trapped inside the bubbles of quarkyonic phase, and at later times it can concentrate in the evaporating droplets of quarkyonic matter. If bubbles of quarkyonic matter have baryon number larger than one, large variation will exist in the spatial baryon number-density. In the following we study how inhomogeneous the corresponding distribution in momentum space will be. The longitudinal position inhomo-geneities will be translated into inhomogeneities in rapidity because of the familiar correlation between momentumand coordinate-space rapidity. In the transverse space there should be a similar correlation due to transverse flow [6].
Baryon correlations due to transverse flow are not likely to be washed out by thermal noise. The momentum imparted to a hadron due to flow depends on the mass and flow velocity p f low ∼ M v f low while the thermal noise depends on the mass and temperature p noise ∼ √ M T . Since p f low /p noise ∼ v f low M/T and since the baryon mass is at least a factor of five larger than the temperature, we expect space inhomogeneities for baryon number to be observable as baryon-baryon correlations in momentum space. We investigate these by constructing a Bubble Monte-Carlo Blast-Wave Model similar to Ref. [10].

Bubble Monte-Carlo Blast-Wave Model
Our model has two components: One is the Bubble Monte-Carlo, in which the spatial distribution of baryon number and bubbles is generated; The other is a Monte-Carlo Blast-Wave which provides the momentum distribution of baryons. We study how spatial correlations in the source function are translated into momentum space and examine the effect of baryonic bubble formation on the two-particle correlation functions.
In our Bubble Monte-Carlo, we first make an assumption for the average number of baryons per bubble and the distribution of the baryon number per bubble. We generate baryon number N B for each event according to a gaussian distribution in a finite window of rapidity within which we know the average baryon number,N B . We then distribute the baryons into bubbles with the number of baryons per bubble n B taking the form We sample this baryon-per-bubble distribution to generate bubbles until we have generated the total number of baryons within the given window. This fixes the number of bubbles. Our treatment of N B and n B is somewhat arbitrary, but in this analysis rather than focusing on number fluctuations sensitive to this implementation, we will study the shape of the two particle correlations in rapidity and azimuth (φ). The bubbles are then distributed in space with two different transverse density profiles: The Center Weighted Source profile, assumes a bubble density given by a Woods-Saxon distribution integrated along the longitudinal z-direction. This distribution yields a density largest at the center of the system. The Surface Weighted Source gives a preference to bubble formation near a surface, motivated by the picture of an expanding system where bubble nucleation happens preferentially at a freezeout surface, near the edge of the system. Transverse projections of both source functions are shown in Figure 3. For both cases, we take the distribution to be flat in z. Baryons are then distributed inside the bubbles using again a Woods-Saxon distribution of the form The bubble radius is taken to be r b = 0.6 fm(n B ) 1/3 , so that the volume per baryon remains constant, the radius of the system R = 12 fm, and the width of the distribution a = r b /8 . In this way we obtain the positions (x, y, z) for each baryon. This discrete set of coordinates is later taken as the source function Ω. Figure 4 shows the baryon-and  bubble-distribution in the transverse plane of one event.
Each bubble is then assigned a radial boost velocity as is done in the blast-wave model. To construct a blast-wave Monte-Carlo we calculate the emission function, S. The emission function is the probability of emitting a baryon of transverse mass m T = p 2 T + m 2 B , transverse momenta p T , and rapidity Y B , at an azimuthal angle φ B from a boosted source of Ω at temperature T . We use the discrete set of points generated with our bubble Monte Carlo as our source function Ω. Following [11], we write the emission function as where: Here we will consider central collisions, which implies the direction of the boost β is identical to the spatial azimuthal angle φ b = φ s . Generalization to finite impact parameter presents no major technical difficulties. We allow the bubbles to decay isotropically in their flow boosted rest frame. Assumption of boost invariance manifests in the equality of the longitudinal flow rapidity to the spacetime rapidity η. For simplicity, we assume that the emission is instantaneous and replace the Gaussian in (4) with δ(τ − τ 0 ). We are left to specify the transverse rapidity boost β = r[β 0 + β 2 cos(2φ b )]. The parameter β 0 is the transverse rapidity in the outward φ s direction. We take β 2 to be zero. We will investigate different values for β 0 taking guidance from the fits in [11], which reproduces well the experimental data.

Baryon-Baryon Correlation Results
We investigate how baryon-baryon correlations in terms of relative angle (azimuth ∆φ and rapidity ∆y) vary with the freeze-out temperature T , average flow velocity β , number of baryons per event N B , number of baryons per bubble n B , time at which bubbles decay τ decay , and source function (Surface vs Center Weighted). We present the correlations in terms of a per particle correlation measure [12] ∆ρ Here ρ is the pair density and ρ ref the product of single particle densities. This normalization of the correlation measure means that for all other parameters kept fixed ∆ρ/ √ ρ ref is independent of the total number of baryons per event N B [12]. This is because normalizing by √ ρ ref instead of ρ ref removes the dilution factor of one over the number of particles analyzed. This was confirmed to be true in our Monte-Carlo (the variation of r b with N B does not seem to influence the correlation function). We generate our reference spectrum by sampling the single particle dN/dφ and dN/dy distributions generated in our Monte-Carlo enough times so that the same number of pairs are created in our reference distribution as in our real distribution. This ensures that ∆ρ/ √ ρ ref will integrate to zero. For this reason we are insensitive to a global offset and focus instead on the shape of the correlation, which is easier to measure experimentally. The variable ∆ρ/ √ ρ ref as a function of the azimuthal angle between baryons from the decayed bubble for a variety of different model inputs is shown in Figure 5. Bubble nucleation and radial flow leads to small angle correlations. We find that all correlation shapes can be reasonably well parameterized by a Gaussian centered at ∆φ = 0 and a constant offset. The correlation can therefore be summarized in terms of the Gaussian amplitude A and root mean squared width σ. By construction the quantity ∆ρ/ √ ρ ref is independent of the total number of baryons per event [12]. We confirmed this in our simulation, and then fixed N B = 100.
The top left panel in Figure 5 shows results with a Surface Source and T = 100 MeV, β = 0.75, and n B = 2, 4, 8. The width of the correlation is narrowest for n B = 2, since this corresponds to the smallest bubble size, but the variation of σ with n B is very weak (only changing from 0.487 for n B = 2 to 0.499 for n B = 8). The dominant effect of varying n B is on the amplitude of the correlation, which changes from 1.3 to 8.9, respectively. This dependence can be understood in the following way: The number of pairs that can be formed from one bubble is n B (n B − 1)/2. So the number of correlated pairs in an event is We note here that this amplitude only applies to the case that the correlation is formed only with baryons produced from bubbles. Any background contribution (for instance from baryons produced outside of bubbles) will dilute the signal.
In the top right panel in Figure 5 we show the dependence of the correlation on the temperature of the bubbles emitting the baryons, for a fixed β = 0.75 and n B = 8. This temperature is not necessarily the temperature at which the bulk of the hadrons freeze-out, but could be closer to the value of the critical temperature T c , where the bubbles actually form. We vary this T from 120 MeV to 180 MeV. It is clear that the shape of the correlation changes appreciably with temperature: higher T lead to a broader correlation function. In this temperature range and for this flow velocity the r.m.s. width is given by the fit σ = 0.694 − 0.00639T + 0.0000435T 2 . Since the random thermal boost competes with the flow boost, a larger bubble decay temperature washes out the correlation. Therefore, if the first order phase transition is at a high T c , correlations from bubbles will be harder to detect.
Since we expect a competing effect of the radial flow and the temperature, we now look at how flow influences the correlation. We illustrate this for T = 100 MeV, n B = 8, and β = 0.75, 0.55, 0.40 on the bottom left panel of Figure 5. The correlation is strongly dependent of the mean flow velocity and is the strongest for bubbles with larger boost. Remember, that flow is required to convert coordinate-space correlations to momentum space. So again, if the phase transition happens early, before flow builds up, bubbles will be hard to detect.
The bottom right panel of Figure 5 displays the flow effects at T = 100 and n B = 8, with either a Centeror Surface Weighted source function. The mean radial position of the source is larger for the Surface Weighted source: 10.1 fm compared to 7.4 fm. This leads to a larger average boost velocity: β = 0.75 compared to β = 0.55 , respectively. A larger β leads to a more pronounced correlation structure with the Gaussian narrowing from σ = 0.61 to σ = 0.50. We also compared a Surface Weighted source with a reduced β 0 such that the β is matched to the Center Weighted case. The Surface Weighted correlation with β = 0.55 is only slightly narrower than the Center Weighted correlation with the same β . The primary difference between Surface-and Center Weighting is caused by the change in β , but some deviation occurs due to the difference in geometry. Our interpretation is that the Surface Weighted source has a smaller range of β values, while the Center Weighted source has a broad distribution of β values. In this case, even for the same β value, the Center Weighted source gives a convolution of a broad range of β values which leads to a broader correlation function. Also, a bubble near the center covers a larger solid angle. We investigated three β values for each geometry. The widths within the ranges studied can be parameterized as follows: σ = 1.317 − 1.103β for the Center Weighted source (see lower right panel) and σ = 1.357 − 1.367β for the Surface Weighted source (see lower left panel). These parameterizations of our model are valid for 0.29 < β < 0.75.  Figure 6: Proton-proton correlation versus azimuthal angle and rapidity difference between the protons for bubble decay at times times τ decay = 1 fm (left) and τ decay = 9 fm (right).
We also looked at correlations between baryons in the rapidity y direction. Figure 6 shows the proton-proton correlation as a function of rapidity difference and azimuthal angle between the protons at different proper times of bubble decay, i.e. τ decay = 1 fm (left panel) and τ decay = 9 fm (right panel). We found that if the time at which bubbles decay τ decay is long then narrow correlations can build up both in ∆φ and ∆y directions. These correlations are increased in amplitude and narrowed in width. We conclude that if bubble decay happens near the beginning of the system evolution (small τ decay ) implying a high decay temperature and low flow, baryon-baryon correlations will be broad in both directions. For bubble decay near freezeout (large τ decay ), on the other hand, strong correlations with narrow widths in both directions can be expected. Our fit for the change in the width in rapidity direction with time is σ y = 0.37 + 0.47τ −1.58 .
Our studies of ∆ρ/ √ ρ ref show that, the amplitude and width of the correlation varies substantially within this bubble Monte-Carlo blast-wave model depending on the parameters assumed: We have shown in Figure 5 That correlation caused by baryonic bubbles emitted from a boosted source can change by a factor of two even with reasonable values for the average boost velocity and freeze-out temperature. The signature of a first order phase transition will be strong if the transition happens near the freeze out since then the flow will be largest and the temperature lowest. We note, that the decay time is linked to the temperature and flow, and with a dynamical model one could make this connection explicit.
We also consider how this correlation can affect measurements of elliptic flow v 2 [13] and it's fluctuations [14]. One method for estimating v 2 = cos(2φ−Ψ) when the reaction plane angle Ψ is not known is to calculate cos(2∆φ) from two-particle correlations. This quantity depends on v 2 2 , σ 2 v2 , and other correlations not related to the reaction plane, called non-flow δ 2 . Since we only consider central collisions in this study v 2 is zero. The calculation of cos(2∆φ) is therefore related to either non-flow correlations or v 2 fluctuations. In either case, the correlations from our model will lead to a deviation between the 2 nd and 4 th order cumulants [13] such that v 2 {2} 2 − v 2 {4} 2 is non-zero. Since v 2 is zero the expected deviation between the cumulants is given by This quantity depends strongly on the variables explored in the discussion above, as well as on the number of uncorrelated particles included in the analysis (misidentifed mesons for example). Larger multiplicities will dilute the correlation leading to smaller differences between v 2 {2} 2 where C can be expected to be between 0.3 and 5.
Measurements of p-p and p-Λ correlations have been carried out at energies ranging from √ s ∼ 2 GeV to 200 GeV [15]. In these HBT analyses correlation functions are defined as C(k * ) = N pair,real /N pair,mixed . The variable k * = Q inv /2 is the relative momentum of the particles in the pair rest frame, N pair,real is the number of pairs in an event with relative momentum k * , and N pair,mixed is the number of pairs expected from random combinatorics and is constructed from mixed events. The region of interest to an HBT analysis is below k * ∼ 100 MeV. The N pair,mixed distribution is usually normalized to the N pair,real distribution at some value of k * above this region. Our Monte Carlo model yields a correlation C(k * ) well described by a Gaussian with an RMS width of 750 MeV. For bubbles with n B = 8 the amplitude of the Gaussian is 8%. If the N pair,mixed distribution is normalized to the region near k * = 500 MeV the correlation induced by bubble formation would show up as a smooth 1-2% variation of C(k * ) between 0 and 500 MeV. In HBT analyses, long-range non-femptoscopic correlations of this kind are treated as background and ignored or parameterized. We therefore conclude that although HBT measurements performed at lower beam energies do not report signs of a strong baryon-baryon correlation as described here, they also do not exclude their presence.

Summary and conclusion
We have shown that bubble nucleation at a first or phase transition coupled with radial flow in heavy-ion collisions can lead to detectable baryon-baryon correlations. We've mapped out the shape of these correlations in rapidity and azimuthal angle and discussed how they will manifest as non-flow or flow fluctuations by contributing to the difference between v 2 {2} and v 2 {4}. We explored how the correlations will depend on the strength of the radial flow, and the temperature of the system at the time when the baryon rich regions decay. We reported the variation of the azimuthal and longitudinal width with temperature, flow, geometry, and decay time. We find that if bubble decay happens late in the evolution, close to the freeze-out when the temperature is lower and the flow is larger, the correlations will be narrow in rapidity and azimuth. We conclude that these correlations would be easy to detect. If the bubbles decay early in the evolution, the flow will be weaker, the temperature will be higher, and the correlation will be wider in both directions and therefore less pronounced in the data. In this later scenario, there will be more hadronic rescattering which may wash out the signal entirely. The effect of the hadronic stage on these correlations is an important study which will be the focus of future more quantitative work.
The observation of baryon-baryon correlations as described in this work will be evidence for the existence of a first-order phase transtion. The lack of such a signal will indicate that either bubble nucleation at a first order phase transition did not occur, or that the transition was sufficiently seperated from freeze-out as to render the correlations unobservable. We propose therefore, that studies of baryon-baryon correlations can be used to answer whether a first-order phase tranition is present in the vicinity of the freeze-out curve in heavy-ion collisions. Future studies of the effect of the hadronic rescattering stage will allow us to more precisely specify the region within which a first order phase transition should be detectable. If the transtion to quarkyonic matter is first order and close enough to hadronic freeze-out, then we conclude that it can be detected through baryon-baryon correlations.