Non-perturbative rheological behavior of a far-from-equilibrium expanding plasma

For the Bjorken flow we investigate the hydrodynamization of different modes of the one-particle distribution function by analyzing its relativistic kinetic equations. We calculate the constitutive relations of each mode written as a multi-parameter trans-series encoding the non-perturbative dissipative contributions quantified by the Knudsen $Kn$ and inverse Reynolds $Re^{-1}$ numbers. At any given order in the asymptotic expansion of each mode, the transport coefficients get effectively renormalized by summing over all non-perturbative sectors appearing in the trans-series. This gives an effective description of the transport coefficients that provides a new renormalization scheme with an associated renormalization group equation, going beyond the realms of linear response theory. As a result, the renormalized transport coefficients feature a transition to their equilibrium fixed point, which is a neat diagnostics of transient non-Newtonian behavior. As a proof of principle, we verify the predictions of the effective theory with the numerical solutions of their corresponding evolution equations. Our studies strongly suggest that the phenomenological success of fluid dynamics far from local thermal equilibrium is due to the transient rheological behavior of the fluid

These findings hint at the existence of a new theory for farfrom-equilibrium fluid dynamics [19][20][21]. Although little is known of this theory, it has been identified that the onset of hydrodynamics depends on the decay of new degrees of freedom, the so-called non-hydrodynamic modes [22][23][24][25][26], whose nature is purely non-perturbative [18,19,22,23,25,[27][28][29][30][31][32][33][34][35]. To this end, it is necessary to study how the nonlinear microscopic dynamics between different modes influences its macroscopic response. But a full-fledged study of nonequilibrium physics is not practically possible unless an effective truncation scheme is instead proposed. Thus, the challenge is to identify a tractable truncation scheme which captures the relevant non-equilibrium dynamics while being simple enough to perform concrete calculations. This task requires creating new computational tools and techniques, which is not otherwise achievable via the usual linear response theory.
In this Letter, we develop a new approach to studying hydrodynamization and nonlinear transport processes, which goes beyond linear response regime. Our treatment is highly nonlinear in the Knudsen and inverse Reynolds numbers, and thus, appropriate for the study of far-from-equilibrium dynamics and understanding the non-perturbative information of short-lived non-hydrodynamic modes. With the help of transasymptotic matching [36], we calculate the sum of short-lived transient non-hydrodynamic modes at a fixed order in the gradient expansion, and the transport coefficients thereby become dynamical quantities evolving with time toward their respective asymptotic equilibrium values obtained by linear response theory. Knowing the late-time behavior of the fluid, we backtrace the flow of transport coefficients all the way into the UV by including as many non-perturbative contributions as possible, to the transient non-hydrodynamic modes, which naturally build up the formal exponential solutions to the fluid evolution equations. This consequently provides renormalization group (RG) equations governing the RG flows of the dynamical transport coefficients, which by construction, converge to the correct values as the fluid reaches the thermal equilibrium. This simple situation is analogous to the case of non-Newtonian fluids [37], where transport coefficients such as shear viscosity, become nonlinear functions of the gradient velocity and thus, their evolution is dictated by the macroscopic deformation history of the fluid [38].
We consider an expanding conformal plasma undergoing Bjorken flow. The microscopic description of this system is governed by the relativistic Boltzmann equation within the relaxation time approximation (RTA-BE). The problem of solving the RTA-BE is transformed into determining the evolution of the moments through their nonlinear kinetic equations. The solutions of the moments can be written in terms of a multi-parameter trans-series [36], which correctly yields the information of the attracting section around the thermal equilibrium fixed point for the gas at late time, otherwise known as the "attractor". Being able to go beyond the attracting region distinguishes the predictions of linear response theory and the physics of non-equilibrated fluid dynamics.
The nonlinear relaxation process described by the Boltzmann equation can be understood in terms of the evolution of the moments of the distribution function [42]. Thus, we use the following ansatz for f k k k [43] f p p p = f eq.
where L (3) n and P 2l denote the generalized Laguerre and Legendre polynomials, respectively. From this expression the moments c nl are directly read as where A nl = 2π 2 (4l + 1) n! (n+3)! and · · · ≡ p · · · f p p p . The hydrodynamic equilibrium is given by c eq. nl = δ n0 δ l0 . Here, we only study the low energy modes, i.e. n = 0, and thus we denote c 0l by c l . Furthermore, we assume that the initial values of c nl are obtained from the RS distribution function ansatz where ξ 0 measures the momentum anisotropy and R(ξ 0 ) is given by Eq. (7a) in [45].
For the ansatz (2) the energy-momentum tensor is T µν = p µ p ν = ε u µ u ν + Ξ µν P T + l µ l ν P L where the energy density ε, the transverse (P T ) and longitudinal (P L ) pressures are given as follows: The Landau matching condition implies c 0 ≡ 1. The only independent component of the normalized effective shear tensor isπ = π ς ς /ε = 2/15 c 1 . The formal trans-series solution. In a scale invariant system governed by (1), an asymptotic series solution can be expanded in terms of the dimensionless parameter w −1 := (τT (τ)) −1 [22]. Then, the underlying non-autonomous dynamical system consisting of evolution equations of the moments and temperature has its dimensions d = L + 1 (with l = 1, . . . , L) reduced by one. Furthermore, w can be interpreted as an 'energy' scale for the renormalization scheme arising naturally from summing all the non-perturbative sectors around the (perturbative) asymptotic expansion of the modes c l . Hence, preparing a dynamical system for the modes c l (w) automatically amounts to having a renormalization group equation on the phase space of the Bjorken flow built out of these modes. These flows emanate from the vicinity of a UV point discussed below and converge to c eq l , alternatively referred to as the IR fixed point due to being reached when w → ∞ (or, equivalently, Kn , Re −1 → 0).
The evolution equations of the moments c l obtained from the conservation laws and the RTA-BE (1) (see Supplementary Material) can be written as 1lc l and c l are the components of c. This system is L-dimensional, rank 1, and level-1 vector differential equation in a neighborhood of an irregular singularity, w = ∞. In the analysis of asymptotic solutions of (5), two terms play a crucial role: (i) the coefficient of order O(w 0c ) term, i.e. the diagonal matrix of Lyapunov exponentsΛ determining the rates at which exponential stability is reached at equilibrium in L directions of the flow space, and (ii) the coefficient of order O(w −1c ), i.e. the (diagonalized) mode-tomode coupling matrixB D that shows how each mode c l is coupled to the neighboring modes. Knowing these two elements paves the way to recasting the dynamical system into the form (5) specifically prepared for asymptotic analysis shown first in generic nonlinear rank-1 systems of ODEs in Ref. [46]. Sincê B is not diagonal as it stands, the prepared form of the dynamical system was written in the basis ofc that diagonalizes it. Notice that in the RTA there is only one parameter θ 0 and thus, the Lyapunov exponents are all identically maximal and given byΛ = 3 2θ 01 .
In building the matrix U it is appropriate to choose the eigenvectors ofB in such a way that U −1 Ll = 1 for any l. Doing is the floor function and b L ∈ R for odd L. Also important is to state that rescaling the eigenvectors immediately corresponds to changing the expansion parameters in the trans-series ansatz defined below and for the sake of consistency, we will stick to the above condition for the eigenvectors from now on. Finally, let us write down the trans-series ansatz forc l : Here, ζ ζ ζ n (w) = exp(−(n · S)w)w n·b = ∏ L l=1 (ζ l (w)) n l where ζ l (w) = exp(−S l w)wb l , σ σ σ n ∈ C L are the product of expansion parameters ∏ L l=1 σ n l l , and n 1 , . . . , n L are non-negative integers. The dot indicates the inner product, A · B = ∑ L l=1 A l B l . It is better to exploit vector expressions for the non-perturbative sector number aka trans-monomial order n = (n 1 , . . . , n L ), the anomalous dimension of the pseudo-modesc l ,b = (b 1 , . . . ,b L ), and the Lyapunov exponents S = (S 1 , . . . , S L ).
By substituting the ansatz (6) in the differential equation (5), one can obtain a recursive relations for the coefficients as where δ n,0 = δ n 1 ,0 · · · δ n L ,0 . With a little bit of algebra, Eq. l,0 = 0 for any l. So physically, c 1 is the slowest mode that contributes to the observables in (4). The zeroth order term in the first non-perturbative sector (|n| = 1) is always a constant normalizing the trans-series expansion parameters σ l . For convenience, we choose to setũ (n) l,0 = 1 for n l = δ l ,l . Other trans-series coefficients can be recursively determined order by order in both k and n in (7).
In addition, at |n| = 1 order Eq. (7) gives S l = 3 2θ 0 and b l = −b l + 1 5 . For example, at first two truncation orders L = 1, 2 the anomalous dimensions are given asb 1 , respectively. It is notable that, physically, all c l must be real-valued butc l are in general complex, justifying the name 'pseudo-modes'. The reality condition for c l is related to the choice of the initial values (expansion parameters) σ l , which imposes the constraints σ 2l−1 = σ * 2l ∈ C for l = 1, . . . , [L/2], and σ L ∈ R if L is odd. Every flow line that evolves to the attractor of the equilibrium fixed point as w → ∞, can be constructed from the set of asymptotic solutions to the dynamical system (5) for some σ l . The allowed set of σ l shapes the basin of attraction at any L. Initial data and UV information. Determining the l-th expansion parameter σ l is the most challenging task related to Eq. (6). This partly goes back to the fact that there are L equations involved in the system of ODEs of (5), whereas the exact numerical expression of (3) is always prepared by only three initial values: τ 0 , T 0 and the parameter ξ 0 . Among these constants, the search for ξ 0 -dependence of σ l is more difficult as the other two explicitly enter the trans-series viaÉcalle time w itself. In order to obtain σ l , we use the least squares optimization to fit the trans-series of all the modes involved at order L simultaneously on the Stokes line. One, however, needs to be aware of "Stokes phenomenon", and the UV data that eventually dictate the fate of flow lines.
First, we have to note that capturing the divergence of original asymptotic expansion (n = 0) in (6) is done by initially Borel transforming the series. This yields some sort of singularity along the physically relevant Stokes line, that is, ℑ(n · Sw) = 0. Then, the Borel resummation of the original asymptotic series picks up a residue or discontinuity around the Stokes line (given more explicitly byÉcalle's bridge equation [47]). Imposing the reality condition, one concludes that the overall constant appearing in the imaginary piece, aka "Stokes constant", has to cancel the imaginary part of the σ l . For L = 1, it is shown that this constant is ℑ(σ 1 ) = −πC 0 θ −b 1 0 with the overall fitting factor of C 0 ≈ 0.4898.
Second, the flow lines heading to the IR should always take values in the basin of attraction, on the boundary of which lies the UV points [18]. These are actually the RG flows that unveil the general behavior of the modes c l in kinetic theory and consequent renormalization of transport coefficients far from equilibrium. In the context of RTA Boltzmann theory, the RG time should appropriately be chosen as t = log(ρ) with ρ = τT 0 /θ 0 , and τ 0 T 0 /θ 0 is a UV cutoff. Then, for L = 1, we get the following autonomous system of equations describing the RG flows in the 3d macroscopic Bjorken phase space: Unless τ 0 = T 0 = 0 (while keeping θ 0 finite), the theory does not have a UV completion. At this limit, there are two UV fixed points: (ρ, T, c 1 ) UV = (0, 0, −3.03), an index-2 (weakly unstable) saddle point, and (0, 0, 6.60) which is an index-1 (strongly unstable) saddle point shown in Fig. 1. The index represents the number of attracting directions that lie on the T -c 1 plane while the orthogonal direction ρ is obviously re- pelling. Considering the full theory L = ∞, the UV points move into the physical domain and become (ρ, T, c) UV = 0, 0, (−2.5, (−1) l (4l + 1) 2l l 4 −l ∞ l=2 ) and 0, 0, (5, {4l + 1} ∞ l=2 ) , of which the former gives the "free-streaming" limit (ξ 0 → ∞). The latter explains a recently found solution in Ref. [48]. At this point, the longitudinal pressure is maximum P L = ε, which in turn causes a huge amount of pressure along the z-direction due to ξ 0 → −1, while there is a slowdown in the transverse direction because P T = 0. For a flow emanating from the neighborhood of this point, there is a local relaxation happening at early times, an indication of more attracting directions in the UV, as is obvious from the larger index of the saddle point A (free-streaming) compared to that of the saddle point B in Fig. 1. Thus, we suitably define Ξ 0 := 1/(1 + ξ 0 ) that measures the proximity of the initial point of flow lines to either UV fixed points. The smaller (or larger) Ξ 0 , the closer flow lines will be to the saddle points A (or B). Non-perturbative contributions renormalize the transport coefficients. To find the renormalized coefficients, we straightforwardly do the sum over n in (7) that leads to the master recursive relatioñ whereC l,k (σ σ σζ ζ ζ(w)) = ∑ ∞ n≥0 σ σ σ n ζ ζ ζ n (w)ũ (n) l,k ,ζ l = ∂ log ζ l . Notice that (9) is composed of L nonlinear PDEs for every k, and it is in general difficult to solve it analytically for L > 1 [49]. Nonetheless, it is exactly solvable for L = 1 at any order k. The exact results for k = 0, 1 are given byC 1,0 (σ 1 ζ 1 ) = −20W ζ andC 1,1 (σ 1 ζ 1 ) = − 8 15 θ 0 (50W 3 ζ + 105W 2 ζ + 36W ζ + 5)/(W ζ + 1) respectively, where W ζ = W (−σ 1 ζ 1 /20) is the Lambert-W function. Fig. 2 compares the bare and renormalized c 1 for different initial conditions for L = 1, and the agreement between the exact numerical and trans-series solutions is overall remarkable. We may directly extract the renormalized transport coefficients from the near-IR asymptotic expansion of nonhydrodynamic modes c l found in [50] all the way into near-UV region by including the exponentially small contributions coming from the nonlinearities introduced in Eq. (5) by modeto-mode coupling terms. Here, we have to pay attention to two important characteristics of an RG flow in connection with this simple picture: (1) The cutoff τ 0 T 0 /θ 0 for a fixed finite θ 0 , is merely a physical bound on the near-UV theory such that, as in Fig. 1, the flows are initiated away from the T -c 1 plane and well inside the basin of attraction of the IR fixed point.
(2) The anomalous dimensionsb l and Lyapunov exponents S l are both the IR data appearing in the trans-series of each non-hydrodynamic mode. Along an RG flow initiated at some near-UV point when log (w 0 /θ 0 ) > −∞ or log ρ > −∞ in the basin of attraction, these two parameters never flow with the RG time. In other words, these data are invariant under an RG transformation of energy-momentum tensor governing the dynamical system (5), which pushes forward the flow lines starting near a fixed point and flowing to another fixed point.
It is, therefore, sought to find a non-perturbative approach to determine the far-from-equilibrium transport coefficients from a rheologic interpretation of the fluid. The deformation history of the fluid is traced in the transport coefficients by considering their nonlinear functional dependence on the velocity gradient tensor, e.g. η ≡ η(σ µν ), which is often studied at a perturbative level in non-Newtonian fluid dynamics. This has given rise to some interesting physical phenomena known as shear thinning and thickening [37].
The second order gradient expansion of the distribution function [50][51][52][53] determines the asymptotic form of c 1 to be Thus, by comparing this result with the one obtained from the trans-series of c 1 (for L = 1), that is c 1 = ∑ ∞ k=0 F 1,k (w)w −k , we see that the first non-perturbatively corrected asymptotic term of order (9). In a similar fashion, one can obtain T (λ 1 /s − τ π η/s) = 9 80 F 1,2 (w) from c 1 . The RG flow equation of these transport coefficients can then be derived from the RG flow of moments. Considering a generic observable O(w) constituting the resummed coefficients F 1 = {F 1,k |k ∈ Z ≥0 }, we can write down its (L = 1, θ 0 < ∞) RG flow equation explicitly as where ζ 1 (w) = e − 3w 2θ 0 wb 1 , F 1,k (w) ≡ F 1,k (σ 1 ζ 1 (w)) and ∆ k (w) = ∂O(F 1 (w))/∂F 1,k (w). As can be seen, the transport coefficients get a fully non-perturbative expression through ζ 1 . By treating ζ 1 as an effectively independent variable from w [36], F 1,k (ζ 1 ) can be regarded as the coefficients of a power series in w −1 , i.e. c 1 (ζ 1 , w) = ∑ ∞ k=0 F 1,k (ζ 1 )w −k . Hence, the pieceζ 1 F 1,k in (11) is calculated by using the Cauchy integral formula to pick out the kth coefficient ofζ 1 c 1 : where the beta function β 1 is defined through (5) as dc 1 d log w = β 1 (c 1 , w). Setting O(w) = − 3 40 F 1,1 and O(w) = 9 80 F 1,2 renders the RG flow equation in terms of η/s and T (λ 1 /s − τ π η/s), respectively. For different initial values of σ 1 , Fig. 3 shows the w-dependence of η/s and T (λ 1 /s − τ π (η/s)). As w → ∞, both observables asymptotically approach to their equilibrium values at the IR fixed point. There is, however, a shear thinning effect in the UV regime due to the early-time free streaming expansion, as seen in Fig. 3, until the interactions kick in and the flow initiates its relaxation onward exponentially in the shear thickening phase. We also observe that in accordance with the second law of thermodynamics, the renormalized η/s remains always positive, and its shape is identical (likewise for T (λ 1 /s − τ π (η/s))) regardless of the initial values. There is an overall tendency to merge to the IR equilibrium faster as Ξ 0 → 0, which is again in line with the argument given at the end of last section.
Conclusions. We outlined a new method of decomposing the distribution function in terms of non-hydrodynamic modes in a suitable basis, which cast the Boltzmann equation into a dynamical system for a non-equilibrium scale-invariant expanding plasma. We found the UV and IR fixed points of this system and examined its flow lines. It was shown that the exponentially small corrections tracing the nonlinear deformation history of the these flows eventually give a prescription for the renormalization of transport coefficients via transasymptotic matching [36]. The short-lived non-hydrodynamic modes give rise to non-Newtonian characteristics of the renormalized transport coefficients, and as a result, hydrodynamics becomes valid in far-from-equilibrium regimes where Kn and Re −1 are large. Therefore, we conclude that the success of hydrodynamics in high energy nuclear collisions is intrinsically related with the transient rheological behavior of the fluid.

SUPPLEMENTAL MATERIAL: EVOLUTION EQUATION OF THE MOMENTS c nl
In this supplemental material we briefly expose the derivation of the evolution equation of the moments c nl in (3). c nl can be conveniently written as c nl (τ) = s n, l d 3 p A n, 2l (τ) f τ, p T , p ς , where s n,l and A n, 2l are respectively given by s n,l = 1 4π Γ(n + 1) (4l + 1) /Γ (n + 4) , (14a) Next, by taking the time derivative of Eq. (13) one obtainṡ c nl = d 3 p s n, 2l A n, 2lḟ + d 3 p s n, 2lȦn, 2l f .