Entropy Production in Dissipationless Hydrodynamics with an order parameter

It is remarkable that relativistic hydrodynamics provides accurate description of bulk evolution of matter produced in heavy ion collisions, which consists of around a thousand particles [1, 2]. With the assumption of local equilibrium, the problem of complicated many particle dynamics is reduced to the conservation of energy, momentum and baryon number as equation of motion for relativistic hydrodynamics. Over the past decade, the framework of relativistic hydrodynamics has been furnished in many aspects for phenomenological application in heavy ion collisions: the inclusion of viscous correction leads to more accurate description of the bulk evolution [3, 4, 5]; the inclusion of noises allows for systematic treatment of fluctuations [6, 7, 8]; the inclusion of particle momentum anisotropy extends the regime of applicability to earlier time [9, 10] etc.

Recently the beam energy scan (BES) program in relativistic heavy ion collider (RHIC) is devoted to pinning down possible existence of critical point in the phase diagram of quantum chromodynamics (QCD) [11]. Hydrodynamic studies of system evolution close to the critical point necessitates the inclusion of critical mode into the present framework of hydrodynamics. This has been pursued both by different group [12, 13, 14, 15]. The coupling of critical mode and the hydrodynamic modes is found to alter the bulk evolution, in particular near the phase transition, see [16] for a recent review.

A key quantity in heavy ion collisions is the entropy production. It is known that in complicated dynamics of fireball evolution, entropy is produced in different stages of the collisions, see [17] and references therein. A significant amount of entropy is produced during hydrodynamic evolution due to dissipation of the hydrodynamic mode. In the presence of critical mode, the dissipation of critical mode can also lead to entropy production [18, 19]. The goal of this paper is to study entropy production when dissipation of either mode is absent. We will show distinct features of boost invariant solution and non-boost invariant solution, with the former preserving entropy and the latter not.

The paper is organized as follows: In Section 2, we describe the hydrodynamics with critical mode based on linear sigma model. In Section 3, we present numerical results for boost invariant and non-boost invariant solutions and discuss physical implications. We give conclusion and outlook in Section 4.

We start with the Lagrangian of linear sigma model [20]:

$\displaystyle{\cal L}={\bar{q}}\left(i{\not{\partial}}-g({\sigma}+i{\gamma}_{5}{\vec{\tau}}{\vec{\pi}})\right)q+\frac{1}{2}\left((\partial{\sigma})^{2}+(\partial{\vec{\pi}})^{2}\right)-U({\sigma},{\vec{\pi}}),$

(1)

with $q$, ${\sigma}$ and ${\vec{\pi}}$ being the quark, sigma and pions fields respectively. The condensation of ${\sigma}$ gives mass to quarks $M_{q}=g\langle{\sigma}\rangle$, breaking chiral symmetry. The symmetry is broken through the potential $U$ given by

$\displaystyle U({\sigma},{\vec{\pi}})=\frac{{\lambda}}{4}\left({\sigma}^{2}+{\vec{\pi}}^{2}-v^{2}\right)^{2}-c{\sigma}.$

(2)

Throughout the paper, we use mean-field approximation for ${\sigma}$ and ${\vec{\pi}}$. In the absence of isospin chemical potential, pions do not condense thus $\langle{\vec{\pi}}\rangle=0$. $\langle{\sigma}\rangle$ is the only order parameter in the model. This is the critical mode to be included in hydrodynamics. The parameters in (1) and (2) are fixed as

$\displaystyle g=M_{q}/f_{\pi},\quad c=M_{\pi}^{2}f_{\pi},\quad{\lambda}=\frac{1}{2f_{\pi}^{2}}\left(M_{\sigma}^{2}-M_{\pi}^{2}\right),\quad v^{2}=f_{\pi}^{2}-M_{\pi}^{2}/{\lambda},$

(3)

with the experimental input $M_{\pi}=138\text{MeV}$, $M_{\sigma}=600\text{MeV}$, $f_{\pi}=93\text{MeV}$. The condensate $\langle{\sigma}\rangle$ is determined dynamically by minimizing the thermodynamic potential ${\Omega}=U+{\Omega}_{q{\bar{q}}}$, with the quark contribution ${\Omega}_{q{\bar{q}}}$ given by

$\displaystyle{\Omega}_{q{\bar{q}}}(T,{\mu})={\nu}_{q}T\int\frac{d^{3}k}{(2{\pi})^{3}}\big{[}\ln\left(1-n_{q}(T,{\mu},k)\right)+\ln\left(1-n_{\bar{q}}(T,{\mu},k)\right)\big{]}.$

(4)

Here ${\nu}_{q}=2N_{c}N_{f}=12$ counts the spins, colors and flavors of the quark field. $n_{q}$ and $n_{\bar{q}}$ are Fermi-Dirac distributions for quark and anti-quark respectively:

$\displaystyle n_{q}(T,{\mu},k)=\frac{1}{e^{\left(\sqrt{k^{2}+M_{q}^{2}}-{\mu}\right)/T}+1},\quad n_{\bar{q}}(T,{\mu},k)=n_{q}(T,-{\mu},k).$

(5)

Following [12, 13], we treat quark fields as hydrodynamic degree of freedom, which is coupled to sigma mean field. The equations of motion (EOM) for ${\sigma}$ and hydrodynamic fields are given by:

		$\displaystyle D_{\mu}D^{\mu}{\sigma}+\frac{{\delta}{\Omega}}{{\delta}{\sigma}}=0,$		(6)
		$\displaystyle D_{\mu}\left(T_{q}^{{\mu}{\nu}}+T_{\sigma}^{{\mu}{\nu}}\right)=0.$		(7)

For simplicity, we choose to work with vanishing quark number, so there is no additional equation for quark number conservation and ${\mu}_{B}=0$. Note that we choose to write the EOM in curved spacetime with $D_{\mu}$ denoting covariant derivative. This will be useful for adapting to Milner coordinates in the next section. $T_{q}^{{\mu}{\nu}}=({\epsilon}+p)u^{\mu}u^{\nu}-pg^{{\mu}{\nu}}$ and $T_{\sigma}^{{\mu}{\nu}}=\partial^{\mu}{\sigma}\partial^{\nu}{\sigma}-g^{{\mu}{\nu}}\left(\frac{1}{2}(\partial{\sigma})^{2}-U({\sigma})\right)$ correspond to contributions to stress energy tensor from quark and sigma fields respectively. Using (6), we can express the divergence of $T_{\sigma}^{{\mu}{\nu}}$ field as

$\displaystyle D_{\mu}T_{\sigma}^{{\mu}{\nu}}=-\frac{{\delta}{\Omega}_{q{\bar{q}}}}{{\delta}{\sigma}}\partial^{\nu}{\sigma}.$

(8)

We can thus rewrite conservation of stress energy tensor as

$\displaystyle D_{\mu}T_{q}^{{\mu}{\nu}}=\frac{{\delta}{\Omega}_{q{\bar{q}}}}{{\delta}{\sigma}}\partial^{\nu}{\sigma}.$

(9)

This has clear interpretation that the stress tensor from the quark field is conserved up to work and force by the sigma field. We will solve (6) and (9) numerically with the hydrodynamic fields $T(x)$, $u^{\mu}(x)$ and mean field ${\sigma}$ in the next section.

Before presenting numerical results, it is instructive to see the implication of (6) and (9) for entropy production. To this end, we define entropy at ${\mu}=0$ from the following form

$\displaystyle dp=sdT+\frac{\partial p}{\partial{\sigma}}d{\sigma}\;\;\Rightarrow s=\frac{\partial p}{\partial T}|_{\sigma}.$

(10)

To derive the evolution of $s$, we contract (9) with $u_{\nu}$ and use $u^{2}=1$ to obtain

$\displaystyle\left({\epsilon}+p\right)D_{\mu}u^{\mu}+u^{\mu}\partial_{\mu}{\epsilon}+\left({\epsilon}+p\right)u^{\mu}u^{\nu}D_{\mu}u_{\nu}=\frac{{\delta}{\Omega}_{q{\bar{q}}}}{{\delta}{\sigma}}u_{\nu}\partial^{\nu}{\sigma}.$

(11)

The last term on the left hand side can be dropped by using $u^{\nu}D_{\mu}u_{\nu}=0$. Using ${\epsilon}+p=Ts$ and (10), we obtain

$\displaystyle d{\epsilon}=Tds-\frac{\partial p}{\partial{\sigma}}d{\sigma}.$

(12)

Plugging (12) into (11), we obtain

$\displaystyle D_{\mu}\left(su^{\mu}\right)=0.$

(13)

This is the conservation of entropy current. Note that the $d{\sigma}$ terms in (10) and (12) corresponding to change of quark mass cancel each other in $s$. In the absence of order parameter, the existence of entropy current in ideal hydrodynamics is expected from the absence of dissipation in the EOM. (13) shows the presence of order parameter does not affect the existence of entropy current with our definition of entropy density.

In this section, we solve (6) and (9) numerically. For comparison, we first study solutions with boost invariance, and then solutions without. In both cases, homogeneity in transverse plane is assumed. In terms of proper time ${\tau}=\sqrt{t^{2}-z^{2}}$ and spacetime rapidity ${\eta}=\tanh^{-1}\frac{z}{t}$, the explicit form of EOM are given by

		$\displaystyle\partial_{\tau}^{2}{\sigma}+\frac{1}{{\tau}}\partial_{\tau}{\sigma}-\frac{1}{{\tau}^{2}}\partial_{\eta}^{2}{\sigma}+\frac{{\delta}{\Omega}}{{\delta}{\sigma}}=0,$		(14)
		$\displaystyle\partial_{\tau}T_{q}^{{\tau}{\tau}}+\frac{1}{{\tau}}T_{q}^{{\tau}{\tau}}+\partial_{\eta}T_{q}^{{\tau}{\eta}}+{\tau}T_{q}^{{\eta}{\eta}}-\frac{{\delta}{\Omega}_{q{\bar{q}}}}{{\delta}{\sigma}}\partial_{\tau}{\sigma}=0,$		(15)
		$\displaystyle\partial_{\tau}T_{q}^{{\tau}{\eta}}+\frac{1}{{\tau}}\partial_{\eta}T_{q}^{{\tau}{\eta}}+\partial_{\eta}T_{q}^{{\eta}{\eta}}+\frac{2}{{\tau}}T_{q}^{{\tau}{\eta}}+\frac{1}{{\tau}^{2}}\frac{{\delta}{\Omega}_{q{\bar{q}}}}{{\delta}{\sigma}}\partial_{\eta}{\sigma}=0.$		(16)

In the boost invariant case, the flow is fixed as $u^{\mu}=(1,0,0,0)$, and the last equation in (14) is trivially satisfied with $T$ and ${\sigma}$ being functions of ${\tau}$ only. To solve (14), we use the following initial conditions:

$\displaystyle T({\tau}={\tau}_{0})=T_{0},\quad{\sigma}({\tau}={\tau}_{0})={\sigma}_{0},\quad\partial_{\tau}{\sigma}({\tau}={\tau}_{0})=0.$

(17)

We fix ${\sigma}_{0}$ to be the equilibrium value of ${\sigma}$ at $T_{0}$. In Fig. 1, we show the ${\tau}$ dependencies of $T$ and ${\sigma}$. As temperature drops due to the expansion, the sigma field rises in an oscillatory fashion. The oscillation persists at very late time in the absence of dissipation. Similar behavior is also observed in [18]. Since the dynamics of ${\sigma}$ is coupled with $T$, oscillation is also present in $T({\tau})$. Fig. 1 also indicates approximate scaling for temperature ${\sigma}$: $T\sim{\tau}^{-{\alpha}}$. We extract the exponent ${\alpha}$ as a function of the initial temperature $T_{0}$. In Fig. 2, we show the $T_{0}$ dependence of ${\alpha}$. ${\alpha}(T_{0})$ is an increasing function, indicating that the dropping of the temperature is faster for higher initial temperature. In the high temperature limit, ${\alpha}\to\frac{1}{3}$. This is an expected limit, as chiral symmetry is essentially restored with ${\sigma}\simeq 0$. The scaling of the temperature is given by Bjorken solution with $T\sim{\tau}^{-c_{s}^{2}}$ and $c_{s}^{2}=\frac{\partial p}{\partial{\epsilon}}\simeq\frac{1}{3}$.

The total entropy in the boost invariant case can be shown to be strictly conserved. To see that, we write out the explicit form of (13) as

$\displaystyle\partial_{\tau}s+\frac{1}{{\tau}}s=0,\Rightarrow\frac{1}{{\tau}}\partial_{\tau}\left({\tau}s\right)=0.$

(18)

Note that ${\tau}sd{\eta}d^{2}x_{\perp}$ is the entropy per unit volume in Milner coordinates. It follows that entropy is strictly conserved. We also verify this numerically from the definition of entropy (10).

Now we move on to solution without boost invariance. In the absence of order parameter, Bjorken solution is supposed to describe well bulk evolution at mid-rapidity. We will see this is not quite true when order parameter is present. Without boost invariance, longitudinal flow velocity $u^{\tau}$ becomes an independent dynamical field. We solve (14) in the interval ${\eta}=[0,{\eta}_{m}=10]$ with the following initial conditions

		$\displaystyle T({\tau}_{0},{\eta})=(T_{0}-T_{\text{min}})\left(e^{-{\eta}^{2}/{\Delta}^{2}}-e^{-{\eta}_{m}/{\Delta}^{2}}\right)+T_{\text{min}},$
		$\displaystyle{\sigma}({\tau}_{0},{\eta})=({\sigma}_{0}-{\sigma}_{\text{max}})\left(e^{-{\eta}^{2}/{\Delta}^{2}}-e^{-{\eta}_{m}/{\Delta}^{2}}\right)+{\sigma}_{0},$
		$\displaystyle\partial_{\tau}{\sigma}({\tau}_{0},{\eta})=0,$
		$\displaystyle u^{\tau}({\tau}_{0},{\eta})=(u_{\text{min}}-u_{\text{max}})\left(e^{-{\eta}^{2}/{\Delta}^{2}}-e^{-{\eta}_{m}/{\Delta}^{2}}\right)+u_{\text{max}}.$		(19)

The initial profile of $T$, ${\sigma}$ and $u^{\tau}$ are taken to be a Gaussian form with a variance ${\Delta}^{2}=5$. The initial maximum/minimum temperatures are at ${\eta}=0$ and ${\eta}={\eta}_{m}$ respectively. ${\sigma}_{0}$ and ${\sigma}_{\text{max}}$ are taken to be equilibrium values for the corresponding temperatures. Motivated by the Bjorken solution, we take $u_{\text{min}}=1$. Early works on longitudinal dynamics choose $u_{\text{max}}=1$ [21, 22], see also [23, 24] with Landau initial condition. In our study, keeping $u_{\text{max}}\neq 1$ is necessary for the stability of solutions. We take $u_{\text{max}}$, together with $T_{0}$, $T_{\text{min}}$ as free parameters. With the finite interval in ${\eta}$, we still need boundary conditions, for which we take

		$\displaystyle T({\tau},{\eta}_{m})=T_{\text{min}},$
		$\displaystyle{\sigma}({\tau},{\eta}_{m})={\sigma}_{\text{max}},$
		$\displaystyle\partial_{\eta}{\sigma}({\tau},0)=0,$
		$\displaystyle u^{\tau}({\tau},0)=u_{\text{min}}.$		(20)

We impose the minimum temperature $T_{\text{min}}$ and equilibrium value of ${\sigma}$ at ${\eta}={\eta}_{m}$. The third condition in (3) follows from the symmetry: ${\eta}\leftrightarrow-{\eta}$. We show in Fig. 3 evolution of ${\sigma}$, $T$ for $T_{0}=240\text{MeV}$, $T_{\text{min}}=40\text{MeV}$ and $u_{\text{max}}=2$. Similar to the case with boost invariance, ${\sigma}$ rises with oscillations. $T$ drops with oscillations correlated with those in ${\sigma}$. As indicated in Fig. 3, the drop of temperature shows an approximate scaling law: $T\sim{\tau}^{-{\alpha}}$ at ${\eta}=0$. However, the exponent ${\alpha}$ takes different values from those in boost invariant case starting with the same $T_{0}$, as shown in Fig. 4. The case without boost invariance shows faster drop in temperature. This is probably caused by extra longitudinal flow, which takes away energy from the higher temperature mid-rapidity region. The behavior of $u^{\tau}$ is unique in the non-boost invariant case. The flow velocity near ${\eta}={\eta}_{m}$ quickly drop to zero. The region near ${\eta}=5$ has persistent extra longitudinal flow. While we impose $u^{\tau}=u_{\text{max}}$ at ${\eta}={\eta}_{m}$ at initial time, we find it quickly drop to $1$ as ${\tau}$ increases, indicating the longitudinal flow quickly approaches Bjorken flow. At around ${\eta}=5$, the extra longitudinal flow persists at late time.

The behavior of entropy is qualitatively different from the boost invariant case. Integrating (13) with $\int{\tau}d{\tau}d{\eta}d^{2}x_{\perp}$, we have

$\displaystyle 0=\int{\tau}d{\tau}d{\eta}D_{\mu}\left(su^{\mu}\right)=\int d{\tau}d{\eta}d^{2}x_{\perp}\big{[}\partial_{\tau}\left({\tau}su^{\tau}\right)+\partial_{\eta}\left({\tau}su^{\eta}\right)\big{]}.$

(21)

The second term on the right hand side (21) is a boundary term. If we ignore it for the moment, the conserved quantity is ${\tau}su^{\tau}$ rather than ${\tau}s$, which is the entropy density per rapidity per transverse area. Therefore we expect entropy is not conserved in this case. Fig. 3 suggests the rapidity averaged $u^{\tau}$ decreases with ${\tau}$, thus we have entropy production in non-boost invariant solution. In Fig. 5, we confirm this by direct numerical integration of $\int d{\eta}{\tau}s$ for different initial $T_{0}$ and $u_{\text{max}}=2$. The expression $\int d{\eta}{\tau}s$ measures total entropy per transverse area along constant proper time slice. This expression depends on our choice of time. For application to fireball evolution in heavy ion collisions, we may also use isothermal contours for integration. This definition is phenomenologically more relevant if we assume the freezeout temperature of the fireball is a constant. Fig. 6 shows isothermal for one particular solution with $T_{0}=280\text{MeV}$. On each isothermal contour, the total entropy per transverse area is given by

$\displaystyle\int d\Sigma{\tau}s=\int{\tau}d{\eta}\left(u^{\tau}-u^{\eta}\frac{\partial{\tau}}{\partial{\eta}}\right).$

(22)

We show total entropy per transverse area at different contour times in Fig. 6. A more pronounced increase with contour time is seen.

The increase of entropy in non-boost invariant case can be understood as follows: while the system we consider has no dissipation, it evolves starting from an off-equilibrium state, with inhomogeneity in rapidity direction. The evolution smoothens out the inhomogeneity and is irreversible. It is analogous to irreversible expansion of ideal gas in free space. This process generates net entropy. For the same reason, the boost invariant solution preserves entropy due to the absence of inhomogeneity.

We studied hydrodynamics with order parameter based on linear sigma model. We obtained numerical solutions both with and without boost invariance. In both cases, we found similar behavior for the order parameter ${\sigma}$ and temperature $T$, with ${\sigma}$ rises in an oscillatory fashion. A correlated behavior is also seen in the temperature: $T$ oscillates around a power law decrease $T\sim{\tau}^{-{\alpha}}$, as shown in Fig. 1 and Fig. 3. However, the exponent for boost invariant solution is smaller than the counterpart of non-boost invariant solution with the same initial temperature at ${\eta}=0$. This is probably due to the extra longitudinal flow in the non-boost invariant case taking away energy from the high temperature mid-rapidity region.

We also studied entropy production for the solutions. The boost invariant solution is shown to preserve entropy in the absence of dissipation. However this is no longer true for non-boost invariant solution. We used two definitions of total entropy, one from integration on constant proper time and the other on isothermal contour, with both showing entropy production. This can be understood as the evolution smoothening out the inhomogeneity of the off-equilibrium state.

While we only studied hydrodynamics at ${\mu}_{B}=0$, where the phase transition is a crossover, it is more interesting to study the hydrodynamics with a first order phase transition. In the latter case, a smooth sigma profile for the initial condition would necessarily be non-monotonic, which contains larger inhomogeneity. By the inhomogeneity induced mechanism discussed in this paper, it is tempting to expect a larger entropy production. We leave more quantitative studies for the future.

S.L. is grateful to Huichao Song for useful discussions. This work is in part supported by NSFC under Grant Nos 11675274 and 11735007.