The unstable mode of the net-baryon density near the spinodal decomposition region

The exploration of the Quantum Chromodynamics (QCD) phase structure is one of the important frontiers in both theoretical and experimental research of strongly interacting matter. Simulations of the lattice QCD revealed a crossover transition from the quark-gluon plasma phase (QGP) to the hadronic phase at vanishing baryon chemical potential ($ \mu_B $) [1−4]. Extensive studies based on the effective theories of QCD conjectured a first-order phase transition at large $ \mu_B $, together with a critical point between them [5−8].

The characteristic feature of the phase transition is the large fluctuations and divergence of variables. After being created in relativistic heavy-ion collisions, the QGP fireball is driven across the QCD phase transition region, and several final observables potentially preserve the non-trivial behavior induced by phase transition. Non-monotonic behavior of net-proton multiplicity fluctuations with respect to the colliding energy is believed to be one of the most important signatures of QCD phase transition in relativistic heavy-ion collisions [9, 10]. Beam Energy Scan (BES) program at RHIC has been dedicated to searching the QCD critical point, and the preliminary non-monotonic behavior of net-proton fluctuations has been observed in the first phase of the program (BES-I) [11, 12]. However, the statistic of net-proton fluctuations is insufficient to conclude the observation of non-monotonic behavior at BES-I, and is looking forward to the higher statistic measurement at the second phase of the BES program (BES-II).

In addition, the theoretical modeling requires further improvement for a better description of the fireball system near the QCD phase transition region. The QGP fireball created in relativistic heavy-ion experiments is a complex system, in which the critical fluctuations substantially deviate from the ones at an ideal and equilibrium context. In particular, the finite time effects induced by the expanding fireball have been revealed to significantly modify the critical fluctuations. For example, the magnitude of the fluctuations is suppressed [13, 14], the sign could be reversed [15], and the largest fluctuations not necessarily correspond to the trajectory closest to the critical point [16]. Consequently, several dynamical models have been built to incorporate the dynamics of the critical effects (see e.g., Refs. [17−22] for recent reviews). Typically, the dynamics of the critical fluctuations are regarded as an additional degree of freedom and couples with the conventional dynamical model (e.g., relativistic hydrodynamics), such as the dynamics of the order parameter field in non-equilibrium chiral hydrodynamics [23−25] and the slow mode in the hydro+ [26]. In particular, the dynamics of the conserved net-baryon density have been investigated [27−31] and it was found the non-monotonic behavior of multiplicity fluctuations with respect to the detector acceptance window [29−31].

On the other hand, the study on the dynamics of the fluctuations near the first-order phase transition side is relatively deficient [20], owing to its complexity. The characteristic property of the system near the first-order phase transition is the existence of the unstable mode, which may result in nucleation or spinodal decomposition. In relativistic heavy-ion collisions, the spinodal decomposition has been investigated within the framework of dissipative hydrodynamics using an approximate equation of state for the coexistence region [32−36]. In these works, the unstable mode in the coexistence region was studied, with the increasing density moments as the consequence. On the other hand, the acceptance dependence of net-proton fluctuations is feasible in heavy-ion experiments [37]. With the fluctuations within a specific detector acceptance window (e.g., $ \Delta y $), the collision energy dependence of the multiplicity fluctuations scans different regions of the QCD phase diagram. From another perspective, varying the acceptance window (e.g., $ \Delta y $) at specific collision energy $ \sqrt{s_{NN}} $ enables us to study the multiplicity fluctuations with increasing particle number involved. For example, Refs. [38, 39] study the multiplicity fluctuations of protons with increasing detector acceptance window (transverse momentum $ p_T $ range and rapidity range $ \Delta y $) near the critical point and find the acceptance dependence as a consequence of long-range correlation. However, the multiplicity fluctuation with respect to the acceptance at the first-order phase transition side has not been addressed yet. This work will focus on the consequence of the unstable mode in the spinodal decomposition region on the acceptance dependence of the multiplicity fluctuations. It is found that the unstable mode results in the correlation function dominating at a much harder mode comparing the ones outside the spinodal decomposition. This effect picks out the harder mode and shows oscillating behavior with increasing acceptance.

To describe the dynamics near the QCD phase transition, it requires the dynamics of the relevant degree of freedom. Relativistic hydrodynamics takes into account the dynamics of several slow degrees of freedom and achieves extraordinary success in relativistic heavy-ion collisions. Unstable mode analysis near the spinodal decomposition region has been performed [32−36] within the framework of relativistic hydrodynamics, with the equation of state constructed by performing a suitable spline between two idealized systems (either a gas of pions and interacting nucleons or a bag of gluons and quarks). In this framework, the thermal noise has not been taken into account, which is essential for the region near the critical point of the coexistence side. As there exist several relevant degrees of freedom in fluctuating hydrodynamics near the coexistence region, it is challenging for the comprehensive study both theoretically and numerically.

As an alternative framework to study the dynamics near the QCD phase transition, model B only evolves the conserved quantities near the critical point and is much easier to implement. Near the QCD phase transition, the order parameter field (chiral condensate $ \sigma\equiv \bar{q}q-\langle \bar{q}q\rangle $) is a slow mode and couples with other slow hydrodynamic modes, such as the conserved baryon number density $ n_B $, as well as densities of energy e and momentum $ \pi^i $ [40, 41]. It was also pointed out that the baryon density mode evolves much slower than the one of σ mode and plays a more important role near the phase transition [40]. This enables us to only study the dynamics of $ n_B $ as an approximation, which is the so-called model B. On the crossover side of the QCD critical point, several phenomenological studies with the conserved quantity in heavy-ion collisions have been implemented [27−31] and the multiplicity fluctuations with respect to the increasing acceptance have been investigated [29−31]. Therefore, this work focuses on the extension of this model to the system near the critical point on the coexistence side and theoretically studies the corresponding multiplicity fluctuations with increasing system size.

The coupling between order parameter field σ and $ n_B $ transmits the effects of phase transition to the baryon number density field and the conventional Landau-Ginzburg free energy can be taken as a good approximation of the effective potential of net-baryon density:

where K is the coefficient of surface tension. As one type of first-order phase transition, spinodal decomposition is an important mechanism for phase separation in many areas of physics. One of the characteristic properties of the spinodal decomposition is the corresponding effective potential exhibits a convex anomaly, $ v<0 $ and $ u>0 $ in this case [42]. The small deviation from uniformity in this convex anomaly region is amplified and eventually results in phase separation.

To study the phenomenological consequence of the spinodal decomposition in relativistic heavy-ion collisions, the coefficients in Eq.(1) are required to be specified from the QCD theory. The pressure P can be obtained from the effective potential (1) and therefore the coefficients in Eq.(1) can be expressed in terms of various orders of baryon number susceptibility, as the way in Refs. [30, 31]. Unfortunately, the surface tension K is essential near the first-order phase transition (we shall see later) and the knowledge of K (as well as other coefficients) is insufficient to make a conclusive comparison with experiments. This work focuses on the theoretical consequence of unstable spinodal decomposition on the multiplicity fluctuations without specifying parameters of strongly interacting matter. A more comprehensive and realistic analysis is left for future study.

The basic theory of spinodal decomposition is the Cahn-Hilliard equation [43], which has been widely studied in metallurgy. The dynamics of conserved net-baryon density $ n_B(t, x) $ at space-time $ (t, x) $ near the phase transition can be described [27, 28]:

where the noise term employed in this work is white noise and its moments have the relations

The coefficient of the second moment of $ \zeta(t, x) $ is determined by the fluctuation-dissipation relation. The notation $ \langle \cdots \rangle $ represents the averaging over the noise configuration. D denotes as the diffusion coefficient and T is temperature.

The diffusion equation in Eq.(2) with non-linear effective potential can not be solved analytically. Nevertheless, it still can be analyzed by assuming the small perturbations $ \tilde{n}_B $ near the minima of free energy $ n_{min, B} $:

and the evolution of $ \tilde{n}_B $ takes the form:

This equation of $ \tilde{n}_B $ has formal solution in Fourier k-space:

with the initial condition $ \tilde{n}(t_0, k) $ at $ t=t_0 $. The notation $ \tilde{D}(t_0, t, k)\equiv \int^t_{t_0}dt'D[(v+3un^2_{min, B})k^2+Kk^4] $ is used for convenience. One can see that the mean value of perturbation for the conserved quantity $ \langle \tilde{n}_B(t, k)\rangle $ decays exponentially from the initial fluctuations $ \langle \tilde{n}_B(t_0, k)\rangle $, without contribution from the noise, as can be seen from Eq.(3). In the limit of zero mode $ k=0 $, corresponding to the entire system, the fluctuations are constant as a function of time, which is the characteristic feature of the dynamics of the conserved quantity (model B).

The initial condition is required to obtain the two-point correlation function. It can be assumed that the system is thermalized and the correlation follows as [29]

With the initial condition of the correlation function and performing the averaging over the noise configuration, the correlation function in Fourier space can be obtained as

Typically, the system is assumed to be homogeneous, namely, the correlation function $ \langle \tilde{n}_B(t, x_1)\tilde{n}_B(t, x_2)\rangle $ only depends on the distance $ x_1-x_2 $, not the abstract value $ (x_1+x_2)/2 $. Therefore, $ \delta(k_1+k_2) $ appears in the correlation function.

As the case of perturbation $ \langle \tilde{n}_B(t, k)\rangle $, the correlation

decays exponentially from the initial value $ \chi(t_0) $. However, the second moment of noise in Eq. (2) is non-zero and contributes to the two-point correlation function $ C(t, k) $ (the second term of Eq.(8)). Because of the noise contribution, the correlation $ C(t, k) $ tends to develop a stronger correlation at larger mode k. Outside the spinodal region ($ v>0 $), $ \tilde{D}(t_0, t, k) $ always positive (the diffusion coefficient D, surface tension coefficient K and fourth order coupling constant u are typically positive), the correlation function decays rapidly for the hard mode k and dominates at small mode k. On the contrary, $ \tilde{D}(t_0, t, k) $ could be negative in the spinodal region ($ v<0 $) due to the convex anomaly, and the correlation function significantly enhances at hard mode k, if the surface tension K is neglected. The correlation function accumulates as a function of time and diverges eventually, namely becomes unstable. Therefore, the surface tension is essential for the stability of the system near the phase transition. The effects of surface tension have been addressed in Refs. [32−34, 36] and introduced as a gradient term in the equation of state to describe the finite range of interaction. In short, the surface tension K stabilizes the system, but the negative contribution of v in the spinodal region results in the largest $ C(t, k) $ at harder mode k, comparing outside spinodal region ($ v>0 $). For the illustrative purpose, the correlation function $ C(t=1, k) $ in k-space is shown in Fig. 1 with $ D=1, \chi(t_0)=\dfrac{1}{2}, T=1, t_0=0, t=1, u=0.1, n_{min, B}=1 $ and $ K=1 $. As shown in Fig. 1, the correlation function has the largest values at harder mode k in the case of $ v<0 $ than the one of $ v>0 $.

Figure 1.  (color online) The correlation function of net-baryon density in Fourier space for the system in the spinodal region ($ v=-1.5 $) and outside the spinodal region ($ v=0 $ and $ v=1 $).

To illustrate the effects of the unstable mode on the multiplicity fluctuations acceptance $ \Delta x $, the correlation function in Fourier space Eq.(8) needs to be transformed into coordinate space. After performing the spatial integration with interval $ (-\Delta x/2, \Delta x/2) $, the correlation function reads:

where the function $ f(k, \Delta x) $ is defined as

The function $ f(k, \Delta x) $ arises from the inverse Fourier transformation and spatial integration with finite range. In the slow mode limit (small k), the function $ f(k, \Delta x) $ is just a quadratic function of $ \Delta x $: $ f(k, \Delta x)\approx 5\Delta x^2 $. While $ f(k, \Delta x) $ is a trigonometric function in the limit of large k: $ f(k, \Delta x)\approx 4\frac{\Delta x}{k}\sin(k\Delta x) $. As shown in Fig.2, the function $ f(k, \Delta x) $ for small k behaves as a quadratic function with increasing $ \Delta x $. On the contrary, the function has oscillating behavior with larger k. This is consistent with the conventional picture that the harder mode varies with smaller wavelengths. In this case, the criteria of the existence of oscillating behavior are determined by the negative second order derivative of $ f(k, \Delta x) $ with respective to $ \Delta x $: $ f''(k, \Delta x)=10\cos(k\Delta x)-4k\Delta x\sin(k\Delta x)<0 $. As shown in Fig.3, the curvature $ f''(k, \Delta x) $ is positive at small spatial interval $ \Delta x $ in the case of slow mode $ k=0.1 $. Namely, no oscillating behavior of $ f(k, \Delta x) $. On the other hand, in the case of hard mode ($ k=1 $ in Fig.3), the second term $ -4k\sin(k\Delta x) $ plays a more important role (oscillating behavior appears), and $ f''(k, \Delta x) $ can be negative with small spatial interval $ \Delta x $.

Figure 2.  (color online) The function related to the inverse Fourier transform as a function of spatial integration interval $ \Delta x $ (see Eq.(10) for the definition of $ f(k, \Delta x) $). The blue dotted dashed curve corresponds to $ 0.2f(k, \Delta x) $ with slow mode $ k=0.1 $, while the red solid one represents the harder mode $ k=1 $.

Figure 3.  (color online) Similar with Fig.2, but for the second derivative of $ f(k, \Delta x) $ with respect to $ \Delta x $: $ f''(k, \Delta x) $.

As it has been stressed, the unstable mode in the spinodal decomposition region ($ v<0 $) derives the correlation function in Fourier space $ C(t, k) $ dominates at large k, while $ C(t, k) $ rapidly decays at large k outside the spinodal region. Due to the property of $ f(k, \Delta) $, it is natural to expect that the correlation function $ C_2(\Delta x) $ exhibits oscillating behavior with increasing system size $ \Delta x $ in the spinodal region, which does not exist outside the spinodal region. As an illustration, the correlation function with increasing system size $ \Delta x $ is shown in Fig.4. Note that the values of v in Fig.4 are chosen just for illustrative purpose. The oscillating behavior continuously smears to the ones of $ v=0 $ with decreasing absolute value of negative v. The oscillating behavior of $ C_2(\Delta x) $ for $ v<0 $ indicates the existence of an unstable mode induced by spinodal decomposition, which is hard to achieve outside the spinodal region ($ v>0 $) at small $ \Delta x $.

Figure 4.  (color online) The correlation function of net-baryon density with increasing spatial integration interval $ \Delta x $ for the system in the spinodal region ($ v=-1.5 $) and outside the spinodal region ($ v=0 $ and $ v=1 $).

In this work, $ \Delta x $ refers to the interval of spatial integration. Such analysis can also be performed in the relativistic frame ($ \tau, x, y, \eta $), where the integration is performed in the space-time rapidity space η. For the system with large collision energy, the space-time rapidity η equals the momentum rapidity y because of the boost-invariance. Detector with wider momentum rapidity window $ \Delta y $ approximately corresponds to the integration of the system with a larger system size. This so-called acceptance dependence with non-monotonic behavior of multiplicity fluctuations has been studied near the critical point on the crossover side [29−31]. As it is addressed in this work, the acceptance dependence with oscillating behavior of second-order multiplicity fluctuations is observed near the first-order phase transition side.

Based on the conserved dynamics near the phase transition, model B, this work focuses on the conserved net-baryon multiplicity fluctuations near the spinodal decomposition region. It is found that the convex anomaly induces unstable mode in the Fourier space, and the correlation function dominates at harder mode, comparing the one outside the spinodal region. As a result, this property picks up the higher mode with oscillating behavior at finite spatial interval. Such oscillating behavior of net-baryon multiplicity fluctuations with increasing detector acceptance could be the possible signature of the convex anomaly of spinodal decomposition.

Note that this work limits its goal to studying the possible consequence of unstable mode in spinodal decomposition with respect to the acceptance. More comprehensive and phenomenological investigation is required in heavy-ion experiments, in a more realistic context. It is essential to evaluate which mode of the correlation function dominates in the realistic setup and estimate the possibility of the survival of the unstable mode. For instance, the QGP fireball created in relativistic heavy-ion collisions is a rapidly expanding system, which derives the system passing through the spinodal decomposition. Such finite time effects have not been considered in this work, which is important for future study to address the possible existence of such a signature of spinodal decomposition.