Exact solution for the non-equilibrium attractor in number-conserving relaxation time approximation

TL;DR

This paper derives an exact numerical solution for the non-equilibrium attractor in a conformal kinetic theory with number conservation, showing faster convergence to the attractor when chemical potential dynamics are included.

Contribution

It introduces a dynamical fugacity to enforce number conservation and derives coupled integral equations for temperature and fugacity, providing the first exact solution for this system.

Findings

Solutions converge to a unique non-equilibrium attractor.

Number conservation accelerates the convergence of moments.

Comparison with hydrodynamic models shows good agreement.

Abstract

We extend previous studies of the conformal 0+1d kinetic non-equilibrium attractor in relaxation time approximation by enforcing number conservation through the introduction of a dynamical fugacity (chemical potential). We derive two coupled integral equations for the effective temperature and fugacity which are then solved numerically to obtain the exact solution. The resulting solutions exhibit convergence to a unique non-equilibrium attractor when the scaled moments of the distribution function are plotted as a function of the rescaled time w = tau/tau_eq. This occurs even though the system is out of chemical equilibrium at late times. In addition, compared to the case where number conservation was not imposed, we find that the moments converge to their respective attractors more quickly, particularly for moments with m=0. Finally, we compare the resulting attractor moments with predictions from different hydrodynamic frameworks.

Equations

$$p^{\mu }{\partial }_{\mu }f(x,p)=C[f(x,p)],$$

$$C[f]=-\frac{p\cdot u}{{\tau }_{\mathrm{eq}}}\left(f-f_{\mathrm{eq}}\right).$$

$$f_{\mathrm{eq}}(\tau ,p)=\Gamma (\tau )\exp \left(-\frac{p\cdot u(\tau )}{T(\tau )}\right).$$

$$n_{\mathrm{eq}}=\frac{\Gamma T^3}{{\pi }^2},{\mathcal{S}}_{\mathrm{eq}}=\frac{4\Gamma T^3}{{\pi }^2},$$

$${\mathcal{E}}_{\mathrm{eq}}=\frac{3\Gamma T^4}{{\pi }^2},{\mathcal{P}}_{\mathrm{eq}}=\frac{\Gamma T^4}{{\pi }^2}.$$

$$w\equiv t{p}_z-zE,$$

$$v\equiv Et-p_{z}z=\sqrt{w^2+\left(m^2+{\vec{p}}_T^2\right){\tau }^2}.$$

$$E=p^0=\frac{vt+wz}{{\tau }^2},p_z=\frac{wt+vz}{{\tau }^2}.$$

$$dP=\frac{d^{4}p}{(2\pi {)}^4}2\pi \delta \left(p^2-m^2\right)2\theta (p^0)=\frac{d{p}_z}{(2\pi {)}^3{p}^0}{d}^2{p}_T=\frac{dw{d}^2{p}_T}{(2\pi {)}^{3}v}.$$

$$\frac{\partial f(\tau ,w,p_T)}{\partial \tau }$$

$$f_{\mathrm{eq}}(\tau ,w,p_T)=\Gamma (\tau )\exp \left[-\frac{\sqrt{w^2+p_T^2{\tau }^2}}{T(\tau )\tau }\right].$$

$$f(\tau ,w,p_T)=D(\tau ,{\tau }_0)f_0(w,p_T)+{\int }_{{\tau }_0}^{\tau }\frac{d{\tau }^{\prime }}{{\tau }_{\mathrm{eq}}({\tau }^{\prime })}D(\tau ,{\tau }^{\prime })f_{\mathrm{eq}}({\tau }^{\prime },w,p_T),$$

$$D({\tau }_2,{\tau }_1)=\exp \left[-\underset{{\tau }_1}{\overset{{\tau }_2}{\int }}\frac{d{\tau }^{\prime \prime }}{{\tau }_{\mathrm{eq}}({\tau }^{\prime \prime })}\right].$$

$${\tau }_{\mathrm{eq}}(\tau )=\frac{5\bar{\eta }}{T(\tau )},$$

$$f_0(w,p_T)={\gamma }_0\exp \left[-\frac{\sqrt{(1+{\xi }_0)w^2+p_T^2{\tau }_0^2}}{{\Lambda }_0{\tau }_0}\right]={\gamma }_0\exp \left[-\frac{\sqrt{w^2{\alpha }_0^{-2}+p_T^2{\tau }_0^2}}{{\Lambda }_0{\tau }_0}\right],$$

$$T$$

$$\Gamma$$

$${\mathcal{H}}^{nm}(y)\equiv \frac{2y^{2m+1}}{2m+1}{}_2{F}_1(\frac{1}{2}+m,\frac{1-n}{2};\frac{3}{2}+m;1-y^2).$$

$${\mathcal{M}}^{nm}[f]\equiv \int dP(p\cdot u{)}^n(p\cdot z{)}^{2m}f(\tau ,w,p_T).$$

$$\mathcal{E}={\mathcal{M}}^{20}=\int dP(p\cdot u{)}^{2}f(\tau ,w,p_T)=T_{\mathrm{LRF}}^{00},$$

$$n={\mathcal{M}}^{10}=\int dP(p\cdot u)f(\tau ,w,p_T),$$

$${\mathcal{M}}^{nm}[f]=\frac{1}{(2\pi {)}^3{\tau }^{n+2m}}\int dw{d}^2{p}_T{v}^{n-1}{w}^{2m}f(\tau ,w,p_T)$$

$${\mathcal{M}}^{nm}(\tau )=D(\tau ,{\tau }_0){\mathcal{M}}_0^{nm}(\tau )+{\int }_{{\tau }_0}^{\tau }\frac{d{\tau }^{\prime }}{{\tau }_{\mathrm{eq}}({\tau }^{\prime })}D(\tau ,{\tau }^{\prime }){\mathcal{M}}_{\mathrm{eq}}^{nm}(\tau ,{\tau }^{\prime }),$$

$${\mathcal{M}}_0^{nm}(\tau )$$

$${\mathcal{M}}_{\mathrm{eq}}^{nm}(\tau ,{\tau }^{\prime })$$

$${\mathcal{M}}^{nm}(\tau )$$

$$\begin{gathered} \Gamma (\tau )T^4(\tau )=D(\tau ,{\tau }_0){\Gamma }_0{T}_0^4\frac{\mathcal{H}\left(\frac{{\alpha }_0{\tau }_0}{\tau }\right)}{\mathcal{H}({\alpha }_0)} \\ +{\int }_{{\tau }_0}^{\tau }\frac{d{\tau }^{\prime }}{2{\tau }_{\mathrm{eq}}({\tau }^{\prime })}D(\tau ,{\tau }^{\prime })\Gamma ({\tau }^{\prime })T^4({\tau }^{\prime })\mathcal{H}\left(\frac{{\tau }^{\prime }}{\tau }\right), \end{gathered}$$

$$\Gamma (\tau )T^3(\tau )=\frac{1}{\tau }\left[D(\tau ,{\tau }_0){\Gamma }_0{T}_0^3{\tau }_0+{\int }_{{\tau }_0}^{\tau }\frac{d{\tau }^{\prime }}{{\tau }_{\mathrm{eq}}({\tau }^{\prime })}D(\tau ,{\tau }^{\prime })\Gamma ({\tau }^{\prime })T^3({\tau }^{\prime }){\tau }^{\prime }\right],$$

$${\overline{\mathcal{M}}}^{nm}(\tau )\equiv \frac{{\mathcal{M}}^{nm}(\tau )}{{\mathcal{M}}_{\mathrm{eq}}^{nm}(\tau )},$$

$${\mathcal{M}}_{\mathrm{eq}}^{nm}(\tau )={\mathcal{M}}_{\mathrm{eq}}^{nm}(\tau ,\tau )=\frac{(n+2m+1)!\Gamma (\tau )T^{n+2m+2}(\tau )}{2{\pi }^2(2m+1)},$$

$${\overline{\mathcal{M}}}_{\mathrm{vHydro}}^{nm}(\tau )=1-\frac{3m(n+2m+2)(n+2m+3)}{4(2m+3)}\bar{\pi }.$$

$${\overline{\mathcal{M}}}_{\mathrm{aHydro}}^{nm}(\tau )=(2m+1)(2\alpha {)}^{n+2m-2}\frac{{\mathcal{H}}^{nm}(\alpha )}{[\mathcal{H}(\alpha ){]}^{n+2m-1}}.$$

Full text

1 Introduction

Based on theory to data comparisons produced over the course of the last decades, there is now a strong body of evidence that the dynamics of the quark-gluon plasma (QGP) created in ultrarelativistic heavy-ion collisions (URHICs) is well described by relativistic dissipative hydrodynamics Averbeck:2015jja ; Jeon:2016uym ; Romatschke:2017ejr ; Florkowski:2017olj ; Alqahtani:2017mhy . Since the phenomenologically extracted values of the shear viscosity to entropy density ratio are finite, this implies that at early times after the nuclear pass through ($0.01-0.1\lesssim\tau\lesssim 1$ fm/c), the QGP possesses large non-equilibrium corrections, e.g. large pressure anisotropy in the local rest frame. Because the system experiences large deviations from local thermal equilibrium, one might expect dissipative hydrodynamics approaches to fail at early times. In practice, however, one finds that dissipative hydrodynamics describes the evolution of the components of the energy-momentum tensor quite well after a rather short time scale $\tau_{\rm hydro}\sim 0.5-1$ fm/c in the center of the overlap region for a central collision. Since dissipative hydrodynamics frameworks perform well after $\tau\sim\tau_{\rm hydro}$, the system is said to hydrodynamize at this time scale Chesler:2008hg . The time scale for hydrodynamization has been extracted by comparing numerical solutions of the underlying microscopic dynamical equations to dissipative hydrodynamics evolution in both the weak and strong coupling limits Chesler:2008hg ; Heller:2013oxa ; Keegan:2015avk ; Romatschke:2017vte ; Strickland:2017kux . From these studies one finds that $\tau_{\rm hydro}\sim 2/T$. At the highest LHC energies and assuming $\eta/s=0.2$, this translates into $\tau_{\rm hydro}\sim 0.5$ fm/c when considering the center of the fireball created in a zero impact parameter collision.

The fact that the system is quickly driven towards dissipative hydrodynamical evolution can be understood using the concept of the hydrodynamical attractor Heller:2015dha . In 0+1d conformal viscous hydrodynamics, one can reduce the two coupled equations for the energy density and the shear pressure correction to a single ordinary differential equation which, subject to the correct boundary conditions, provides a universal “attractor” solution for the scaled shear correction $\bar{\pi}=\pi^{\eta}_{\eta}/\varepsilon$ as a function of the scaled time $\bar{w}=\tau/\tau_{\rm eq}$, for example. If one solves the hydrodynamic equations with different initial conditions and plots the results versus $\bar{w}$, one finds that the solutions with different initial conditions converge to the universal attractor solution on a very short time scale (in the sense of small $\bar{w}$). This observation is not restricted to second-order viscous hydrodynamics and has been shown to hold in numerical solutions to Einstein’s equations obtained in the strong coupling limit of ${\cal N}=4$ supersymmetric Yang-Mills in the large $N$ limit Heller:2013oxa ; Keegan:2015avk ; Romatschke:2017vte , QCD effective kinetic theory simulations Keegan:2015avk ; Kurkela:2015qoa ; Romatschke:2017vte , third-order viscous hydrodynamics JaiswalForth , anisotropic hydrodynamics Strickland:2017kux , and exact solutions to the Boltzmann equation in relaxation time approximation (RTA) subject to both Bjorken and Gubser flows Romatschke:2017vte ; Strickland:2017kux ; Behtash:2017wqg ; Behtash:2018moe ; Denicol:2018pak ; Strickland:2018ayk ; Behtash:2019txb .

Recently, using the exact solution of the RTA Boltzmann equation subject to Bjorken flow, it was demonstrated that the idea of the non-equilibrium attractor can be extended beyond the low-order moments of the one-particle distribution function typically considered in hydrodynamic approaches Strickland:2018ayk . In Ref. Strickland:2018ayk it was demonstrated that the full one-particle distribution exhibits attractor-like behavior and that higher moments, ${\cal M}^{nm}$ of the one-particle distribution function converge more quickly to their respective attractors, with the exception being moments with $m=0$, which are more sensitive to the squeezed free-streaming part of the exact solution. For moments with large $m$ and $n$, Ref. Strickland:2018ayk showed that there is a parametrically large separation between the scaled time at which solutions converge to the non-equilibrium attractor $\bar{w}_{c}$ and the time at which the moment approaches to within 10% of its equilibrium value $\bar{w}_{\rm therm}$. Finally, In Ref. Strickland:2018ayk comparisons were made between the exact attractor moments and various dissipative hydrodynamics frameworks including relativistic Navier-Stokes (NS) Eckart:1940te ; ldlandau2013 ; Weinberg:1971mx , second order viscous hydrodynamics Muller:1967zza ; Israel:1976tn ; Israel:1979wp ; Muronga:2001zk ; Muronga:2003ta ; Muronga:2004sf ; Heinz:2005bw ; Baier:2006um ; Romatschke:2007mq ; Baier:2007ix ; Dusling:2007gi ; Luzum:2008cw ; Song:2008hj ; Heinz:2009xj ; Schenke:2010rr ; Schenke:2011tv ; Bozek:2011wa ; Niemi:2011ix ; Denicol:2011fa ; Niemi:2012ry ; Bozek:2012qs ; Denicol:2012cn ; Denicol:2012es ; Denicol:2014vaa ; Denicol:2014mca ; Jaiswal:2014isa , third-order viscous hydrodynamics Jaiswal:2013npa ; Jaiswal:2013vta , and anisotropic hydrodynamics Florkowski:2010cf ; Martinez:2010sc ; Ryblewski:2010ch ; Florkowski:2011jg ; Martinez:2012tu ; Ryblewski:2012rr ; Bazow:2013ifa ; Tinti:2013vba ; Nopoush:2014pfa ; Florkowski:2014bba ; Tinti:2015xwa ; Bazow:2015cha ; Bazow:2015zca ; Nopoush:2015yga ; Alqahtani:2015qja ; Molnar:2016vvu ; Molnar:2016gwq ; Bluhm:2015raa ; Bluhm:2015bzi ; Alqahtani:2017jwl ; Alqahtani:2017tnq ; Alqahtani:2017mhy ; Almaalol:2018ynx ; Almaalol:2018gjh . It was found that in all cases anisotropic hydrodynamics provided the best approximation to the exact attractor irregardless of the moment considered.

Importantly, it was shown that, when $m$ or $n$ are large, both the Navier-Stokes and second order viscous hydrodynamics results for the attractor failed to describe the exact solution. The fact that a subset of the exact moment solutions converge to something that is not well-described by traditional viscous hydrodynamics treatments means that we must refine our terminology a bit: instead of calling the convergence to the attractor “hydrodynamization”, we should instead call it pseudo-thermalization to emphasize that the attractor has a non-hydrodynamic nature reflected in the behavior of higher moments of the one-particle distribution function. In addition, we can associate the loss of memory of the precise initial conditions used with the pseudo-thermalization of the system. This is similar to the loss of memory which occurs if a system fully thermalizes, but with the universal state which emerges after pseudo-thermalization being far from equilibrium.

In this paper, we extend Ref. Strickland:2018ayk to study the effect of imposing number conservation on the dynamics and underlying non-equilibrium attractor. In RTA, one can enforce number conservation by introducing a fugacity (chemical potential) in both the dynamical and equilibrium distribution functions Florkowski:2012as ; Florkowski:2015cba ; Florkowski:2017ovw ; Almaalol:2018jmz . Requiring both energy and number conservation, one can derive two coupled integral equations which can be solved iteratively in order to obtain the effective temperature $T$ and fugacity $\Gamma$ as a function of proper-time. We demonstrate that for classical statistics, the integral equations can be written solely in terms of the rescaled variables $\hat{T}=T/T_{0}$ and $\hat{\Gamma}=\Gamma/\Gamma_{0}$ and the initial momentum space anisotropy $\xi_{0}$. As a result, one can construct the exact solution from any initial temperature and fugacity from a trivial scaling of the solution obtained for $\hat{T}$ and $\hat{\Gamma}$. We then determine the attractor solution to the coupled integral equations numerically by finding solutions which obey $\lim_{\tau\rightarrow 0}{\cal P}_{L}/{\cal P}_{T}\rightarrow 0$. We find that, in general, the resulting attractor solutions do not reach chemical equilibrium at late times, i.e. $\lim_{\tau\rightarrow\infty}\Gamma(\tau)\neq 1$. Despite the existence of a finite chemical potential at late times, we still observe attractor behavior in all moments and the full distribution function itself.

The structure of this paper is as follows. In Sec. 2, we briefly review how to rewrite the 0+1d RTA Boltzmann equation using boost invariant variables. In Sec. 3 we present the integral equations obeyed by the one-particle distribution function and all moments of the one-particle distribution function. In Sec. 4, we present our numerical results and discussion of the results. In Sec. 5, we present our conclusions and an outlook for the future.

2 Setup

Our starting point is the Boltzmann equation for massless particles

[TABLE]

in RTA,

[TABLE]

The relaxation time $\tau_{\rm eq}$ above can depend on proper time, however, since the system is conformal (massless) it must be proportional to the inverse effective temperature. The equilibrium distribution function $f_{\rm eq}$ may be taken to be a Bose-Einstein, Fermi-Dirac, or Boltzmann distribution. Here we will assume that $f_{\rm eq}$ is given by a Boltzmann distribution

[TABLE]

where $\Gamma(\tau)=\exp(-\mu_{\rm eff}(\tau)/T(\tau))$ is the effective fugacity with $\mu_{\rm eff}(\tau)$ and $T(\tau)$ being the local effective chemical potential and temperature, respectively.

The effective temperature $T$ and fugacity $\Gamma$ will be obtained via matching conditions which demand that the energy and number densities calculated from the dynamical distribution function $f$ be equal to the energy and number densities determined from the equilibrium distribution, $f_{\rm eq}$. The quantity $u^{\mu}$ represents the four-velocity of the local rest frame of the matter which herein we assume to be given by the transversally homogenous and boost invariant Bjoken flow (0+1d).

In equilibrium, for massless particles obeying classical statistics the particle density, entropy density, energy density, and pressure are

[TABLE]

2.1 Boost-invariant variables

For one-dimensional boost-invariant expansion, all scalar functions of time and space depend only on the longitudinal proper time $\tau\equiv\sqrt{t^{2}-z^{2}}$. In addition, the hydrodynamic flow $u^{\mu}$ has the following form $u^{\mu}=\left(\frac{t}{\tau},0,0,\frac{z}{\tau}\right)$ Bjorken:1982qr . One can introduce a space-like vector that is orthogonal in all frames and corresponds to the z-direction in the local rest frame of the matter $z^{\mu}=\left(\frac{z}{\tau},0,0,\frac{t}{\tau}\right)$. The requirement of boost invariance implies that $f(x,p)$ can depend only on three variables: $\tau$, $w$ and $\vec{p}_{T}$ Bialas:1984wv ; Bialas:1987en ; Florkowski:2013lza ; Florkowski:2013lya . The boost-invariant variable $w$ is defined by

[TABLE]

where $z$ is the spatial coordinate, not to be confused with the basis vector $z^{\mu}$. Using $w$ and $\vec{p}_{T}$ one can define

[TABLE]

From (5) and (6) one can easily find the energy and the longitudinal momentum of a particle

[TABLE]

The momentum integration measure in phase-space is

[TABLE]

In the following we shall consider massless partons, $m=0$.

2.2 Boost-invariant form of the kinetic equation

Making use of the boost-invariant variables introduced in the previous subsection, one finds $p^{\mu}\partial_{\mu}f=\frac{v}{\tau}\frac{\partial f}{\partial\tau}$, $p\cdot u=\frac{v}{\tau}$, and $p\cdot z=-\frac{w}{\tau}$. With this, Eq. (1) becomes simply Florkowski:2013lza ; Florkowski:2013lya

[TABLE]

with the finite chemical potential equilibrium distribution function (3) given by

[TABLE]

Note also that $f(\tau,w,\vec{p}_{T})$ is an even function of $w$ and depends only on the magnitude of the transverse momentum $\vec{p}_{T}$.

3 Exact solution for the distribution function

The formal solution of the kinetic equation (9) is Florkowski:2013lza ; Florkowski:2013lya

[TABLE]

where we have introduced the damping function

[TABLE]

For $\tau=\tau_{0}$ the distribution function $f$ reduces to the initial distribution function, $f_{0}$. For the conformal RTA solution we use the relation

[TABLE]

where ${\bar{\eta}}\equiv\eta/{\cal S}$ is the ratio of the shear viscosity $\eta$ to entropy density ${\cal S}$.

3.1 Initial distribution

Herein, for the initial condition we take the Romatschke-Strickland form Romatschke:2003ms with a classical Boltzmann distribution as the underlying isotropic distribution

[TABLE]

where, in the second equality, we have introduced the elliptical anisotropy parameter $\alpha\equiv 1/\sqrt{1+\xi}$ for convenience. The distribution function above reduces to an isotropic Boltzmann distribution if the anisotropy parameter $\alpha_{0}=\alpha(\tau_{0})=1$. If $\alpha_{0}=1$, the transverse momentum scale $\Lambda_{0}$ is equal to the system’s initial temperature $T_{0}$ and the initial microscopic fugacity $\gamma_{0}$ is equal to the initial effective fugacity $\Gamma_{0}$. In general, one must use Landau matching of the initial energy and number densities to fix $\Lambda_{0}$ and $\gamma_{0}$ in terms of $T_{0}$ and $\Gamma_{0}$. The resulting “matching conditions” are Almaalol:2018jmz

[TABLE]

where ${\cal H}(\alpha)\equiv{\cal H}^{20}(\alpha)$. The special functions ${\cal H}^{nm}$ needed are

[TABLE]

3.2 General moments of the distribution function

In order to solve Eq. (11), one can take a general moment of both sides using

[TABLE]

For $n=2$ and $m=0$, one obtains the energy density

[TABLE]

and, for $n=1$ and $m=0$, one obtains the number density

[TABLE]

Using the mass shell constraint, one can always rewrite the transverse momentum squared in terms of the energy and longitudinal momentum, so that any moment containing $p_{T}^{2\ell}$ can be written as a linear combination of the ${\cal M}^{nm}$ moments above. As a result, in the general case, we need to compute

[TABLE]

3.3 Integral equation obeyed by a general moment

Taking a general moment of Eq. (11) one obtains

[TABLE]

where

[TABLE]

One can rewrite the first term, which involves the initial values of the microscopic parameters $\gamma_{0}$ and $\Lambda_{0}$, in terms of the initial effective fugacity $\Gamma_{0}$ and temperature $T_{0}$ using Eqs. (15) and (16). Putting the pieces together, our final result for the general moment equation is

[TABLE]

Final equations

From Eq. (25) we can obtain two integral equations by evaluating the $n=1$ and $m=0$ and $n=2$ and $m=0$ moments which map to the number density and energy density, respectively, with the results being

[TABLE]

and

[TABLE]

where we used the matching conditions ${\cal E}={\cal E}_{\rm eq}(T,\Gamma)$ and $n=n_{\rm eq}(T,\Gamma)$ on the left-hand-side and the fact that ${\cal H}^{10}(\alpha)=2\alpha$ to simplify the second integral equation. Note, importantly, that one can divide the left- and right-hand sides of Eqs. (26) and (27) by the initial fugacity $\Gamma_{0}$ and rewrite them entirely in terms of $\hat{\Gamma}\equiv\Gamma/\Gamma_{0}$. As a result, one can solve these coupled integral equations with a given value of $\Gamma_{0}$, e.g. $\Gamma_{0}=1$, and then obtain solutions with different initial fugacity by scaling the result by the initial fugacity.111A similar scaling can be done with $T_{0}$, however, in addition to scaling Eqs. (26) and (27) by $T_{0}^{4}$ and $T_{0}^{3}$, respectively, one must also make a change of variables in the proper-time integrations by introducing $\hat{\tau}\equiv\tau T_{0}$.

4 Results

For our results, we solve Eqs. (26) and (27) numerically to obtain $T(\tau)$ and $\Gamma(\tau)$ given a set of initial values at $\tau_{0}$: $T_{0}$, $\Gamma_{0}$, and $\alpha_{0}$. For this purpose we wrote a CUDA-based GPU code which allows one to efficiently solve Eqs. (26) and (27) on very large proper-time lattices using a logarithmically-spaced grid. For all results presented herein we used a temporal lattice size of 4096 points and iterated the coupled integral equations until the effective temperature and fugacity converged to sixteen digits at all proper times. The code used to produce our results is included in the arXiv bundle for this paper and is also publicly available for download using the link provided in Ref. MikeCodeDB . Once the solutions for the effective temperature and fugacity are obtained, one can use Eq. (25) to obtain the proper-time dependence of any moment required. One can also use Eq. (11) to reconstruct the full one-particle distribution in a grid in momentum-space.

We will present the resulting exact solutions for the scaled moments

[TABLE]

where

[TABLE]

are the moments associated with an equilibrium Boltzmann distribution function. The scaled moments approach one at late times by construction and the rate at which they approach one provides a quantitative measure of how quickly the system thermalizes. Note that higher moments are sensitive to higher average momenta where the hydrodynamics assumption of small gradients could fail.

4.1 Attractor moments

In Fig. 1, we present sixteen panels containing our numerical results for the scaled moments ${\overline{\cal M}\hskip 1.42262pt}^{nm}(\tau)$ with $m,n\in\{0,1,2,3\}$. In all panels, the horizontal axis of Fig. 1 is the scaled proper-time $\overline{w}\equiv\tau/\tau_{\rm eq}=\tau T/5\bar{\eta}$, the black line is the exact solution for the attractor, and the various dashed/dotted curves are exact solutions with different values of $\alpha_{0}$. To obtain the attractor solution, we solved the coupled integral equations (26) and (27) with $\tau_{0}=10^{-3}$ fm/c, $T_{0}=1$ GeV, $\alpha_{0}=2.5\times 10^{-2}$, and $\Gamma_{0}=1$. To obtain the specific solutions (dashed/dotted lines), we solved the same coupled integral equations (26) and (27) with $\tau_{0}=10^{-1}$ fm/c, $T_{0}=1$ GeV, $\alpha_{0}\in 0.1\ldots 1.5$, and $\Gamma_{0}=1$. The range of initial anisotropies considered covers both very oblate and very prolate initial momentum-space anisotropy.

As can be seen from this figure, for generic initial conditions, all moments with $m>0$ visually converge to their respective attractors after a short rescaled time $\bar{w}\sim 2$. For ${\overline{\cal M}\hskip 1.42262pt}^{10}$ (scaled number density) and ${\overline{\cal M}\hskip 1.42262pt}^{20}$ (scaled energy density), we see that the constraints are properly enforced, resulting in these moments being equal to their equilibrium values at all proper times. For moments with $m=0$, we see a somewhat slower approach to the attractor. This is similar to what was found when not enforcing number conservation Strickland:2018ayk , however, herein we see smaller deviations from equilibrium. Despite these smaller deviations from equilibrium compared to the prior studies, moments with $m=0$ still converge more slowly than other moments. The slow convergence of moments with $m=0$ is related to the fact that they contain no powers of $p_{z}$ in their integrands and are, therefore, more sensitive to the free streaming term (first term) in Eq. (11). Free streaming results in momentum modes from the initial distribution being squeezed to smaller and smaller $|p_{z}|$ as a function of proper time (see Ref. Strickland:2018ayk for details).

Note, however, that although all scaled moments approach one in the large $\bar{w}$ limit, the system generically possesses a finite fugacity at late times.222One can adjust the late-time fugacity by changing the initial fugacity $\Gamma_{0}$. To demonstrate this, we plot the effective temperature, effective fugacity, and the pressure anisotropy associated with the attractor in Fig. 2. As can be seen from this figure, the scaled temperature evolution (left panel) obtained using attractor initial conditions shows characteristics of early time free streaming, for which the temperature scale (average momentum scale) is constant Mauricio:2007vz ; Martinez:2008di ; Martinez:2009mf ; Martinez:2010sc , followed by a power law decrease at late time indicative of hydrodynamic evolution. The attractor’s scaled effective fugacity (middle panel) decreases as a power law at early times and eventually saturates at late times. Finally, we see that attractor’s pressure anisotropy is large (highly oblate) at early times with $P_{L}\ll P_{T}$ and then slowly relaxes towards isotropy at late times.

4.2 Pseudo-thermalization time

In order to quantitatively assess the convergence of generic exact solutions to the attractor for each moment, one can compute the scaled time at which all solutions collapse to the attractor by requiring that $\max|{\overline{\cal M}\hskip 1.42262pt}^{nm}_{\rm i}(\overline{w}_{c})-{\overline{\cal M}\hskip 1.42262pt}^{nm}_{\rm attractor}(\overline{w}_{c})|<\delta_{c}$ for $i$ in the entire set of trial runs. In Ref. Strickland:2018ayk , $\delta_{c}=10^{-6}$ was chosen in order to require that the solutions were extremely well converged to the attractor. Herein, we will also consider the weaker convergence criteria of $\delta_{c}=10^{-2}$, which should correspond more closely to the time that one extracts when visually checking for convergence in Fig. 1.

In Figs. 3 and 4 we plot the convergence or pseudo-thermalization time $\bar{w}_{c}$ with $\delta_{c}=10^{-6}$ and $\delta_{c}=10^{-2}$, respectively. Figure 3 uses the strong convergence criteria of $\delta_{c}=10^{-6}$ which was the condition used in Ref. Strickland:2018ayk . The two top panels show $\bar{w}_{c}$ as a function of $m$ and $n$ and the bottom panel shows the case $m=0$ as a function of $n$. As can be seen from the top left panel, for $m\geq 2$, $\bar{w}_{c}$ is a decreasing function of $m$ and $n$. From the top right panel we see that for $m=1$ the pseudo-thermalization time increases at large $n$, but moments with $m>1$ have a pseudo-thermalization which decreases as $n$ increases. Turning to the bottom panel ($m=0$), we see that the $n=1$ and $n=2$ moments thermalize “instantly” since these are enforced by conservation laws333In this case we set $\bar{w}_{c}$ to be the smallest scaled time in the temporal grid. and we see a strong increase in $\bar{w}_{c}$ as $n$ increases. In the range of $n$ and $m$ shown, the maximum pseudo-thermalization time is $\bar{w}_{c}^{\rm max}\simeq 14$. This can be compared with the maximum pseudo-thermalization time obtained without enforcing number conservation (see Fig. 8 in Ref. Strickland:2018ayk ), which was $\bar{w}_{c}^{\rm max}\simeq 28$.

Turning to Fig. 4 we see the same three panels but now for the weaker convergence criterium of $\delta_{c}=0.01$. As can be seen from these figures, one obtains a shorter pseudo-thermalization time with the weaker convergence criteria, as expected. In addition, the scaled times extracted for moments with $m\neq 0$ are in the range of $\bar{w}_{c}\sim 1-3$. For a RHIC energy heavy-ion collision with a typical initial central temperature of 500 MeV at $\tau_{0}=0.1$ fm/c this translates into a physical pseudo-thermalization time of $\tau_{c}\sim 0.5-3$ fm/c with the precise value depending on the mode considered. We emphasize that the higher $n$ and $m$ moments converge more quickly and have pseudo-thermalization times on the low side of this window, while the lower $n$ and $m$ moments converge more slowly to the attractor.

4.3 Comparison with Navier-Stokes, vHydro, and aHydro

Finally, we compare the exact results for the scaled attractors moments with results obtained from anisotropic hydrodynamics (aHydro) and second-order viscous hydrodynamics (vHydro). For vHydro, we use the complete second-order viscous hydrodynamics equations of Denicol, Niemi, Molnar, and Rischke (DNMR) Denicol:2010xn ; Denicol:2011fa . For aHydro, we use the moments method introduced originally by Florkowski and Tinti Tinti:2013vba . For both vHydro and aHydro, the attractor is determined from the solution of a one-dimensional ordinary differential equation subject to the appropriate initial condition. For details concerning the determination of the attractor for both aHydro and vHydro, we refer to the reader to Ref. Strickland:2017kux .

For vHydro, one extracts $\bar{\pi}=\pi/\epsilon$ and, using this, one can reconstruct the solution for any moment using Strickland:2018ayk

[TABLE]

One can obtain the Navier-Stokes (NS) result by taking $\bar{\pi}=16\bar{\eta}/(9\tau T)$ in (30).

For aHydro, one extracts the anisotropy parameter $\alpha(\tau)$ associated with the attractor solution. Once this is determined one can use Eq. (4.4) of Ref. Strickland:2018ayk , modified to take into account finite fugacity, to obtain a compact expression for any moment

[TABLE]

The aHydro dynamical equations taking into account number conservation can be found in Ref. Almaalol:2018jmz . In Fig. 5, I compare the exact attractor (black solid lines) with the aHydro attractor (red dashed lines), DNMR attractor (blue long dashed lines), and the NS limit (green dot-dashed lines) for each moment. In all cases shown, aHydro provides a better approximation to the exact moments than the vHydro or NS frameworks. This is similar to what was found in the case where number conservation was not imposed Strickland:2018ayk . Finally, we note that, for aHydro, the moment with the best agreement is the $n=1$ and $m=1$ moment, for which the two results are virtually indistinguishable. This can be contrasted with Ref. Strickland:2018ayk which found that it was the $n=0$ and $m=1$ moment which was best described when number conservation was not imposed. It is not clear to us why this would be the case.

5 Conclusions

In this paper we extended previous studies of the conformal 0+1d kinetic non-equilibrium attractor in relaxation time approximation by imposing number conservation through the introduction of a dynamical fugacity (chemical potential). We derived two coupled integral equations for the effective temperature and fugacity which were then solved numerically to obtain the exact solution. We demonstrated that the resulting solutions exhibited convergence to a unique non-equilibrium attractor even though the system is out of chemical equilibrium generically ($\lim_{\tau\rightarrow\infty}\Gamma(\tau)\neq 1$). We found that, compared to the case where number conservation was not imposed, the moments converge to their respective attractors more quickly. Overall, however, we found that the behavior in the two cases is qualitatively very similar, providing further evidence that the non-equilibrium attractor is ubiquitous. We also compared the resulting attractor moments with predictions of different hydrodynamic frameworks. We found that anisotropic hydrodynamics provided the best approximation to the exact results for all moments.

Looking forward, herein we used a RTA collisional kernel and enforced number conservation by introducing a dynamical fugacity. It would be interesting to look at leading-order scalar field theory, in which case one only has 2 $\leftrightarrow$ 2 collisions and hence a theory which automatically conserves number. Such comparisons have been made in the context of aHydro in Ref. Almaalol:2018jmz where the authors studied both number-conserving RTA and scalar collisional kernels. Therein, it was shown that one could numerically extract the aHydro attractor for both RTA and scalar kernels, with the two being qualitatively similar. It would be very interesting to consider the 2 $\leftrightarrow$ 2 scalar kinetic theory using Monte-Carlo-based transport in order to compare with the exact results obtained herein using number-conserving RTA. It would also be interesting to make comparisons with the attractor extracted from the effective kinetic theory framework of Kurkela et al, particularly in the case that baryon number conservation is at play Kurkela:2015qoa ; Keegan:2015avk ; Kurkela:2018oqw ; Kurkela:2018xxd .

Acknowledgements.

We thank D. Almaalol for comments and suggestions. M.S. and U.T. were supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award No. DE-SC0013470.