# Spinodal instabilities of baryon-rich quark matter in heavy ion collisions

###### Abstract

Using the test-particle method to solve the transport equation derived from the Nambu-Jona-Lasino (NJL) model, we study how phase separation occurs in an expanding quark matter like that in a heavy ion collision. To test our method, we first investigate the growth rates of unstable modes of quark matter in a static cubic box and find them to agree with the analytical results that were previously obtained using the linear response theory. In this case, we also study the higher-order scaled density moments in the quark matter, which have values of one for a uniform density distribution, and they are found to increase with time and saturate at values significantly larger than one after the phase separation. The skewness of the quark number event-by-event distribution in a small sub-volume of the system is also found to increase, but this feature disappears if the sub-volume is large. For the expanding quark matter, two cases are considered with one using a blast-wave model for the initial conditions and the other using initial conditions from a mulple-phase transport (AMPT) model. In both cases, we find the expansion of the quark matter is slowed down by the presence of a first-order phase transition. Also, density clumps appear in the system and the momentum distribution of partons becomes anisotropic , which can be characterized by large scaled density moments and non-vanishing anisotropic elliptic and quadrupolar flows, respectively. The large density fluctuations further lead to an enhancement in the dilepton yield. In the case with the AMPT initial conditions, the presence of a first-order phase transition also results in a narrower rapidity distribution of partons after their freeze out. These effects of density fluctuations can be regarded as possible signals for a first-order phase transition that occurs in the baryon-rich quark matter formed in relativistic heavy ion collisions.

## I Introduction

Studying the properties of baryon-rich quark-gluon plasma (QGP) is the main focus of the beam energy scan (BES) experiments Nayak (2009); Aggarwal et al. (2010); McDonald (2015) at the Relativistic Heavy Ion Collider (RHIC) as well as at the future Facility for Antiproton and Ion Research (FAIR). These experiments are expected to shed light on whether the phase transition from the baryon-rich QGP to the hadronic matter is a first-order one and the location of the critical end point in the QCD phase diagram if the phase transition is first-order. To help understand what could happen in a baryon-rich QGP, we have recently used the Polyakov-Nambu-Jona-Lasinia (PNJL) model to study its spinodal instability Li and Ko (2016). We have found via the linear response theory that the spinodal boundary in the temperature and density plane of the QCD phase diagram shrinks with increasing wave number of the unstable mode and is also reduced in the absence of the Polyakov loop. In the small wave number or long wavelength limit, the spinodal boundary coincides with that determined from the isothermal spinodal instability in the thermodynamic approach. We have further found that the quark vector interaction suppresses unstable modes of all wave numbers. For the wave number dependence of the growth rate of unstable modes, it initially increases with the wave number but decreases when the wave number is large. For the collisional effect from quark scattering, we have included it via the linearized Boltzmann equation and found it to decrease the growth rate of unstable modes of all wave numbers. In the present study, we continue the above work by investigating how unstable modes would grow if one goes beyond the linear response or small amplitude limit. This is carried out by using the transport equation derived from the NJL model to study the time evolution of density fluctuations in a confined as well as in an expanding quark matter. Specifically, we study the time evolution of higher-order density moments in the quark matter, the distribution of quark number in a sub-volume of the quark matter, the quark momentum anisotropy, and dilepton production rate from quark-antiiquark annihilation. As shown below, these observables could serve as signatures for a first-order phase transition of the baryon-rich quark matter produced in heavy ion collisions.

The paper is organized as follows: In the next section, we give a brief review on the NJL model and the transport equations based on its Lagrangian. The transport equations are solved by the test-particle method in Section III to study both the short and long time behavior of the spinodal instability of a quark matter in a periodic box. The same method is applied in Section IV to an expanding quark matter to study how density fluctuations are affected by the expansion of the system as in heavy ion collisions. Finally, a summary is given in Section V. In the Appendix, we describe in detail the effect due to the finite grids used in the numerical calculations on the growth rate of unstable modes.

## Ii The NJL Lagragian And The Transport Model

The NJL Lagrangian containing only the scalar interaction for three quark flavors has the form Bratovic et al. (2013):

(1) |

where , and are the Gell-Mann matrices for , with being the identity matrix multiplied by . The Lagrangian preserves symmetry but breaks the axial symmetry, which is broken in QCD by the axial anomaly, by the Kobayashi-Masakawa-t’Hooft (KMT) interaction given by the last term in Eq. (II) ’t Hooft (1976). The in this term denotes the determinant in the flavor space Buballa (2005), that is

(2) |

where denotes either a Dirac gamma or the identity matrix. The determinantal term is responsible for obtaining the correct splitting in the masses of and mesons.

Because the NJL model is not renormalizable, a regularization scheme is required to remove infinities in the momentum integrations. In this study, we assume that all interactions are among quarks of 3-momenta with magnitudes below the cutoff momentum . Taking , the values of the scalar coupling and the KMT interaction can be determined from fitting the pion mass, the kaon mass, and the pion decay constant, and their values are , and if the current quark masses are taken to be and MeV Holstein (1990).

A flavor-singlet vector interaction can be added to the NJL Lagrangian as follows:

(3) |

where the coupling strength is assumed to be independent of the temperature and the net quark chemical potential . The value of affects the order of quark matter phase transition. If is large, the first-order phase transition induced by the attractive scalar interaction could disappear Li and Ko (2016). In the present study, we treat it as a parameter to change the equation of state of quark matter.

For describing the quark matter produced in a heavy ion collision, we use the Boltzmann (or transport) equations that can be derived from the NJL Lagrangian in terms of the non-equilibrium Green’s functions for quarks and antiqaurks Klevansky et al. (1997), and they are:

(4) |

where is the phase-space distribution function of quarks or antiquarks of flavor . In the above, is the kinetic momentum with the subscript referring to quarks and referring to antiquarks, is the vector potential, and is the effective quark mass with being the scalar potential with .

The right hand side of Eq.(II),

(5) | |||||

is the collisional term that describe the scatterings among quarks and antiquarks, with the subscripts , , , and now denoting not only the flavor but also the spin and color, and baryon charge (quark or anti-quark) of a parton. The above equation can be solved using the test particle method Wong (1982) by expressing the distribution function in terms of the density of test particles, whose equations of motions are determined by the left hand side of Eq.(II), and they are

(6) | |||||

(7) |

The second equation in the above can also be written as

(8) |

where is the strong magnetic field and is the strong electric field.

Besides the mean fields, test particles are also affected by collisions, which can be treated geometrically by generalizing the method of Ref. Bertsch and Das Gupta (1988) to use the particle scattering cross section in the quark matter frame to check whether the impact parameter between two colliding particles is smaller than and if the two colliding particles pass through each other at the next time step during the evolution of the system. For two particles of masses and , momenta and , and energies and , this cross section is related to the cross section in their center-of-mass frame with being the square of their invariant mass, which is the one used in Ref. Bertsch and Das Gupta (1988), by

(9) |

In the above, and are the velocities of the two particles. The 3-momenta of the two particles after the scattering are taken to be isotropic in their center-of-mass frame. Because of the high quark baryon chemical potential considered in the present study, the Pauli blocking effect on scatterings is also included by checking the available phase space for the final states Bertsch and Das Gupta (1988). We have checked that the above treatment of parton scattering reproduces the expected scattering rate evaluated via direct numerical integrations.

## Iii Quark Matter in a Box

This section serves as a bridge between the studies of the spinodal instabilities in the small and large amplitude limits. Although the case of small amplitude has already been discussed in Ref. Li and Ko (2016), we can develop an intuitive picture for how an initial sinusoidal fluctuation in a baryon-rich quark matter grows during the early stage of its time evolution from solving the Boltzmann equation in the test particle method as discussed in the previous section. For the large amplitude case, which also includes the growth of instabilities during the late stage, solving the Boltzamnn equation allows us to follow the whole phase separation process to see how dense clusters develop inside a box of initially uniform quark matter and finally lead to the formation of large scale structures. It also provides the possibility to find the appropriate observables to characterize these structures.

### iii.1 Small amplitude density fluctuations

We consider a quark matter that is confined in a cubic box with periodic boundary conditions. The system is prepared by distributing many test particles inside the box according to the density of the system with their momenta given by the Fermi-Dirac distribution at certain temperature. We then study the growth of density fluctuations from an initial distribution with density and temperature corresponding to that inside the spinodal region. Results obtained from solving the Boltzmann equation by following the classical motions of these test particles are compared with those obtained from the linear response theory in Ref. Li and Ko (2016). Specifically, we introduce an initial density fluctuation that has a sinusoidal oscillation in the z direction, , where is the average initial density and is the length of the box with corresponding to wave numbers , respectively. As an example, Fig. 1 shows how the amplitude of the sinusoidal wave grows with time in the case of , the average density , and an initial temperature MeV. Since the amplitude of density fluctuation at early times is expected to grow exponentially, it can be approximated by a hyperbolic cosine function of time, i.e.,

(10) |

where is the growth rate and can be extracted directly from the numerical results, and they are shown in Fig. 2 by solid circles. They are seen to agree very well with those obtained from an analytical calculation based on the linearized Boltzmann equation Li and Ko (2016) after including the finite grid size effect as described in the Appendix, shown by the solid and dashed lines for the cases with and without the collision term in the Boltzmann equation, respectively.

### iii.2 Large amplitude density fluctuations

To study how density fluctuations emerge and grow, we compare results from two calculations based on the same initial conditions but with and without the spinodal instability in the equation of state. This is achieved by introducing a vector interaction in the NJL model, which is known to move the state of a quark matter from inside the spinodal region to the outside if its strength is sufficiently large Li and Ko (2016). For example, for a quark matter of temperature 20 MeV and net quark density , the spinodal region disappears if the vector coupling has the same value as the scalar coupling , although the state of the quark matter is well inside the spinodal instability region for .

Figure 3 shows the time evolution of the density distribution in a box of size for the two cases of (left column) and (right column), with the darker color denoting the high density regions and the lighter color denoting the low density regions. Although the system is initially uniform in space, some dense spots are present due to statistical fluctuations as a result of finite number of test particles used in the calculation. In the case of without a first-order phase transition or spinodal instability, the density distribution in the box remains unchanged with time as shown in the right column. This changes dramatically, however, for the case of . Due to the spinodal instability, the initial dense spots act like ”seeds”, which create several small low pressure regions that attract nearby partons and lead to the formation of many clusters at fm. These clusters further grow in size by connecting with each other and form stable large structures at fm, when the system clearly separates into two phases of matter with one of high density and the other of low density.

A clearer picture can be obtained by taking a cross sectional view on the plane as shown by the density distribution contours in Fig. 4. The two phases are now distinguishable with the dilute phase having a density of about and the dense phase having a density of about . According to the phase diagram in Fig. 5, the initial location of the system is indicated by the circle inside the spinodal region. During the phase separation, the location of most part of the system moves towards the left boundary of the spinodal instability region that has a density of about , while that of the small part of the system moves towards the right boundary of the spinodal instability region that has a density of about , consistent with the picture shown by the density evolution.

As the large scale structure forms, we expect the density-density correlation to get stronger and the correlation length to become larger. This is indeed the case as shown in Fig. 6, where it is seen that both the amplitude of the correlation function and the correlation length increases with time.

The density fluctuations can also be quantified by the scaled density moments Steinheimer and Randrup (2013), where

(11) |

This quantity is scale invariant since its value remains unchanged under a scale transfomation , where can be any positive number. The scaled density moments are all equal to one for a uniform density distribution but become greater than one as the density fluctuations grow. In Fig. 7, we show by dotted, dashed, and solid lines the scaled density moments for , 4 and 6, respectively. Our results show that the scaled moments increase during the phase separation and reach their saturated values at about fm, when the phase separation almost ends. Also, moments with larger increase faster and saturate at larger values. The final saturation values can be estimated as follows. For a system of an initial density that separates into two phases of density and with volumes and , respectively, the scaled density moments are then

(12) |

Using the condition of particle number conservation

(13) |

the scaled density moments after the phase separation is thus

(14) |

For our case of , , and , we have ,, and , which are close to the final saturation values shown in Fig. 7.

Other quantities of interest are the skewness and kurtosis of the particle multiplicity distribution, which were proposed as possible signals for the critical phenomena Stephanov (2009) and have been studied in the beam energy scan experiments at RHICNayak (2009); Aggarwal et al. (2010). They are defined as follows:

(15) |

Both quantities characterize how far an event-by-event multiplicity distribution deviates from a normal distribution. A positive skewness means a long tail on the right side of the distribution, i.e., most events have the net quark number below the mean value, while some events have an extreme high net quark number. A positive kurtosis implies a sharper peak than the peak in a normal distribution, while a negative kurtosis corresponds to a flatter one. Theoretical calculations based on the grand canonical picture predict that both quantities diverge with the correlation length when a system approaches its critical point Stephanov (2009), with the kurtosis diverging faster than the skewness. Therefore, they have thus been suggested as the signals for the existence of a critical end point in the QCD phase diagram.

To be consistent with the grand canonical picture, we consider quarks in a sub-volume of the box in our study, such as its central cell, and treat the remaining part as the reservoir. When the system is initially inside the spinodal instability region, quarks in the reservoir can sometimes move into the sub-volume, but in most of the times quarks would leave from the sub-volume to the reservoir. The number of quarks inside this sub-volume thus varies drastically from event to event, leading to large values for the skewness and kurtosis in its event-by-event distribution. In Figs. 8, we show the event-by-event distribution of the number of quarks in the central cell from 1000 events at , , and fm by the solid, dashed and dotted lines, respectively, for the two cases of sub-volume of size (upper window) and (lower window). The upper window of Fig. 8 clearly shows that the distribution for the small sub-volume becomes asymmetric as time increases, starting with an initial skewness of 0.11 and increasing to 0.60 at 20 fm and 0.75 at 40 fm/. This feature is absent in the lower window of Fig. 8 for the larger sub-volume, where the distribution remains essentially symmetric with increasing time, with the skewness changing slowly from -0.001 (t=0) to 0.086 (t=20 fm) and 0.132 (t=40 fm), and there is no apparent increase or decrease in the kurtosis.

## Iv Expanding Quark Matter

### iv.1 Blast wave initial conditions

To study how large density fluctuations due to the spinodal instability as a result of a first-order phase transition obtained from the box calculation in the previous section are affected by the expansion of the system as in a heavy ion collision, we carry out a dynamical calculation using the transport model that includes parton scatterings besides the mean-field potentials described in Section II. For the initial parton distributions, their positions are taken to follow that of a spherical Wood-Saxon form:

(16) |

with a radius and a surface thickness parameter , similar to that expected from a central Au+Au collisions. The momenta of these patons are again taken to be that of a Fermi-Dirac distribution at certain temperature. Calculations are then carried out with two different equations of state with and without a first-order phase transition, which can be realized by adjusting the coupling strength for the vector interaction.

To see how the expanding system goes into the spinodal region in the QCD phase diagram, we first study the time evolution of the temperature and net quark density in the central volume of 42.875 fm, which has an initial density and temperature MeV, and trace its phase trajectory as shown in Fig. 9 for the two cases with (solid line) and without (dashed line) a phase transition. Although the quark matter described by the transport model may not always be in perfect thermal equilibrium, we approximate its temperature by that of an equilibrated one that has the same energy density and net quark density in the NJL model. As expected, the quark matter with a first-order phase transition (solid curve) enters the spinodal instability region, which is shown by the gray color, at about and leaves the region at about after spending about inside this region. How the central density decreases with time is shown by the solid line in Fig. 10, which is seen to decrease slower than in the case without a first-order phase transition shown by the dashed line obtained with

The density fluctuations can be seen from the density distribution on a plane such as the one at shown in Fig. 11. The left window shows the density distribution at fm/ for the case with a first-order phase transition, while the right window shows that at fm/ for the case without a first-order phase transition, when the density of the central cell is about fm in both cases. Although density clumps appear in both cases, those in the one with a first-order phase transition are significantly larger. As in the case of quark matter in a box, we can quantify the density fluctuations by the scaled density moments Randrup (2010). They are shown in Fig. 12 by the black and red lines for the cases with and without a first-order phase transition, respectively. The dotted, dashed, and solid lines are for , , and , respectively. In both cases, the scaled density moments first increase and then decrease with time. In the case without a first-order phase transition, this is caused by the fast increase of the surface of the quark matter and the quick deviation from its initial smooth Wood-Saxon density distribution. To the contrary, the scaled density moments in the case with a first-order phase transition becomes much larger with time and only decreases slightly afterwards, reflecting the effect due to density clumps that distribute randomly inside the expanding quark matter. Therefore, the saturated scaled density moments, which are larger for larger , can be regarded as signals for a first-order phase transition in a baryon-rich quark matter Steinheimer and Randrup (2013).

Since density fluctuations can lead to spatial anisotropy even in central heavy ion collisions, it has been suggested that they may affect the anisotropic flows in the transverse plane Herold et al. (2014); Steinheimer et al. (2014). The latter are defined by the coefficients in the expansion of the transverse momentum distribution as a Fourier series in the azimuthal angle ,

(17) |

where is the event plane angle Alver et al. (2010). To calculate the anisotropic flow coefficients, we use the two particle cumulant method Wang et al. (1991); Borghini et al. (2001), namely, by averaging over all particle pairs in an event. We have calculated and for 100 events of an expanding quark matter with the same blast wave initial conditions, and their final event distributions are shown, respectively, in the upper and lower windows of Fig. 13 with the solid and dashed lines for the cases with and without first order phase transition, respectively. Both distributions peak at a larger value for the case with a first-order phase transition, particularly for , thus providing a plausible signal for the first-oder phase transition. However, the values of the fluctuation induced and are much smaller than those in non-central heavy ion collisions.

We have also studied the effect of density fluctuations on dilepton production from a quark matter. Since the dilepton production rate is proportional to the square of parton density, more dileptons are produced when the density fluctuation is large. Also, a longer partonic phase as a result of a first-order phase transition would increase the depletion yield as well. As usually done in studying dilepton production in heavy ion collisions Xiong et al. (1990), we use the perturbative approach to calculate the dilepton yield from the quark-antiquark scattering by neglecting its effect on the dynamics of the expanding quark matter. Using the dilepton production cross section,

(18) | |||||

where is the square of the dilepton invariant mass, we have calculated the dilepton invariant mass spectrum from the expanding quark matter, and they are shown in Fig. 14 by the solid and dashed lines for the cases with and without first-order phase transition, respectively. As expected, more dileptions are produced from the quark matter with a first-order phase transition. We note the dilepton invariant mass spectrum peaks at GeV with the peak value being about GeV, which is comparable with the result obtained from a hadronic transport model Galatyuk et al. (2015). This enhancement in dilepton production may thus be detectable in experiments. We also note that most dileptons are produced from quark-antiquark annihilation as very few pions are present in the system due to the low phase transition temperature in the SU(3) NJL model.

### iv.2 AMPT initial conditions

In this subsection, we use a more realistic initial parton distribution for heavy ion collisions. Specifically, the initial partons are obtained from a multiphase transport (AMPT) model with string melting Lin et al. (2005) that uses the heavy ion jet interaction generator (HIJING) Wang (1991); Wang and Gyulassy (1991); Gyulassy and Wang (1994) as the input. This model includes not only the mini-jet partons from initial hard collisions but also hadrons produced from excited strings, which are projectile and target nucleons that have suffered interactions, by converting them to partons according to the flavor and spin structures of their valence quarks. In particular, a meson is converted to a quark and an anti-quark, while a baryon is first converted to a quark and a diquark, and the diquark is then decomposed into two quarks. The quark masses are taken to be , , and as in the PYTHIA program SjÃ¶strand (1994). The above two-body decomposition is isotropic in the rest frame of the parent hadron or diquark. These partons are produced after a formation time of , with and denoting, respectively, the energy and transverse mass of the parent hadron. We obtain these partons as the initial conditions for our study of an expanding quark matter by running the AMPT program with vanishing parton scattering cross sections in Zhang’s parton cascade (ZPC)Zhang (1998) and with the hadronic afterburner based on a relativistic transport (ART) Li and Ko (1995); LI et al. (2001) turned off. Using the partons from Au+Au collisions at zero impact parameter and a center-of-mass energy GeV as the initial distribution, we have found that some parts of the system go through the spinodal region when the SU(3) NJL model with is used in the Boltzmann equation and in constructing the phase diagram.

As shown by the solid line in Fig. 15, the trajectory of the central part of the system goes into the spinodal instability region at about after expansion, and moves out of this region at about . Although is too short for the spinodal instability to develop in the central part of the quark matter, its other parts may stay longer in the spinodal instability region due to both the spatial distribution of initial partons and the correlations between their rapidities and longitudinal () coordinates.

Figure 16 shows the rapidity and longitudinal coordinate correlation of initial partons from a typical AMPT event for central Au+Au collisions at . This correlation can be quantified as follows:

(19) |

This positive correlation indicates that partons initially on the right side of the quark matter are more likely to have momenta pointing to the right or forward direction, while partons initially in the left of the quark matter are more likely to have momenta pointing to the left or backward dirction. This correlation helps the initially disc-shaped quark matter to expand, leading to a fast decrease of the density in the center of the quark matter as shown in the upper row of Fig. 17. Here, the quark matter is initially largely confined in a thin disk of thickness less than 0.5 fm. When it is allowed to free streaming without any interactions, there appear two high density clumps that fly apart in the opposite directions. This feature becomes less prominent after the inclusion of quark scattering and mean-field potentials but without a phase transition in the quark matter, i.e., taking , as shown in the middle row of Fig. 17. With a first-order phase transition in the quark matter by setting , the lower row of Fig. 17 shows that the initial central disk evolves into three disks of dense matter with one in the middle due to the strong attractions that keep some partons from moving away, besides the two forward and backward moving disks. As the quark matter expands, these disks transform into rings and finally turn into disjointed clumps. Furthermore, the density distribution of the quark matter in the reaction plane () shown in Fig. 18 indicates that the quark matter with a first-order phase transition expands twice as slow as that without a first-order phase transition.

Because of the non-trivial spatial distribution even in the case of free-streaming quark matter, the scaled density moments are no longer useful quantities to characterize the density fluctuations of an expanding quark matter due to its spinodal instability or a first-order phase transition. On the other hand, the different density variations along the beam () axis shown in Fig. 18 are expected to affect the parton rapidity distribution. This is because partons in the middle disc, which is present only in the case with a first-order phase transition, have a small rapidity and due to the attractive quark interactions, they attract partons from the other two discs and slow down their expansion in the longitudinal direction, thus restricting their rapidities to a narrow region around the midrapidity. As shown by the solid line in Fig. 19, the parton rapidity distribution in the case with a first-order phase transition is indeed much narrower than that in the case without a first-order phase transition, shown by the dashed line. This effect can be regarded as a possible signal of a first-order phase transition and is worth studying in experiments.

We have also studied the dilepton invariant mass spectrum from an expanding quark matter with initial conditions from the AMPT model. This is shown in Fig. 20 by the solid and dashed lines for the cases with and without a first-order phase transition, respectively. As in the previous section using the blast-wave initial conditions, the presence of a first-order phase transition enhances the dilepton yield as a result of density fluctuations and a longer partonic phase. However, the dilepton yield is lower than that obtained from the calculation with the blast wave initial condition by two orders of magnitude because there are very few antiquarks in the partonic matter produced in heavy ion collisions at such a low energy and also because we have not included the bremsstrahlung contribution to dilepton production from the quark-quark scattering.

## V conclusions

The spinodal instability is a thermodynamic feature of a first-order phase transition in a many-body system. It occurs when its pressure in some parts decreases with increasing density. This can amplify the density fluctuations and lead to a phase separation in the system. We have studied this phenomenon by solving the Boltzmann equations using the test particle method. The calculations are based on the NJL model, which has been shown to give good a description of the vacuum properties of the hadrons and also predicts the existence of a first-order phase transition in baryon-rich quark matter. We have obtained some intuitive pictures on the phase separation in a quark matter that is either in a static box or undergoes expansion. For the case of a static box, we have found that the growth rates extracted from the early growth of a sinusoidal density fluctuation agree with the analytical results obtained from the linearized Boltzmann equation. We have also calculated the higher-order density moments of the quark matter and found them to increase and saturate at large values after phase separation, making them possible signals for the first-order phase transition. The skewness of the quark number event-by-event distribution in a small sub-volume of the quark matter is also found to increase, but this feature disappears if the sub-volume is large. As for the expanding quark matter, two cases have been studied. One is based on the blast-wave initial conditions, while the other using the AMPT initial conditions, which are disc-like as a result of the strong correlations between the parton rapidity and longitudinal coordinate. In both cases, we have found that the expansion of the quark matter is slowed down by the presence of a first-order phase transition. Density clumps are found to appear and lead to an anisotropy in the momentum space, which can be characterized by the scaled density moments and the anisotropic flows and , respectively. An enhancement in the dilepton yield is also observed. The expansion of the quark matter with the AMPT initial conditions is more complex. Normally, the initial disc-like quark matter splits into two discs, moving along the beam axis in opposite directions. If the expanding quark matter undergoes a first order-phase transition, a third disc appears in the middle and pulls the other two discs towards it, resulting in a narrower rapidity distribution.

In the future, we plan to develop a more consistent transport model, in which all cross sections are calculated self-consistently from the NJL model, so that the temperature and density dependence of the collisional effect can be taken into account. The dilepton production through the process will also be included, since it could be the main contribution to the dilepton yield from a quark matter of high baryon chemical potential. We also plan to extend the transport model using the PNJL modelFukushima (2008), which is more realistic and agrees better with the lattice results for a quark matter with low baryon chemical potential. We hope that our study will help to understand the phase transition in the baryon-rich matter by comparing theoretical predictions with available and future experimental data.

## Acknowledgements

This work was supported by the US Department of Energy under Contract No. DE-SC0015266 and the Welch Foundation under Grant No. A-1358.

*

## Appendix A Finite grid size effects

Counting partons in a grid of finite size in evaluating the mean fields effectively allows the partons in the grid interact with each other, thus modifying the contact interactions in the NJL model to finite-range ones. To study this effect, we need to calculate the probability for two partons in the same grid to have a separation . Given a parton located at in a 1-dimensional grid , the probability to find another parton located at in the same grid is

(20) | |||||

where

(21) |

The above expression can be straightforwardly generalized to the 3-dimensional case to give

(22) |

where are the grid lengths. The interaction between two partons at and is then replaced by

(23) |

Transforming Eq. (A) from -space to -space gives

(24) |

Note that in the limit that for all the , and , which means the modification does not affect the long wavelength modes.

Replacing and in Eq. (35) in Ref. Li and Ko (2016) with and , respectively, and solving the resulting equation, we obtain the modified dispersion relation, and they are shown in Fig. 2 for a grid size by the solid and dashed lines for the cases with and without the collision term, respectively. As expected, the growth rate is not much affected in the small region but is significantly suppressed in the large region. The finite grid size effect is thus similar to the quantum effect shown in Ref. Li and Ko (2016). Using a finite grid size essentially allows partons to interact at finite separation, resulting in an effective finite-range interaction.

## References

- Nayak (2009) T. K. Nayak (STAR Collaboration), Nucl. Phys. A 830, 555C (2009).
- Aggarwal et al. (2010) M. M. Aggarwal et al. (STAR), Phys. Rev. Lett. 105, 022302 (2010).
- McDonald (2015) D. McDonald (STAR), EPJ Web Conf. 95, 01009 (2015).
- Li and Ko (2016) F. Li and C. M. Ko, Phys. Rev. C 93, 035205 (2016).
- Bratovic et al. (2013) N. M. Bratovic, T. Hatsuda, and W. Weise, Phys. Lett. B 719, 131 (2013).
- ’t Hooft (1976) G. ’t Hooft, Phys. Rev. D 14, 3432 (1976).
- Buballa (2005) M. Buballa, Phys. Rept. 407, 205 (2005).
- Holstein (1990) B. R. Holstein, Phys. Lett. B 244, 83 (1990).
- Klevansky et al. (1997) S. P. Klevansky, A. Ogura, and J. Hufner, Ann. Phys. 261, 37 (1997).
- Wong (1982) C.-Y. Wong, Phys. Rev. C 25, 1460 (1982).
- Bertsch and Das Gupta (1988) G. F. Bertsch and S. Das Gupta, Phys. Rept. 160, 189 (1988).
- Steinheimer and Randrup (2013) J. Steinheimer and J. Randrup, Phys. Rev. C 87, 054903 (2013).
- Stephanov (2009) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
- Randrup (2010) J. Randrup, Phys. Rev. C 82, 034902 (2010).
- Herold et al. (2014) C. Herold, M. Nahrgang, I. Mishustin, and M. Bleicher, J. Phys. Conf. Ser. 509, 012065 (2014).
- Steinheimer et al. (2014) J. Steinheimer, J. Randrup, and V. Koch, Phys. Rev. C 89, 034901 (2014).
- Alver et al. (2010) B. H. Alver, C. Gombeaud, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 82, 034913 (2010).
- Wang et al. (1991) S. Wang, Y. Z. Jiang, Y. M. Liu, D. Keane, D. Beavis, S. Y. Chu, S. Y. Fung, M. Vient, C. Hartnack, and H. Stöcker, Phys. Rev. C 44, 1091 (1991).
- Borghini et al. (2001) N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, Phys. Rev. C 64, 054901 (2001).
- Xiong et al. (1990) L. Xiong, Z. G. Wu, C. M. Ko, and J. Q. Wu, Nucl. Phys. A512, 772 (1990).
- Galatyuk et al. (2015) T. Galatyuk, P. M. Hohler, R. Rapp, F. Seck, and J. Stroth (2015).
- Lin et al. (2005) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005).
- Wang (1991) X.-N. Wang, Phys. Rev. D 43, 104 (1991).
- Wang and Gyulassy (1991) X.-N. Wang and M. Gyulassy, Phys. Rev. D 44, 3501 (1991).
- Gyulassy and Wang (1994) M. Gyulassy and X.-N. Wang, Comput. Phys. Commun. 83, 307 (1994), ISSN 0010-4655.
- SjÃ¶strand (1994) T. SjÃ¶strand, Comput. Phys. Commun. 82, 74 (1994), ISSN 0010-4655.
- Zhang (1998) B. Zhang, Comput. Phys. Commun. 109, 193 (1998), ISSN 0010-4655.
- Li and Ko (1995) B.-A. Li and C. M. Ko, Phys. Rev. C 52, 2037 (1995).
- LI et al. (2001) B.-A. LI, A. T. Sustich, B. Zhang, and C. M. Ko, Int. J. Mod. Phys. E 10, 267 (2001).
- Fukushima (2008) K. Fukushima, Phys. Rev. D 77, 114028 (2008), [Erratum: Phys. Rev.D78,039902(2008)].