Analytical attractor for Bjorken expansion

We present an analytic solution of a simple set of equations that govern the expansion of boost-invariant plasmas of massless particles. These equations describe, approximately, the early time, collisionless regime, and the transition to hydrodynamics at late time. Their mathematical structure encompasses all versions of second order hydrodynamics when applied to Bjorken flows. The analytic solution provides an explicit expression for the attractor solution that connects the collisionless regime to hydrodynamics. It also constitutes a neat example of an application of the theory of resurgence in asymptotic series.

The success of the fluid dynamical description of ultrarelativistic heavy ion collisions has triggered recently a lot of theoretical work [1,2]. While hydrodynamics is often viewed as an effective theory for long wavelength modes with microscopic degrees of freedom near local equilibrium, implying small gradients and small mean free paths, a number of recent studies suggest that it may work even when these conditions are not fully satisfied. It was found for instance that viscous corrections can account for the large pressure anisotropy in the longitudinally expanding plasma formed in heavy ion collisions, for either a strongly coupled system in holography [3], or a weakly coupled system in kinetic theory [4]. Thus the general question of how hydrodynamical behavior emerges in various systems has become a prominent one. In this context, an important progress has been the identification of attractor solutions in dynamical equations whose late time behavior is hydrodynamical [5][6][7][8]. Besides its physical relevance, this issue also connects with beautiful mathematical developments in the theory of resurgence in asymptotic series [9,10].
This letter aims to contribute to this general discussion by providing an explicit analytical solution for a simple set of equations that describes the transition from kinetics to hydrodynamics for a rapidly expanding plasma of massless particles. The system that we consider is an idealization of the early stage of a high-energy heavy-ion collision, where matter expands longitudinally along the collision axis in a boost invariant fashion, the so-called Bjorken expansion [11]. Our starting point is a simple kinetic equation, p µ ∂ µ f (t, x, p) = C[f ], where f denotes a distribution function for massless particles, and C[f ] is a collision term treated in the relaxation time approximation [12], a widely used approximation [2].
It has been shown in previous works [13] that the solu- * jean-paul.blaizot@ipht.fr † cliyan@fudan.edu.cn tion of the kinetic equation can be accurately reproduced by a set of coupled equations for the two independent components of the energy-momentum tensor, the energy density ε = L 0 and the viscous pressure, expressed as the difference between the longitudinal and transverse pressures P L − P T = L 1 . Owing to the symmetries of the Bjorken expansion, L 0 and L 1 depends only on the proper time τ , and obey the coupled equations where τ is the proper time. The coefficients, a 0 = 4/3, a 1 = 38/21, b 1 = 8/15, and c 0 = 2/3, are pure numbers whose values are fixed by the geometry of the expansion. The last term, proportional to the collision rate 1/τ R , isolates in a transparent way the effect of the collisions. We shall allow τ R to depend on τ (see below). Eqs. (1) are the first in an infinite hierarchy of equations for moments defined by L n ≡ d 3 p pP 2n (p z /p)f p (t, x, p), with P n (x) a Legendre polynomial and p = |p| [13] 1 . The goal of this paper is to present an exact analytical solution of Eqs. (1). Without the last term in Eq. (1b), Eqs. (1) describe (approximately) free streaming: the moments evolve as power laws governed by the eigenvalues of the linear system. The collision term in Eq. (1b) produces a damping of L 1 and drives the system towards isotropy, a prerequisite for local equilibrium. When L 1 = 0, the system behaves as in ideal hydrodynamics L 0 ∼ τ −a0 .
However there is more in the approach to hydrodynamics than the simple vanishing of L 1 . Indeed, the term −L 1 /τ R in Eq. (1b) suggests the possible presence of exponential terms in the solution, in addition to the power laws originating from the expansion. Such exponential contributions would spoil the hydrodynamic gradient expansion, which, in the present context, is an expansion in inverse powers of τ . However, these contributions cancel thanks to the relation L 1 /τ R −b 1 L 0 /τ valid at late time [13]. This relation, in fact, happens to be the leading order constitutive equation relating the viscous pressure to the viscosity, with η = (b 1 /2)ετ R being the standard value of the viscosity. We shall return to these exponential contributions later in this letter.
The connection between Eqs. (1) and hydrodynamics is tight. Eq. (1a) translates the conservation of the energy momentum tensor, ∂ µ T µν = 0, for Bjorken flow. An equation similar to Eq. (1b) was introduced by Israel and Stewart [17] to overcome limitations of the Navier Stokes equation in the relativistic context. In fact, as discussed in [18,19], if one adjusts properly the coefficients a 1 , b 1 and the relaxation time τ R (a 0 and c 0 are fixed by energy-momentum conservation, as just mentioned), all the known formulations of second order hydrodynamics (e.g. [20,21]) can be mapped, in the context of Bjorken flows, into Eqs. (1), and share therefore the same mathematical structure. Indeed, the authors of Ref. [7] (see also [22]), using some version of Israel-Stewart hydrodynamics, obtained a solution that is a particular case of that presented here.
The use of a time dependent relaxation time τ R (τ ) allows us to capture the qualitative features of more realistic calculations. The solution presented in this letter holds for a time dependence of the form τ R ∼ τ 1−∆ , with ∆ constant 2 [23]. Commonly used are a constant τ R (∆ = 1), or τ R ∼ T −1 with T the effective temperature (∆ ≈ 2/3). Note that as long as ∆ > 0, the expansion rate decreases faster than the collision rate, and the system is driven to hydrodynamics at late time. When ∆ < 0, the expansion eventually overcomes the collisions and the system evolves towards the free streaming regime. The limiting case ∆ = 0 mimics the system of hard sphere scatterings recently studied in [24], and is also close to more realistic QCD kinetics with 2-to-2 scatterings [25]. In this case, collisional effects perfectly balance those of the expansion, resulting in a stationary state that differs from hydrodynamics. More sophisticated time dependence have also been considered [26]. A detailed description of these physical situations and the corresponding solutions as a function of ∆ will be presented in a forthcoming publication [19]. In this letter we focus on the solution for ∆ > 0.
To proceed, it is convenient to measure the time in units of the instantaneous relaxation time and define the dimensionless variable w ≡ τ /τ R (τ ) ∼ τ ∆ , which plays here the role of an inverse Knudsen number. We also define This quantity g(w) may be viewed as the exponent of the power laws obeyed by the energy density at early or late times. It is also a measure of the pressure asymmetry. In particular, the second relation, which follows easily from Eq. (1), shows that in the free streaming regime where P L = 0, g = −1, while in the hydrodynamical regime where P L = ε/3, g = −4/3. In terms of g(w), Eqs. (1) become a first order nonlinear ODE, (3) In the absence of collisions, this non linear equation has two fixed points, that we refer to as unstable (g − ) and stable (g + ) free streaming fixed points, whose values coincide with the eigenvalues of the linear system (1). Numerically, g + = −0.929, g − = −2.213 3 . This fixed point structure continues to play a role when collisions are switched on [18]: The unstable fixed point moves to large negative values, while the stable fixed point g + evolves adiabatically 4 [27] to the hydrodynamic fixed point, g * = −4/3. The location of this "pseudo fixed point" as w runs from 0 to ∞ corresponds (approximately) to what has been dubbed "attractor" [5].
Here the attractor 5 is understood as the particular solution of Eq. (3), g att (w), that connects g + as w → 0 to g * as w → ∞. As was shown in [18] the fixed point structure at small w is only moderately affected by the higher moments (L n , with n ≥ 2), which explains why the two-moment truncation (1) is a good approximation to the solution of the full kinetic equation for the lowest two moments 6 . Finally, in view of the close relation between the coupled equations (1) and hydrodynamics, the same fixed point structure plays a prominent role in second order hydrodynamics.
By setting g(w) + a 0 = ∆(b − a) − w + wy (w)/y(w), with b and a constants, andw = w/∆, one transforms Eq. (3) into a linear, second order ODE: The solutions y(w) of this equation (Kummer's equation) are confluent hypergeometric functions M (a, b, w/∆) and U (a, b, w/∆) [29]. The parameters a and b can take one of two possible sets of values, a ± = 1 − g∓+a0 ∆ , and b ± = 1 ± g+−g− ∆ . However, it can be shown that the physically meaningful solution corresponds to the choice a + , b + 7 . Thus, from now on we shall simply set a + = a, b + = b and write the general solution for g(w) in the form where A is a constant to be fixed by the initial condition 8 . Eq. (5) is the exact solution to Eq. (3) for all ∆ > 0. Since the structure of the solution does not change much as ∆ varies, we set from now on ∆ = 1. 9 An illustration of the solution for various initial conditions is provided in Fig. 1 for five values of A. As can be seen, the various solutions eventually merge onto the solid line, the attractor solution which connects g + to g * . This attractor solution is obtained for A = 0, viz.
5 For a definition of the notion of attractor more in line with that in use in the literature on dynamical systems we refer the reader to [16,28] and references therein. 6 The main effect of the higher moment is to bring g + to its correct value, thereby eliminating the regions of the attractor corresponding to negative longitudinal pressure. 7 One can show for instance that only this choice leads to an increase of entropy in the late time hydrodynamic regime [19]. 8 Note that a single initial condition is required for the three versions of the equations of motion, Eqs. (1,3,4). 9 Solutions for values of ∆ > 0 other than 1 can in fact be deduced from the solution for ∆ = 1 via a suitable rescaling of the parameters of the equation. See [19] for details.
This attractor solution possesses a well behaved expansion at small w. Physically, such an expansion represents perturbative corrections to free-streaming due to collisions, and it emerges from Eqs. (1) as one expands in powers of the last term in Eq. (1b). When the parameters a and b differ from non-negative integers, which is the case here, the function M (a, b, w) is analytic and expressible as a convergent series in w. However, in contrast to the function M itself, the series g att (w) = n=0 γ att n w n , with γ att 0 = g + , has a finite radius of convergence determined by the location of the zero of M (a, b, w) in the denominator of Eq. (6) that is closest to the origin. The Borel sum of the series, after analytic continuation, converges to the exact function, as it has been verified numerically in [8].
When A = 0, the solution (5) involves the function U (a, b, w) which contains a singular contribution ∼ w −b at the origin [29]. The small w expansion becomes then more intricate and takes the form of a trans-series [16] where the expansion coefficient γ (m) n for m > 0 depends on A. We note that the ratio γ This implies that all solutions (but the attractor) start from the unstable freestreaming fixed point at w → 0 + , even though at any small but finite w one can find a value of A such that a given solution can take any desired value.
At large w, a solution g(w) of Eq. (3) can be obtained as an expansion in powers of 1/w, giving rise to the hydrodynamic gradient expansion, In early works, the coefficients f n were determined iteratively, their observed factorial growth reflecting the asymptotic character of the expansion [30]. The knowledge of these coefficients could then be exploited in a Borel summation. The Borel sum is the inverse Laplace transform,g where B[g(t)] is the Borel transform of the original series, viz.
and C is a contour in the complex plane that joins 0 and ∞. The singularities of the Borel transform, whose knowledge is required to calculate the integral in (9), are typically determined by obtaining numerically a large number of terms in Eq. (10), and using a Padé approximant extrapolation [5,31]. Using such a procedure g(w) was obtained in the form of a trans-series [5] where σ is a complex expansion parameter and is itself an asymptotic series. Note that g (0) (w) = g hydro (w). The quantity ζ(w), which captures exponential corrections, can be determined by perturbing the solution around the hydrodynamic gradient expansion [10]. When applied to the solution of Eq. (3), this yields [18] The same quantity ζ(w) plays a role in the Borel summation. Indeed the Borel transform in Eq. (10) exhibits a branch cut on the positive real axis, which forces the integration contour C in Eq. (9) to be deformed into the complex plane. This generates an imaginary part in the Borel sum, which can be shown to be proportional to ζ(w) [10]. For the Borel sum to be real, this imaginary part has to cancel against the next order in the transseries in Eq. (11), which involves the imaginary part of the complex expansion parameter σ. Such cancellations between coefficients of different members of the transseries is an illustration of the well-known resurgence phenomenon (see. e.g. [9,10]), to which we shall return.
Since we have the analytical solution (5) at hand, together with the known asymptotic properties of the confluent hypergeometric functions [29], one can verify very explicitly the features of the solution that we have just recalled, without going through a numerical Borel summation (while at the same time allowing for a precise test of this procedure -see Fig 2). The asymptotic forms of the confluent hypergeometric functions M (a, b, w) and U (a, b, w) contain both an exponential factor e w and a non-analytic power factor w α , with α determined by the parameters a and b [29]. In the analytical solution Eq. (5) these factors naturally combine into ζ(w), yielding the asymptotic expression where and F(a, b, z) is the asymptotic series This series is Borel summable. The Borel transform is a Gauss hypergeometric function, a, b, 1, z), which has a branch cut running from z = 1 to ∞. The Borel sumF is the inverse Laplace transform of 2 F 1 (a, b, 1, z). Its analytical expression reads [32] Clearly, the expansion of Eq. (14) in powers of σ generates a trans-series of the form of Eq. (11). The leading order term (O(σ 0 )) is the hydrodynamic gradient expansion, whose Taylor expansion coefficients can be verifed to coincide with those obtained by the method indicated after Eq. (8). Moreover, the Borel summability of F, and the availability of the analytical expression of the Borel sum F, suggests the following analytical expression for the Borel sum of the hydrodynamic gradient expansioñ obtained by substituting F byF in Eq. (18), withF given in Eq. (17). In fact this is an exact procedure. Indeed one can verify that the substitution of F byF in Eq. (14) reconstructs the exact solution. The higher order asymptotic series g (m) in the trans-series (11) can be determined from Eq. (14) in terms of the function F(a, b, z) as well, and their Borel sumsg (m) obtained similarly by substituting F byF in the corresponding expressions. For instance, the analytical Borel sum of g (1) (w) reads The Borel sum of the hydrodynamic expansion has an imaginary part, related, as mentioned above, to the branch cut of the Borel transform. This imaginary part is cancelled by the real part of the first correction, proportional to the imaginary part of σ. More precisely, we have where S 1 , a so-called Stokes constant, is proportional to the imaginary part of σ (and hence independent of the initial condition). This cancellation of imaginary parts is one manifestation of the resurgence phenomenon [33].
Resurgence also implies relations between the coefficients of the trans-series, which can be easily verified. For instance, the n th order term in the hydrodynamic gradient expansion is related to higher order terms in the transseries by (see e.g. Eq. (39) in Ref. [10]). Already the first two terms in this expansion yield an excellent approximation to f (0) n for not too small n. But perhaps the most powerful test of the resurgence properties is to consider the complete trans-series solution [33,34], where σ R = Re σ is the real part of the σ. This equation is non-trivial. It encompasses all cancellations of imaginary parts, according to the resurgence relations. Fig. 2 displays the trans-series solution for the attractor. The convergence towards the exact solution is illustrated by adding to the Borel sum of the hydrodynamic gradient expansiong (0) , the contributions proportional tog (1) (NLO), and tog (2) (NNLO). Including all orders in the trans-series amounts to substitute F byF in Eq (14), which, as we have mentioned, reproduces the exact solution (see Fig. 2).
Further examples of mathematical relations that can be explored with the analytic solution will be discussed in [19]. Let us just mention here the rearrangement of the trans-series, leading to the so-called "transasymptotic matching" (see e.g. [10,28] and references therein). In essence, the procedure consists in summing the exponential corrections in front of each power of w −n , thereby increasing the accuracy of the trans-series in the region of not so large w where these corrections become of the same order of magnitude as the corresponding terms in the gradient expansion. One then finds that the resulting trans-series has here a simple structure. Indeed, the coefficient of each order n is a polynomial of order n + 1 in ζ, as can be easily deduced from Eq. (14). For instance at order 1/w we get We recognize in the first term the leading order in the gradient expansion, with the coefficient b 1 proportional to the viscosity (see the discussion following Eqs. (1)). Following [35], one may then interpret the whole collection of terms within the brackets in Eq. (25) as an effective viscosity. This effective viscosity, normalized to the standard value, decreases as w decreases. This is the same behavior as that observed in [13], but for a completely different reason: there the renormalization of the viscosity originates from the contribution of the higher L n moments that are neglected in the two-moment truncation. This observation points to the fact that various factors, of different dynamical origins, may contribute to modify the value of the effective viscosity extracted in heavy ion collisions. As a final remark, we note that the attractor solution (6) receives exponential corrections from all the "nonhydrodynamic" sectors of the trans-series (11), an observation also made in [5]. From that perspective, we find the often used terminology "hydrodynamic attractor" somewhat misleading. As we have emphasized in this letter, we prefer to view the attractor, in the present context, as the particular solution that connects two fixed points of the non linear equation (3) for the pressure asymmetry. The two fixed points are associated with different physics: one corresponds to hydrodynamics, the other to the collisionless regime. The requirement of connecting these two different regimes is reflected in the intricate mathematical structure of the solution. As we have emphasized, the particular fixed point structure which emerges naturally from kinetic theory turns out to be shared by all versions of second order hydrodynamics applied to Bjorken flows.