Chiral Phase Transition in an Expanding Quark System

We investigate the influence of chiral symmetry which varies along the space-time evolution of the system by considering the chiral phase transition in an non-equilibrium expanding quark-antiquark system. The chiral symmetry is described by the mean field order parameter, whose values is the solution of a self-consistent equation, and affects the space-time evolution of the system through the force term in the Vlasov equation. The Vlasov equation and the gap equation are solved concurrently and continuously for a longitudinal boost-invariant and transversely rotation-invariant system. This numerical framework enables us to carefully investigate how the phase transition and collision affect the evolution of the system. It is observed that the chiral phase transition gives rise to a kink in the flow velocity, which is caused by the force term in the Vlasov equation. The kink is enhanced by larger susceptibility and tends to be smoothed out by non-equilibrium effect. The spatial phase boundary appears as a"wall"for the quarks, as the quarks with low momentum are bounced back, while those with high momentum go through the wall but are slowed down.


I. INTRODUCTION
One of the major motivations in the study of finite-temperature Quantum ChromoDynamics (QCD) is to shed light on the phase structure of the strong interaction matter. Lattice-QCD data has confirmed the chiral and deconfinement phase transition is a crossover for small baryon chemical potential [1,2]. The sign problem at large baryon chemical potential [3,4] prevents lattice-QCD from giving precise predictions about the phase transition at finite density. Model calculations based on NJL model, quark meson model, and various beyond mean field frameworks have all predicted that the chiral phase transition at finite density is a first order phase transition [5][6][7][8]. Thermodynamic theorem then predicts a critical end point (CEP) between the crossover and the first order phase transition, which is a second order phase transition. However, due to various approximations adopted in the model calculations, there is not an agreement on the location of the CEP on the phase diagram.
The exploration of the QCD phase diagram is also one of the most important goals for relativistic heavy-ion experiments. Through a systematic measurement over a range of beam energies, the beam energy scan (BES) program makes it possible to search for the CEP in the QCD phase diagram [9][10][11]. Besides the ongoing BES program at RHIC, several other programs at other facilities such as FAIR and NICA have also contributed to the searching of QCD critical end point. The related experiments are mainly driven by measurements of net-proton or net-charge multiplicity fluctuations [12][13][14] which are expected to show characteristic non-monotonic behavior near the phase transition and especially near the CEP [15,16]. The fast dynamics in the fireball renders it difficult to bridge the gap between the experiment and the theories. Since from the theoretical aspect, the QCD phase structure is investigated in an equilibrium, long-lived, extremely large and homogeneous system. On the contrary, the fireball in the heavy ion collision is a highly dynamical system, characterized by very short lifetime, extremely small size and fast dynamical expansion. The finite-size effect and off-equilibrium effect prevent a divergent correlation length, and thus weaken the critical phenomenon which takes place in the equilibrium system. On one hand, the fast expansion and cooling during the evolution of the fireball tends to drive the system out of local equilibrium. On the other hand, the relaxation time diverges around the critical point [17], the critical slowing down renders it harder for the system to reach local equilibrium. A thorough understanding of phase transitions in the dynamical environment is thus of fundamental necessary to make profound predictions from the BES project.
Various models have been applied to study the chiral phase transition and related critical phenomenon in an out-of-equilibrium system. In Ref. [18,19], the authors investigate the chiral phase transition in a free-streaming quark-antiquark system by solving the Vlasov equation through the test particle method. The chiral fields are included considering their mean field values and their equations of motion. The Vlasov equation for the quarkantiquark system is also analytically solved in Ref. [20], assuming constant quark mass, and a shell-like structure at late evolution times in the center-of-mass (CM) frame is discovered. In Ref. [21] the nonequilibrium and collision effects on the deconfinement phase transition is investigated by solving the Vlasov equation assuming Bjorken symmetry. The elastic two-body collisions for the quarks and antiquarks is included by simulating a Vlasov-type of equation with MonteCarlo test-particle approach [22][23][24]. Among the aforementioned works, although the force term is also considered in the test particle method [18,19,[22][23][24], its effect has not been carefully examined, which we find plays an important role in the evolution of the system around the phase transition.
In order to study the chiral phase transition in the non-equilibrium state, we investigate an expanding quarkanitquark gas system by solving the coupled Vlasov equation as well as the gap equation. This paper is arranged as follows. In section III, we introduce and analyze the coupled Vlasov equation and gap equation which we are going to solve, and give the related thermal quantities that can be obtained as momentum integrals of distribution function. In section III, we consider a longitudinal boost invariant and transverse rotational symmetric system, and derive the transport equation under such condition. In section IV, we present our numerical process to solve the coupled equations and then analyze the numerical result. In section V, we summarize this work and give a brief outlook.

II. VLASOV EQUATION AND THERMAL QUANTITIES
The partons in an off-equilibrium systems with background field and collisions can generally be described by the Vlasov equation The above equation is applicable when the external fields and the interactions between the (quasi-)particles are sufficiently weak, so each particle can be considered to be moving along a classical trajectory, punctuated by rare collisions. As an example of an evolving global symmetry in an expanding parton system, we here consider the chiral symmetry in an off-equilibrium quark-antiquark system. At classical level, the velocity is v = p/E p , the energy of the quasi-particle is E p = p 2 + m(x) 2 , where m is the effective mass of the quark and antiquark, which is space-time dependent and is determined by the evolving chiral symmetry. The chiral mean field acts as a background field, and affects the motion of quarks through the gradient of the field energy F = ∇ x E p , which is a continuous force on the quarks. The mass of the quasi-particles m(x) is no longer a free parameter but is determined by the spacetime dependent chiral symmetry. The constituent quark mass serves as the order parameter of chiral symmetry. Its temperature dependence in equilibrium can be found in lattice-QCD simulation [1,2] and other model calculations.
We here consider the SU (2) Nambu-Jona-Lasinio Lagrangian (NJL) model [25][26][27], where ψ = (u, d) T is the two-component quark field in the flavor space, m 0 is the degenerate current mass of the quarks, τ is the Pauli matrix in the isospin space. The transport equations can be derived from a first principle theory or an effective model in the framework of Wigner function [28][29][30][31][32][33]. At the classical level, the quarks are treated as quasi-particles, and the chiral field is approximated by mean field, hence the Vlasov equation is coupled to the gap equation [33]. The disoriented chiral condensate (DCC) [34][35][36] is negligible here because it appears as a quantum effect, and the σ condensate appears as the mean field and couples to the Vlasov equation through the inhomogeneous quark mass. In order to investigate the collisions and non-equilibrium effect, we take the relaxation time approximation for the collision terms. The distribution function of the quark/antiquark number density f ± (t, x, p) satisfies the coupled Vlasov equation and gap equation, where the +(−) sign stands for quark (antiquark). In this work, we take the relaxation time approximation for the collision term, C[f ] = −(f ± − f ± eq )/τ θ , with f ± eq represents the corresponding local-equilibrium distribution function, while τ θ is the relaxation time. It is worth noticing that one can make the substitution f − (t, x, p) ≡ 1 − f − (t, x, −p) which follows the same equation of motion as f + . The transport equations can be further simplified by adopting the recombination f = f + + f − and g = f + − f − . In such way, the evolution of f and g can be separated. Since f − and f + satisfy the same transport equation, f (t, x, p) and g(t, x, p) also satisfy the same transport equation, while the gap equation depends only on f but not g. One thus solve the coupled transport equation of distribution function f (t, x, p) and g(t, x, p) as well as the gap equation for a finite density system, and solve transport equation of distribution function f (t, x, p) together with gap equation for a system with vanish baryon density, In the numerical procedures, we consider zero density system in this paper and solve both the transport equation of f (t, x, p) and the gap equation, and eventually get the time and space dependence of the distribution function and the quark mass. The second integral in the gap equation is the vacuum part, which has ultra-violet divergence and needs to be regularized. Here we take the hard cut-off regularization only for the vacuum part, The same cutoff Λ = 496 MeV is adopted for the longitudinal momentum and the transverse momentum in the integral. The NJL coupling constant is set to be G = 1.688/Λ 2 , so as to guarantee the quark mass m ∼ 300 MeV in the vacuum. The momentum integral of the finite temperature part is free from divergence, and is left unregularized. Under such choice of parameters, the temperature dependence of the quark mass in an equilibrium system can be directly calculated from the gap equation (7) by taking Fermi-Dirac distribution. The quark mass for chiral limit m 0 = 0 and real case m 0 = 3.7 MeV are presented in Fig. 1, the critical temperature is about 156 MeV at vanishing baryon chemical potential. The chiral phase transition at vanishing baryon chemical potential is a crossover in the real case, and is a second order phase transition in the chiral limit. For a crossover, the phase transition can be defined at the maximum susceptibility dm/dT ; for a second order phase transition, the susceptibility diverges at critical point. In an equilibrium system, the order parameter displays critical scaling m ∝ ((T c − T )/T c ) β in the vicinity of a second order phase transition which can be described by critical exponent β. The divergence of the susceptibility is expected to affect the transport phenomenon of the system through the force term ∇ r m 2 · ∇ p f , in which the force is provided by the ingredient of mass ∇ r m 2 . In an equilibrium system with zero baryon density, the mass is determined by temperature alone, and the force can be expressed as ∇ r m 2 = 2m(dm/dT )∇ r T . Since the susceptibility exhibits a peak around the phase transition, the temperature gradient ∇ r T in a realistic system is nonzero, the force term is expected be very large at the phase transition point in the time-space. However, when the system is close to the second order phase transition, the relaxation time may diverge [17], the critical slowing down takes place, makes it harder for the system to reach local equilibrium. When the system has not yet reached local equilibrium, the temperature is not well-defined, and neither for the expression 2m(dm/dT )∇ r T of the force term. If the phase transition takes place in an out-of equilibrium system, the aforementioned effects of the force term may have been overestimated. In the following, we study the chiral phase transition in both local equilibrium and out-of-equilibrium systems, controlled by the relaxation time. Then we analyze the influence of the phase transition and the force term on the evolution of the system.
The thermodynamic quantities can be constructed from the distribution function. For single component medium, positive particle f + for example, the current and energy-momentum stress tensor are defined by: which should satisfy that the conservation laws ∂ µ J µ = 0 and ∂ ν T µν = 0. Since the total energy density and entropy is the sum of that of positive particles and negative particles, while the net-quark number density is the difference of positive particles and negative particles, the above equations can be rewritten using the redefined distribution Using the Landau frame definition, the fluid velocity can be determined as the time-like eigenvector (u µ u µ > 0) of the stress tensor T µ ν u ν = u µ , with energy density being the corresponding eigenvalue. One could further obtain the particle number density and entropy density as n = u µ J µ and s = S µ u µ . The time-space evolution of (t, x), n(t, x) and u z (t, x) would give us quantitative idea about how the system evolves.
By taking the relaxation time approximation for collision kernel, we also need the corresponding local-equilibrium distribution function for any given time and space point, In the above distribution, the temperature and chemical potential are determined by matching the energy and number density, i.e. = eq and n = n eq , where with E p = m(T eq , µ eq ) 2 + p 2 . Similarly, we define the equilibrium limit of the entropy density as The difference between s eq and the actual entropy density s quantifies how close the system is to the equilibrium state. In a zero chemical potential limit considered in this paper, one would automatically find g = 0, hence J = 0, n = n eq = 0, µ eq = 0.

III. SYMMETRY AND SIMPLIFICATION
For numerical simplicity, we will focus on the longitudinal boost-invariant and transversal rotational-symmetric systems, which is a good approximation for ultra-central relativistic heavy-ion collisions. Under such symmetries, two constraints are applied to the system, and the distribution function f five degrees of freedom. We introduce a new set of coordinates (τ, η, ρ, φ, p ⊥ , ξ, θ), the original coordinates and the new ones can be transformed through the following relation, where ρ ∈ [0, +∞), p ⊥ ∈ [0, +∞) and ξ ∈ (−∞, +∞). Under the new set of coordinate, the longitudinal boostinvariance and transversal rotational-symmetry of the distribution function can be translated into its independence of φ and η, namely ∂ φ f = ∂ η f = 0. The phase space of the distribution function becomes (τ, ρ, p ⊥ , ξ, θ). The transport equation (6) is then reduced to where f eq is the equilibrium distribution. To clearly analyze the θ-dependence of the distribution function, and to simplify the calculation, we take the Fourier expansion of f (τ, ρ, p ⊥ , ξ, θ) with respect to θ, and also the equilibrium distribution function f eq (τ, ρ, p ⊥ , ξ, θ), where the Fourier coefficients are obtained by definition a n (τ, ρ, −π f (τ, ρ, p ⊥ , ξ, θ) sin(nθ)dθ, and similarly for A n and B n . In a realistic system, one can further expect its symmetry under reflection along eitherx-,ŷ-, orẑ-direction. Under such condition, one can show the θ-odd components of f and f eq vanish, b n ≡ B n ≡ 0, as well as a n (−ξ) = a n (ξ), A n (−ξ) = A n (ξ). Substituting the Fourier expansion (17) back into the Vlasov equation (16), we then reduce the Vlasov equation to the transport equations of the corresponding Fourier components a n , where the energy of quasi-particle is E p (τ, ρ) = m 2 (τ, ρ) + p 2 ⊥ cosh(ξ), with mass m(τ, ρ) obtained from the gap equation, which in the new coordinate becomes, The momentum integral in the gap equation only relates to the zeroth component a 0 , while the first and second order Fourier components contribute to the currents and energy-momentum tensor. The energy-momentum tensor is with the non-vanishing components are where p is the abbreviation for p ⊥ dp ⊥ dξ 4(2π) 2 . For this energy-momentum tensor, one can explicitly write down the flow velocity and energy density, being the time-lime eigenvector and the corresponding eigenvalue: Noting that the velocity has only two non-zero components, u τ and u ρ , the equilibrium distribution function (11) can be expressed as

IV. NUMERICAL PROCEDURE AND RESULT
In the numerical procedures, we solve the finite difference versions of transport equations. The distribution function is discretized on a fixed grid in the calculation frame. The phase space is discretized as follows, take 200 points for ρ within the range ρ/ρ 0 ∈ [−3, 3], 100 points for p T within p T /T 0 ∈ [0, 8], 100 points for ξ within the range ξ ∈ [0, 6]. The Fourier expansion of distribution function f with respect to θ is taken with maximum n = 7 to guarantee the convergence. To eschew the numerical instability around τ = 0, we take the initial time as τ 0 = 0.5 fm, the time step in the evolution is taken to be dτ = 0.0005 fm to guarantee the stability. The calculation at the discrete time step n+1 involves only quantities at the previous time step. At each time step, we solve both the transport equation and the gap equation, and eventually get the time and space dependence of the distribution function and the quark mass. This numerical framework is verified to be reliable by comparing the result with the analytical solution in a spherical symmetric system [20], as well as checking the conservation of the particle number and energy-momentum.
The initial state of the fireball is a highly off-equilibrium system, the evolution towards local equilibrium quark gluon plasma is also an interesting problem. However, this is not our concern in this paper, since the temperature is not well-defined in such initial stage, the discussion of phase transition is also questionable. We here discuss the evolution of the system from local equilibrium towards off-equilibrium state after the formation of quark gluon plasma. A local equilibrium initial state is adopted, and a local temperature could be assigned. The large gradient in the initial condition drives towards off-equilibrium state, while the collisions drives the system towards local equilibrium. In order to describe the hot chiral restored medium in the inner part and the cold chiral symmetry broken medium in the outer part, we choose a Gaussian temperature profile T (ρ) = T 0 exp −ρ 2 /ρ 2 0 for the initial state, with T 0 = 300 MeV and ρ 0 = 2 fm. For the local equilibrium initial state, the distribution function is the Fermi-Dirac distribution, a 0 (τ 0 , ρ, p ⊥ , ξ) = 2(e Ep/T (ρ) + 1) −1 , where energy is E p = m 2 + p 2 ⊥ cosh ξ, and T (ρ) is the above initial temperature profile. One can easily check that all other Fourier components vanish, a i (τ 0 ) = 0 for i ≥ 1. In the real case, the current mass is chosen to be m 0 = 3.7 MeV, and in the chiral limit, we have the current mass m 0 = 0.
The chiral phase transition is characterized by the chiral order parameter σ, or the quark constituent mass. The constituent mass is generated by the gap equation at each space-time point, and enters the transport equation through three ways: the energy E p = m 2 + p 2 ⊥ cosh ξ, the evolution rate of constituent mass ∂ τ m, and through the spatial gradient of mass ∂ ρ m. In order to illustrate the influence of the phase transition and the force term in the transport equation, we consider the following three different conditions. First, when solving the transport equation and the gap equation concurrently, the effect of the phase transition and the force terms are both taken into consideration. Second, for a comparison, we solve the transport equation alone and keep the quark mass as a constant, for instance m = 150 MeV. In this case, there is no phase transition nor force term. Third, in order to further illustrate the influence of the force term, we solve the transport equation and the gap equation at each step, but ignore the force term in the transport equation, namely assuming ∂ τ m = 0 and ∂ ρ m = 0. We also study the influence of out-ofequilibrium effect by comparing the results of different relaxation time, for small relaxation time the system stays close to local equilibrium; while for large relaxation time, the system is away from equilibrium. First, we self-consistently solve the coupled transport equation as well as the gap equation with finite current mass, and adopt different relaxation time. At each time step, the distribution function is obtained by evolving the transport equation, the energy density and entropy density are calculated by definition Eq. (??) and Eq. (21). The temperature is obtained by matching the realistic energy density to that of the equilibrium state. The energy density, entropy and temperature are presented in Fig. 2. The left panel corresponds to the evolution with small relaxation time, and the right panel corresponds to those of large relaxation time. The lines are rainbow colored which represents different evolution time, from the red line to the purple line represent the initial distribution to the distribution at later time.
The initial condition of the system is chosen to be an equilibrium distribution, with a gaussian distribution temperature profile, the inner part (small ρ) has higher temperature and the outer part (large ρ) has lower temperature. With the expansion of the system, the temperature of the core area gradually decrease, and the temperature of the outer area increases. The initial large gradient of the energy density drives the system away from the local equilibrium, while the collisions bring the system back to equilibrium. The entropy tells whether the system has reached local equilibrium, the solid line represent the realistic entropy density and the dashed lines represent the entropy of the equilibrium state. Since the equilibrium state takes the maximum entropy. When the collisions are not strong enough, the system takes longer time in the out-of-equilibrium state, where the realistic entropy is smaller than that of the equilibrium state. If the relaxation time is small, the collisions are strong enough to keep the system at local equilibrium, the entropy density of the state is the same as that of the equilibrium state. In the equilibrium state, the quark constituent mass m ψ serves as the order parameter of the chiral symmetry, which tells to what extent the chiral symmetry is broken. Here in an expanding quark-antiquark system, the space-time dependent quark mass m ψ (τ, ρ) signals the chiral symmetry in the space-time. In the initial state, the inner part of the system has higher temperature and is in the chiral symmetry restored phase, the constituent mass is small; the outer part of the system has lower temperature and is in the chiral symmetry broken phase with large constituent mass. By solving together the coupled transport equation and gap equation, the evolution of constituent quark mass in the time-space can be obtained self-consistently. The evolution of both the real case and chiral limit are investigated for large and small relaxation time, the results are presented in Fig. 3.
The upper two figures present the quark mass m ψ (τ, ρ) and transverse velocity v ρ (τ, ρ) ≡ (u ρ /u τ ) in the real case, the lower two figures correspond to those of chiral limit. With the expansion of the system, the quark mass in the inner area grows with time, indicating the gradually restoring of chiral symmetry; the quark mass of outer area decreases, indicating the breaking of chiral symmetry. In the equilibrium chiral phase transition, the phase transition point in the chiral limit is well-defined, while that of crossover does not has a strict definition, one of the usually used definition is the maximum of susceptibility dm ψ /dT . In the expanding system, the phase transition point in the chiral limit can still be defined by the time-space point where that quark mass reached zero. For the real case, we here take phase boundary as the space-time point where the equilibrium temperature is around the critical temperature and dm/dx takes the maximum. The phase transition points are marked out by the dotted lines in the Fig. 3. The phase diagram of equilibrium strong interaction matter in the T − µ plane is a map for the chiral symmetry. For non-equilibrium quark matter, the chiral symmetry breaking and restoration can be described by a phase diagram in space-time, with the phase boundary hypersurface indicates where the symmetry changes in the time-space. In the chiral limit, the phase transition hypersurface can still be defined by x µ (m ψ = 0) no matter whether the system is in the local equilibrium or not. Fig. 4 presents the hypersurface obtained by self-consistent solution of a chiral limit system, with different colored lines represent different relaxation time. The collisions cast the kinetic energy into internal energy, thus decelerates the expansion of the system. In the initial state (τ − τ 0 )/τ 0 = 0, the core area ρ/ρ 0 < 0.8 is in the chiral symmetry restored phase. With the expansion of the system, the core area gradually cools, the area of chiral symmetry restored phase shrinks. The free streaming system expands the fastest, after (τ − τ 0 )/τ 0 > 2.2, the chiral restored phase disappears. While the system with small relaxation time expands slower, it takes longer time for the chiral restored phase to disappear. It appears from the evolution of mass (Fig. 3) and the phase boundary (Fig. 4) that the various relaxation time does not have obvious influence on the quark mass. Since the order parameter describes the long-range correlation and the overall property of the system, while the collision is the local process in a system and is a short range correlation, the collision does not have big impact on the global symmetry of the system.

C. Kink in velocity
When solving together the transport equation and the gap equation, kinks in the velocity u τ and u ρ are discovered. The velocity along ρ-direction v ρ ≡ u ρ /u τ for different current mass and relaxation time are also presented in Fig. 3. As we have mentioned above, the quark mass enters the transport equation through both the energy E p and the force term ∇E p · ∇ p f . The gradient of the field energy acts as a continuous force on the quarks, describing the interaction between the quarks and the mean fields. This interaction changes the quarks' momenta. Although this force term is also considered in the simulations such as test particle method [18,19,[22][23][24], its influence on the phase transition and the expansion has not been carefully investigated. For a chiral phase transition at low density, it is either a crossover or a second order phase transition, in both case the order parameter changes continuously. In comparison, the gradient of order parameter diverges at a second order phase transition, hence the force could be extremely strong.
We first present the velocity v ρ in the scenario of constant mass so as to illustrate the influence of the phase transition on the evolution of the system, see Fig. 5. We take constant homogeneous mass m ψ (τ, ρ) = 150MeV and solve only the transport equation for a free streaming system, the force term vanishes since ∂ τ m = 0 and ∂ ρ m = 0. In this scenario, the kink does not appear, indicating the kink arises from the inhomogeneous mass distribution and thus the force term. In order to further illustrate the influence of the force term, we now concurrently solve the transport equation and the gap equation, but removes the force terms in the transport equation by fixing ∂ τ m = 0 and ∂ ρ m = 0. The evolution of mass and thermodynamic quantities has no obvious difference compared to those of self-consistent solutions, however, the velocity v ρ is quite different, which is presented as dashed lines in Fig. 6, with the solid lines are the velocity of self-consistent solution. As shown in each figures Fig. 6, when the force terms are ignored, the velocity has a bump around the phase transition. Namely if there is phase transition but no force term, the quark near the phase boundary are accelerated. In comparison, when the force term is considered, the quark near the phase boundary are slowed down. The appearance of kinks is closely related to the phase transition and the spatial distribution of the quark mass. The velocity kinks are more obvious in scenarios with small current mass or small relaxation time. This phenomenon can be understood as follows, since the force term is related to the susceptibility as well as the gradient of temperature ∇m = (dm/dT )∇T , the phase transition of the chiral limit has divergent susceptibility while that of real case is finite, thus the kinks in the chiral limit panels are more obvious compared with those in real case panel. The expansion starts with equilibrium distribution and gradually becomes out-of-equilibrium due to the huge pressure gradient. With large relaxation time, the system spends longer time in the out-of-equilibrium state, the critical effect is further washed out.

D. Distribution function
The phase transition hypersurface affects the motion of quarks around it, and has an influence on the p T spectrum of the quarks in the system. This can be directly revealed from the distribution function. Under the given symmetry in section III, the distribution function is defined on a 4 + 1d phase space, (ρ, p ⊥ , θ, ξ) as well as τ , where θ is the angle between the transverse coordinate ρ and the transverse momentum p T . The θ-dependence has been expanded to a series of Fourier coefficients. The zeroth component a 0 , which is obtained by integrating the distribution function f (ρ, p T , θ, ξ) over θ is related to the density distribution. While the first component a 1 obtained by integrating f (ρ, p T , θ, ξ) cos θ over θ gives hint to the direction of the particle velocity. In Fig. 7, we present the zeroth and the first Fourier components in the scenario with weak collisions, namely with large relaxation time τ θ = 10/T eq . The solid lines correspond the self-consistent distribution function, while the dashed lines are the distribution function where the force term is ignored. Coefficients a 0 and a 1 as a function of transverse momentum p T at various given transverse coordinates ρ/ρ 0 are presented in the left panel. The right panel shows the coefficients a 0 and a 1 as a function of transverse coordinates ρ/ρ 0 at various given transverse momentum p T .
At some evolution time τ = 1.0 fm or (τ − τ 0 )/τ 0 = 1, the phase transition takes place around the position ρ/ρ 0 ∼ 0.7, which is also the location of the kink in the velocity, see Fig. 3. From the left panel of Fig. 7, in the self-consistent solution (solid lines), the distribution function a 1 is negative at low p T for ρ/ρ 0 around 0.6 to 0.8, which means the particles with low momentum are bounced back by the phase transition "wall"; while the dashed lines reveals that without force term, the particles cannot see the "wall". This effect is also obvious from the a 1 in the right panel, for small momentums, the distribution function changes sign around ρ/ρ 0 = 0.7; for large momentums, the distribution function stays positive but has smaller values. This indicates that the particles with small momentum are bounced back by the phase transition wall around ρ/ρ 0 = 0.7, particles with large momentum go through the "wall", but have been slowed down. The integral of a 0 (ρ, p T , ξ) over p T at fixed ρ corresponds to the number density of particles somewhere in the transverse plane, while integral of a 0 (ρ, p T , ξ) over ρ at fixed p T is the number density of particles of some fixed momentum. It can be observed from the left panel, whether the force term is ignored or not, the number density away from the "wall" (ρ/ρ 0 = 0.7) is similar, while the force term would collect more particles around the "wall". From the right panel, there are more low momentum particles kept inside the wall because of the force term. FIG. 8: Zeroth and the first Fourier components in the case the collisions are strong τ θ = 0.005/T eq , the solid lines are the distribution function from the self-consistent solution, the dashed lines are the distribution function when the force term is ignored. The left panel shows a 0 and a 1 as functions of transverse momentum p T at various given transverse coordinates ρ/ρ 0 , the right panel shows a 0 and a 1 as functions of transverse coordinates ρ/ρ 0 at various given transverse momentum p T .
As is presented in Fig. 8, the effect of the phase transition "wall" is smoothed by the collision. Since the particles are relaxed into thermal distribution, the direction of a single particle is aligned along the collective velocity. With strong collisions, the influence of the force term on number density almost disappear, giving the same a 0 whether the force term is taken into consideration. The first order components a 1 is still affected by the force term. Inside the "wall" (ρ/ρ 0 < 0.7), the force term does not has any influence. Outside the "wall" (ρ/ρ 0 > 0.7), the distribution a 1 is smaller when the force term is considered, which means the particles are slowed down by the force term.

V. SUMMARY
Quite different from the equilibrium thermodynamics, the realistic chiral phase transition in heavy ion collision takes place in a highly inhomogeneous, fast evolving dynamical system. This requires the understanding of the chiral phase transition in an expanding, out-of-equilibrium system.
In this work, we investigate the evolution of an expanding quark-antiquark system with self-consistently taking into account the dynamical quark constituent mass. In order to reduce the dimension of the phase space, we con-sider a longitudinal boost invariant and transversal rotational symmetric system, which is a good approximation for ultra-central heavy ion collisions. In the numerical process, both Vlasov and gap equations are solved concurrently, giving a self-consistent evolution of both the quark-antiquark distribution function and the quark constituent mass. The spacetime-dependent constituent mass serves as the chiral order parameter and affects the evolution of the quark distribution function through the force term. In order to investigate the off-equilibrium effects, we introduce relaxation time approximation for the collision term, and compare the local equilibrium and out-of-equilibrium result by considering small and large relaxation time.
The evolution of the quark mass illustrates the chiral phase transition, and defines the phase diagram in the spacetime. A kink in the transverse velocity in observed around the phase transition boundary, which appears because the large force term around the phase transition. The kink is more obvious for smaller current mass and smaller relaxation time, which means that the crossover and non-equilibrium effect tends to smooth out the kink. The influence of the phase transition hypersurface is further investigated by directly analyzing the distribution function. It is observed that, the phase transition wall would bounce back the low momentum particles, thus gives a kink in the velocity and may enhance the low p T part of the momentum spectrum. The kink indicates that the parton decelerates when propagating from hot to cold region because of the increase in its effective mass. This effect may have impact on the produced hadrons and other observables in experiment, such as reducing the mean p T and the enhancement of particle production yield, which will be investigated in future work.