Spinodal decomposition in Bjorken flow

Mapping the QCD phase diagram has been a priority for the heavy-ion collision community. Lattice calculations with non-zero quark masses have shown that the QCD matter changes phase from a hadron gas to the quark gluon plasma (QGP) by a smooth crossover when the baryon chemical potential is near zero Aoki:2006we . At high baryon densities, the phase change is expected to proceed via a first-order phase transition Fukushima:2010bq . The first-order phase transition curve ends in a critical point before reaching the zero chemical potential axis.

The Beam Energy Scan (BES) program at the Relativistic Heavy-Ion Collider (RHIC) and the upcoming experiments at the Facility for Antiproton and Ion Research (FAIR) aim to explore different regions of the QCD phase diagram with the goal of quantitatively mapping its features. Heavy nuclei are collided at a range of collision energies which cover a wide region of the phase diagram. This program relies on identifying the observable signatures of the critical point and the phase transition curve and doing extensive model-to-data comparisons An:2021wof . We expect the phase transition curve to cover a bigger region of the phase diagram than the critical point and hence it is crucial that we understand its dynamics.

The phase transition in these collisions is expected to proceed via spinodal decomposition. Here we extend the equation of state to metastable and unstable regions and report on the equations describing spinodal decomposition in a relativistic QCD medium. We write and solve these equations in a fluid with Bjorken flow.

An important step towards describing the spinodal decomposition is extending the equation of state to the coexistence phase. This phase is often not described by the first principle approaches as it is composed of metastable and unstable regions. Given an equation of state for the stable region, we can interpolate the pressure in this phase as Kapusta:2024nii

$$P_{{\rm int}}(n)=P_{X}(T)+\sum_{i=1}^{4}c_{i}(n-n_{G})^{i}.$$

(1)

The subscript $X$ denotes the pressure along the line where the phases coexist. This line is often defined by Maxwell construction. The baryon density is denoted by $n$ while $n_{G}$ denotes the baryon density in the hadronic gas phase along the coexistence line. Similarly, the chemical potential can be parameterized as

$$\mu_{{\rm int}}(n)=\mu_{X}(T)+d_{0}\ln(n/n_{G})+\sum_{i=1}^{3}d_{i}(n-n_{G})^{i},$$

(2)

The coefficients $c_{i}$ and $d_{i}$ can be obtained by applying the thermodynamic relations at $n=n_{G}$ and $n=n_{L}$ where the equation of state is stable. Further constraints are obtained by applying the continuity conditions on pressure and chemical potential and their derivatives on the border between the coexistence and the QGP and hadronic gas phases. The thermodynamic relations also hold at the boundaries and provide additional constraints. The interpolated pressure and baryon chemical potential for a temperature of 100 MeV are shown in figure 1. The stable QCD equation of state is taken from reference Kapusta:2021oco .

The Helmholtz free energy can be written as

$$F\{n({\bf x},t)\}=\int d^{3}x\left[{\textstyle{\frac{1}{2}}}K(\bm{\nabla}n)^{2}+f(T,n)\right],$$

(3)

where first function on the right accounts for the free energy due to a gradient in baryon density across the phase boundary. The second term $f(T,n)$ is the bulk free energy. Cahn and Hilliard CahnHilliard1 ; CahnHilliard2 showed in their seminal papers that this leads to the stress tensor

$$T_{ij}=\tilde{P}\delta_{ij}+K(\partial_{i}n)(\partial_{j}n).$$

(4)

For a relativistic fluid, the covariant energy-momentum tensor can be given by Kapusta:2024nii

$$T^{\mu\nu}=\tilde{P}(u^{\mu}u^{\nu}-g^{\mu\nu})+\tilde{\epsilon}u^{\mu}u^{\nu}+K(D^{\mu}n)(D^{\nu}n),$$

(5)

where the local thermodynamic variables are pressure $\tilde{P}=P+KnD^{2}n+{\textstyle{\frac{1}{2}}}K(D^{\mu}n)(D_{\mu}n)$ and energy density $\tilde{\epsilon}=\epsilon-{\textstyle{\frac{1}{2}}}K(D^{\mu}n)(D_{\mu}n)$. Here $D^{\mu}=\partial^{\mu}-u^{\mu}u^{\alpha}\partial_{\alpha}$ is the gradient orthogonal to the fluid velocity $u^{\mu}$. We are using the mostly negative metric. In the Landau-Lifshitz frame baryon current has the form

$$J^{\mu}=nu^{\mu}+\sigma_{B}TD^{\mu}\left(\frac{\tilde{\mu}}{T}\right),$$

(6)

with baryon chemical potential $\tilde{\mu}=\mu+KD^{2}n$ and baryon conductivity $\sigma_{B}$.

We solve equations for boost invariant hydrodynamics with the simplification that system is homogeneous in transverse direction. We start the evolution at proper time $\tau_{x}=10$ fm when energy density is kept boost invariant. Baryon density is initialized as $n(\tau_{x})=n_{0}(\tau_{x})+n_{fl}\sin(k\xi)$ where $n_{0}$ is the average value of baryon density. Here $n_{fl}$ is the amplitude of sinusoidal fluctuations in baryon density in the rapidity direction $\xi$. The values of $n_{0}$ and $n_{fl}$ are chosen such that the lowest point of baryon density distribution just touches the coexistence phase region. Energy density gradient terms are neglected for simplicity. We verified that the violation of boost invariance in energy density is less than 3%.

The evolution of different phases and that of baryon density can be seen in figure 2. The baryon density reduces as the system expands and cools with time. It is interesting to note that unlike the steady-state systems undergoing spinodal decomposition, which have been extensively studied, this is a rapidly expanding system. So the system always dilutes and ends up entirely in the hadronic gas phase. If the phase change is via a first-order transition, then it encounters the coexistence phase.

The baryon density chosen here is smoothly varying. This is a simplification and ensures that the fourth-order derivatives encountered in the Cahn-Hilliard model are simplified. We chose a relatively small value of $K=5\times 10^{-5}$ MeV^-4 which ensured that this simplification holds for the duration of this evolution. In the future, this simplification can be relaxed for a more realistic baryon density distribution.

The energy density evolution is shown in figure 3. We show the results with $K=5\times 10^{-5}$ MeV^-4 and $K=0$. The later choice means that the phase surface has no energy contribution and is shown for reference. The energy increases immediately as the system enters the coexistence region. The energy ends up reducing in the negative rapidity when it goes in the coexistence phase. Whether the energy density increases or decreases depends on the sign of baryon density gradients. These surface energy contributions have the net effect of enhancing the fluctuations already present in the system.

The QCD phase diagram is postulated to have a first-order phase transition curve at large baryon densities. Experimental verification of such a phase transition relies on successful model-to-data comparisons with the models incorporating the phase transition dynamics. In heavy-ion collisions, such a phase transition is expected to proceed by spinodal decomposition. We wrote the equations of relativistic hydrodynamics with spinodal decomposition. The covariant energy-momentum tensor incorporating the surface effects is explicitly shown. We also provide an interpolation prescription to extend the QCD equation of state to metastable and unstable phases. We solved the equations for Bjorken flow.

The formalism can be included in the existing models of heavy-ion collisions to simulate the QCD phase transition. This will help us identify the observables containing the signatures of the phase transition and help map the QCD phase diagram.

This work was supported by the U.S. Department of Energy Grant Nos. DE-FG02-87ER40328 (JK, MS and TW) and DE-SC-0024347 (MS).