Dynamical systems and nonlinear transient rheology of the far-from-equilibrium Bjorken flow

TL;DR

This paper develops a non-perturbative dynamical renormalization scheme for far-from-equilibrium relativistic flows, revealing new non-hydrodynamic modes and providing insights into the nonlinear evolution of distribution functions.

Contribution

It introduces a transseries-based dynamical renormalization framework for kinetic theory, extending beyond linear response and uncovering new high-energy modes.

Findings

Existence of a Lyapunov function for RG flows

Transport coefficients are dynamically renormalized with non-perturbative corrections

Identification of a new high energy non-hydrodynamic mode

Abstract

In relativistic kinetic theory, the one-particle distribution function is approximated by an asymptotic perturbative power series in Knudsen number which is divergent. For the Bjorken flow, we expand the distribution function in terms of its moments and study their nonlinear evolution equations. The resulting coupled dynamical system can be solved for each moment consistently using a multi-parameter transseries which makes the constitutive relations inherit the same structure. A new non-perturbative dynamical renormalization scheme is born out of this formalism that goes beyond the linear response theory. We show that there is a Lyapunov function, aka dynamical potential, which is, in general, a function of the moments and time satisfying Lyapunov stability conditions along RG flows connected to the asymptotic hydrodynamic fixed point. As a result, the transport coefficients get dynamically renormalized at every order in the time-dependent perturbative expansion by receiving non-perturbative corrections present in the transseries. The connection between the integration constants and the UV data is discussed using the language of dynamical systems. Furthermore, we show that the first dissipative correction in the Knudsen number to the distribution function is not only determined by the known effective shear viscous term but also a new high energy non-hydrodynamic mode. It is demonstrated that the survival of this new mode is intrinsically related to the nonlinear mode-to-mode coupling with the shear viscous term. Finally, we comment on some possible phenomenological applications of the proposed non-hydrodynamic transport theory.

Figures

Equations

$${\partial }_{\tau }f\left(\tau ,p_T,p_{\varsigma }\right)=-\frac{1}{{\tau }_r(\tau )}\left[f\left(\tau ,p_T,p_{\varsigma }\right)-f_{eq.}\left(-u\cdot p/T\right)\right],$$

$${\tau }_r=\frac{{\theta }_0}{T^{1-\Delta }}.$$

$$f\left(\tau ,p_T,p_{\varsigma }\right)=f_{eq.}\left(\frac{p^{\tau }}{T}\right)\left[\sum _{n=0}^{N_n}\sum _{l=0}^{N_l}{c}_{nl}(\tau ){\mathcal{P}}_{2l}\left(\frac{p_{\varsigma }}{\tau p^{\tau }}\right){\mathcal{L}}_n^{(3)}\left(\frac{p^{\tau }}{T}\right)\right],$$

$$\displaystyle c_{nl}(\tau)=\,2\pi^{2}\,\frac{(4l+1)}{T^{4}(\tau)}\,\frac{\Gamma(n+1)}{\Gamma(n+4)}\,\biggl{\langle}\,\left(p^{\tau}\right)^{2}\,P_{2l}\left(\frac{p_{\varsigma}}{\tau\,p^{\tau}}\right)\,\mathcal{L}_{n}^{(3)}\left(\frac{p^{\tau}}{T}\right)\,\biggr{\rangle}\,,$$

$$T^{\mu \nu }=(ϵ+P_0)u^{\mu }{u}^{\nu }+g^{\mu \nu }{P}_0+{\pi }^{\mu \nu },$$

$$\begin{array}{rl}{T}^{\mu \nu }& =ϵu^{\mu }{u}^{\nu }+P_L{l}^{\mu }{l}^{\mu }+P_T{\Xi }^{\mu \nu },\end{array}$$

$$ϵ=⟨{\left(-u\cdot p\right)}^2⟩$$

$$\displaystyle P_{T}=\biggl{\langle}\,\frac{1}{2}\Xi^{\mu\nu}p_{\mu}p_{\nu}\,\biggr{\rangle}$$

$$\displaystyle P_{L}=\biggl{\langle}\,\left(l\cdot p\right)^{2}\,\biggr{\rangle}$$

$$-2.5\le c_{01}\le 5.$$

$$\bar{\pi }$$

$$\frac{1}{T}\frac{dT}{d\tau }+\frac{1}{3\tau }=-\frac{c_{01}}{30\tau },$$

$$\frac{d{c}_{nl}}{d\tau }+\frac{1}{\tau }\left[{\alpha }_{nl}{c}_{nl+1}+{\beta }_{nl}{c}_{nl}+{\gamma }_{nl}{c}_{nl-1}-n\left({\rho }_l{c}_{n-1l+1}+{\psi }_l{c}_{n-1l}+{\varphi }_l{c}_{n-1l-1}\right)\right]+\frac{1}{{\tau }_r(\tau )}\left(c_{nl}-{\delta }_{n,0}{\delta }_{l,0}\right)=0,$$

$$\begin{array}{rl}& {\alpha }_{nl}=\frac{(2+2l)(1+2l)(n+1-2l)}{(4l+3)(4l+5)},{\beta }_{nl}=\frac{2l(2l+1)(5+2n)}{3(4l+3)(4l-1)}-\frac{(4+n)}{30}{c}_{01},{\varphi }_l=\frac{(2l)(2l-1)}{(4l-3)(4l-1)},,\\ & {\gamma }_{nl}=(2l+2+n)\frac{(2l)(2l-1)}{(4l-3)(4l-1)},{\psi }_l=\frac{1}{3}\left(\frac{4l(2l+1)}{(4l+3)(4l-1)}\right)-\frac{c_{01}}{30},{\rho }_l=\frac{(2l+1)(2l+2)}{(4l+3)(4l+5)},\end{array}$$

$$\frac{d\mathbf{c}}{dw}=\mathbf{f}(w,\mathbf{c}),$$

$$\mathbf{f}(w,\mathbf{c})=-\left[\hat{\Lambda }\mathbf{c}+\frac{1}{w}\left(\mathfrak{B}\mathbf{c}+\mathbf{A}\right)\right]+\mathcal{O}({\mathbf{c}}^2,\mathbf{c}/w^2),(w\to \infty ,\mathbf{c}\to 0)$$

$$\mathbf{c}=(c_{01},\dots ,c_{0L},c_{10},c_{11},\dots ,c_{1L},\dots ,c_{N0},\dots ,c_{NL}),$$

$$\tilde{\mathbf{c}}(w)=U\mathbf{c}(w),\hat{\mathfrak{B}}=U\mathfrak{B}{U}^{-1}=\mathrm{diag}(b_1,\dots ,b_L)\in {\mathbb{C}}^L,\tilde{\mathbf{A}}=U\mathbf{A},\hat{\Lambda }=U\hat{\Lambda }{U}^{-1}.$$

$$\frac{d\tilde{\mathbf{c}}}{dw}=\tilde{\mathbf{f}}(w,\tilde{\mathbf{c}}),$$

$$\tilde{\mathbf{f}}(w,\mathbf{c})=-\left[\hat{\Lambda }\tilde{\mathbf{c}}+\frac{1}{w}\left(\hat{\mathfrak{B}}\mathbf{c}+\tilde{\mathbf{A}}\right)\right]+\mathcal{O}({\tilde{\mathbf{c}}}^2,\tilde{\mathbf{c}}/w^2).(w\to \infty ,\tilde{\mathbf{c}}\to 0)$$

$${\tilde{c}}_i(w)=\sum _{|\mathbf{m}|\ge 0}^{\infty }\sum _{k=0}^{\infty }{\tilde{u}}_{i,k}^{(\mathbf{m})}{E}_k^{(\mathbf{m})}(w),$$

$$E_k^{(\mathbf{m})}(w)={\mathbf{\sigma }}^{\mathbf{m}}{\mathbf{\zeta }}^{\mathbf{m}}(w)w^{-k},$$

$${\mathbf{\zeta }}^{\mathbf{m}}(w)=e^{-(\mathbf{m}\cdot \mathbf{S})w}{w}^{\mathbf{m}\cdot \tilde{\mathbf{b}}}=\prod _{i=1}^I{\left[{\zeta }_i(w)\right]}^{m_i},$$

$${\zeta }_i(w)=e^{-S_{i}w}{w}^{{\tilde{b}}_i},$$

$${\mathbf{\sigma }}^{\mathbf{m}}=\prod _{i=1}^I{\sigma }_i^{m_i},$$

$${\sigma }_i={\sigma }_j^{\ast },$$

$${\sigma }_i\in \mathbb{R},$$

$$\frac{d\log T}{d\log \tau }$$

$$\frac{d{c}_{0l}}{dw}$$

$$\frac{d\mathbf{c}}{dw}=\mathbf{f}(w,\mathbf{c}),$$

$$\mathbf{f}(w,\mathbf{c})=-\frac{1}{1-\frac{c_1}{20}}\left[\hat{\Lambda }\mathbf{c}+\frac{1}{w}\left(\mathfrak{B}\mathbf{c}-\frac{c_1}{5}\mathbf{c}+\mathbf{A}\right)\right],$$

$$\mathbf{c}=(c_{01},c_{02},\dots ,c_{0L-1},c_{0L}{)}^{\top },$$

$$\mathbf{A}=\frac{3}{2}{\left({\gamma }_{01},0,\dots ,0\right)}^{\top },$$

$$\hat{\Lambda }=\mathrm{diag}\left(\frac{3}{2{\theta }_0},\dots ,\frac{3}{2{\theta }_0}\right),$$

$$\mathfrak{B}=\frac{3}{2}\left(\begin{array}{cccccc}\frac{2}{3}{\Omega }_1& {\alpha }_{01}& & & & \\ {\gamma }_{02}& \frac{2}{3}{\Omega }_2& {\alpha }_{02}& & & \\ & {\gamma }_{03}& \frac{2}{3}{\Omega }_3& {\alpha }_{03}& & \\ & & \ddots & \ddots & \ddots & \\ & & & {\gamma }_{0L-1}& \frac{2}{3}{\Omega }_{L-1}& {\alpha }_{0L-1}\\ & & & & {\gamma }_{0L}& \frac{2}{3}{\Omega }_L\end{array}\right),$$

$${\Omega }_l=\frac{5l(2l+1)}{(4l+3)(4l-1)}.$$

$$\frac{d\tilde{\mathbf{c}}}{dw}=\tilde{\mathbf{f}}(w,\tilde{\mathbf{c}}),$$

$$\tilde{\mathbf{f}}(w,\tilde{\mathbf{c}})=-\frac{1}{1-\frac{c_1}{20}}\left[\hat{\Lambda }\tilde{\mathbf{c}}+\frac{1}{w}\left(\hat{\mathfrak{B}}\tilde{\mathbf{c}}-\frac{c_1}{5}\tilde{\mathbf{c}}+\tilde{\mathbf{A}}\right)\right],$$

$$c_i=\sum _{l^{\prime }=1}^L{U}_{i{i}^{\prime }}^{-1}{\tilde{c}}_{i^{\prime }}.$$

Full text

dynamical systems, resurgence, transseries, rheology, non-Newtonian fluids, relativistic kinetic theory

I Introduction

Recent measurements of multiparticle cumulants in $d+Au$ , $p+Pb$, $p+Au$, and ${}^{3}He+Au$ systems have provided a compelling evidence for collective flow in heavy-light ions collisions Chatrchyan et al. (2013); Adare et al. (2013); Aad et al. (2013); Abelev et al. (2013); PHE (2018). Similar observations had been made previously in nucleus-nucleus collisions (cf. Refs. Gale et al. (2013); Heinz and Snellings (2013) and references therein). Altogether, these measurements seem to imply that the collective evolution of a nuclear fireball created in these experiments can be understood in terms of hydrodynamics.

The phenomenological success of hydrodynamic models in producing an accurate description of these extreme experimental situations have raised some questions as per the validity of the assumptions that hydrodynamics relies on. For instance, the local thermal equilibrium hypothesis has been usually thought of as a necessary and sufficient criterion for the applicability of the fluid-dynamical equations of motion. Different theoretical studies Kurkela and Zhu (2015); Critelli et al. (2017); Denicol et al. (2014a); Florkowski et al. (2013a, b); Denicol et al. (2014b); Chesler and Yaffe (2010); Heller et al. (2012); Wu and Romatschke (2011); van der Schee (2013); Chesler (2016); Casalderrey-Solana et al. (2013); Martinez and Strickland (2010) have shown that the equilibrium condition seems to be too restrictive, and that surprisingly no inconsistency arises when hydrodynamics is used for interpreting a system sitting far from equilibrium. These findings have repeatedly pointed out to the possibility of generalizing relativistic hydrodynamic theories for systems with the local thermal equilibrium removed Romatschke (2018a); Florkowski et al. (2018); Romatschke and Romatschke (2017). Although the idea of employing fluid dynamics in non-thermal equilibrium physics is not entirely new Laghaei et al. (2005); Eu (1983), little progress has been made to formulate it from first principles.

In recent years we have learned more about certain generic properties of far-from-equilibrium hydrodynamics that follows suitably from the nonlinear nature of fluid dynamics. The breakdown of the naive perturbative asymptotic expansion is, therefore, imminent, meaning that a series ansatz fails to converge and accordingly it cannot be a full-blown solution to the underlying nonlinear differential equations. This has been long known in mathematics Costin (1998). An example of this fact in hydrodynamics was given in Heller and Spalinski (2015) where the gradient expansion of the energy-momentum tensor of a conformal fluid was shown to be divergent. The divergences of the energy-momentum tensor ought to beg for the existence of additional degrees of freedom called non-hydrodynamic modes. Since then more examples have appeared in generic setups in both strong and weakly coupled regimes  Casalderrey-Solana et al. (2018); Noronha and Denicol (2015); Bazow et al. (2016a); Strickland et al. (2018); Denicol and Noronha (2018a, b); Romatschke (2018a); Buchel et al. (2016); Spaliński (2018); Kurkela and Mazeliauskas (2018a, b); Heller et al. (2016); Heller and Svensson (2018); Heller et al. (2014); Heller and Spalinski (2015).

A more relevant and important question to ask in this context has to do with how non-hydrodynamic modes can affect the transport properties of the system and what their possible phenomenological significance is. We addressed this question in our previous work Behtash et al. (2018a) by investigating the nonlinear dynamics of a far-from-equilibrium weakly coupled plasma undergoing Bjorken expansion. Our work first treated the equations governing the time evolution of momentum moments of the one-particle distribution function as a coupled nonlinear dynamical system, the solutions of which are found to be multi-parameter transseries with real integration constants $\sigma$. These transseries carry exponentially small (non-perturbative) information in terms of the Knudsen $Kn$ and inverse Reynolds $Re^{-1}$ numbers which both play the role of expansion parameters in the transseries, whereas the relaxation time $\tau_{r}$ characterizes the size of exponential corrections. A linear combination of these moments defines the constitutive relations at every order in the perturbative expansion whose coefficients are indeed related to the transport coefficients. Then summing over all exponentially small factors and other monomials of the expansion parameters in the transseris renders an effective dynamical renormalization of the transport coefficients away from the values obtained at equilibrium using linear response theory. It should be noted that in this approach, ‘dynamical renormalization’ means that there exists a positive-definite monotonically decreasing function(al) also known as Lyapunov function in the stability theory, along every flow line in the space of moments, that can be easily derived from our dynamical system.

One advantage that comes with studying dynamical systems is the ability to have full control over the late-time (or IR) as well as early-time (or UV) behavior of the flow lines. We should keep in mind that the alternative techniques such as Borel resummation heavily rely on the asymptotics of the divergent series. However, in a dynamical system of evolving moments, the resurgent multi-parameter transseries is directly found as a formal solution without the need for going to Borel plane where UV information is completely lost. The Borel singularities in this sense only encode IR data such as the size of exponential corrections, which in our approach corresponds to the eigenvalues of the linearization matrix commonly known as Lyapunov exponents around the asymptotic hydrodynamic fixed point. Therefore, the theory after Borel resummation is only asymptotically equivalent to the original theory, and drawing any conclusions on the fate of solutions of the original theory in the UV from a resumed series is, in general, not possible.

The organization of this paper is as follows. We extend and generalize our previous findings Behtash et al. (2018a) to include higher non-hydrodynamic modes in Sects. II-III while providing technical details of our derivations in the Appendices. In Sect. IV, we explain further the deep relationship between dynamical systems and the resurgence theory in the time-dependent (non-autonomous) system at hand, which was first discussed in the context of relativistic hydrodynamics in Ref. Behtash et al. (2018b). Finally, in Sect. V it will be reported that a new high energy non-hydrodynamic mode exists whose first-order perturbative asymptotic term goes like $\sim(\tau T)^{-1}$, which obviously displays the same decay as the Navier-Stokes (NS) shear viscous tensor component. The survival of this new non-hydrodynamic mode is due to the nonlinear mode-to-mode coupling with the effective shear viscous tensor. Some of the quantitative and qualitative properties of this new mode will also be discussed. As a last couple of remarks, we should remind the readers that for the sake of keeping the paper self-contained, we have included a small number of mathematical definitions from dynamical systems in App. B that will make the read easier. Also, as mentioned above, we will adopt the field theory language when talking about the early-time and late-time regimes, which then can be interchangeable with UV and IR, respectively.

II Kinetic model

The statistical and transport properties of the weakly coupled systems are usually described by kinetic theory. Within this approach, it becomes important to understand the evolution of one-particle distribution function from the Boltzmann equation de Groot et al. (1980a); Cercignani and Medeiros Kremer (2002); Cercignani (1988). Fluid-dynamical equations of motion arise as a coarse-graining process of the distribution function where the fastest degrees of freedom are integrated out. The approximations made in the coarse-graining procedure do not necessarily lead to the same macroscopic evolution equations Denicol et al. (2012a, b), and nor do they in general draw an accurate picture of the physics when compared against exact numerical solutions Denicol et al. (2014b, a); Florkowski et al. (2013a, b); Florkowski et al. (2014); Martinez et al. (2017). The common approaches to finding approximate solutions of the Boltzmann equation are the Chapman-Enskog approximation  Chapman et al. (1990) and Grad’s moment method Grad (1949); de Groot et al. (1980a). In the Chapman-Enskog method, the distribution function is expanded as a power series in the Knudsen number around its equilibrium value, which eventually leads to a divergent series Santos et al. (1986); Denicol and Noronha (2016). In Grad’s moment method, the distribution function is expanded about the equilibrium state in terms of an orthogonal set of polynomials, and the coefficients of such expansion are average momentum moments Grad (1949). The second approach, however, turns out to be more convenient when dealing with far-from-equilibrium backgrounds Bazow et al. (2013); Molnar et al. (2016); Bluhm and Schäfer (2015); Martinez and Strickland (2010, 2011).

In this work, we also adopt the latter method and study the full nonlinear aspects of the moments’ evolution equations. We carry out our research plan by considering a weakly coupled relativistic system of massless particles which undergoes Bjorken expansion Bjorken (1983) with vanishing chemical potential. We simplify our problem by assuming the Boltzmann equation within the relaxation time approximation (RTA-BE).

For the Bjorken flow we use the Milne coordinates $x^{\mu}=(\tau,x,y,\varsigma)$ (with the longitudinal proper time $\tau=\sqrt{t^{2}-z^{2}}$ and pseudorapidity $\varsigma={\rm arctanh}\left(z/t\right)$), and the metric is $g_{\mu\nu}=\text{diag.}\left(-1,1,1,\tau^{2}\right)$. The time-like normal vector identified with the fluid velocity is taken to be $u^{\mu}=(1,0,0,0)$ (with $u^{\mu}u_{\mu}=-1$) and the spatial-like normal vector pointing along the $\varsigma$-direction is $l^{\mu}=(0,0,0,1)$ (with $l^{\mu}l_{\mu}=1$). The Bjorken flow symmetries reduce the RTA-BE to the following relaxation type equation Baym (1984)

[TABLE]

where $p_{T}$ is the transverse momentum, $p_{\varsigma}$ denotes the momentum component along $\varsigma$-direction, and $\tau_{r}$ represents the relaxation time scale. Also, $f_{eq.}$ is the local equilibrium Jüttner distribution which without loss of generality is chosen to be of Maxwell–Boltzmann type, i.e., $f_{eq}(x)=e^{-x}$. We recall that $\tau_{r}$ sets the time scale at which the system relaxes to its thermal equilibrium. We shall consider models whose relaxation time is a power law in the effective temperature, namely

[TABLE]

For $\Delta=1$ means that the constant $\theta_{0}$ is dimensionful, while for $\Delta=0$ the theory is conformally invariant.111A slightly more general class of models for the relaxation time approximation have been studied in Ref. Dusling et al. (2010). For pedagogical purposes we perform explicit calculations for the case of $\Delta=0$ in the main body of the paper, while the gist of results for the general case are discussed in Appendices  D, E and F.

The RTA-BE (1) can be solved exactly Baym (1984). In what follows, we will take a distinct approach in which the mathematical problem of solving the Boltzmann equation is recast into seeking solutions to a set of nonlinear ODEs for the moments. For the Bjorken flow we propose the following ansatz for the single-particle distribution function Behtash et al. (2018a)

[TABLE]

where $p^{\tau}=\sqrt{p_{T}^{2}+(p_{\varsigma}/\tau)^{2}}$ is the energy of the particle in the comoving frame, $\mathcal{L}_{n}^{(3)}$ and $P_{2l}$ denote the generalized Laguerre and the Legendre polynomials, respectively. The ansatz (3) allows us to study hydrodynamization processes Behtash et al. (2018a). Furthermore, the nonlinear relaxation of the low energy ($n=0$) as well as the high energy tails ($n>0$) of the distribution function are better understood when $f(\tau,p_{T},p_{\varsigma})$ is expanded in terms of orthogonal polynomials Bazow et al. (2016b, 2015).

The moments $c_{nl}$ are read directly from Eq. (3) as follows 222Blaizot and Li studied the time evolution of similar moments for a constant relaxation time  Blaizot and Yan (2017a) and a more general nonlinear collisional kernel in the small angle approximation Blaizot and Yan (2017b). In their case, the authors were interested in the details of the longitudinal momentum anisotropy which in our notation this corresponds to those moments with $n=0$. Up to some normalization factor, the main difference between the moments $\mathcal{L}_{l}$ (see Eq. (2.8) in Ref. Blaizot and Yan (2017a)) and our moments $c_{0l}$ is that the former are dimensionful.

[TABLE]

where $\Gamma(n)$ is the Gamma function and the momentum average of any observable $\mathcal{O}(x^{\mu},p^{\mu})$ weighted by an arbitrary distribution function $f_{X}$ is denoted as $\langle\,\mathcal{O}\,\rangle_{X}\equiv\int_{\bf p}\mathcal{O}(x^{\mu},p^{\mu})f_{X}(x^{\mu},p_{i})$ with $\int_{\bf p}\equiv\int d^{2}p_{T}dp_{\varsigma}/[(2\pi)^{3}\,\tau\,p^{\tau}]$. If $f_{eq.}(x)=e^{-x}$, then have that the hydrodynamic equilibrium (asymptotic IR fixed point) is given by $c_{nl}^{eq.}=\delta_{n0}\delta_{l0}$.

For the Bjorken flow the energy momentum tensor $T^{\mu\nu}=\langle\,p^{\mu}\,p^{\nu}\,\rangle$ is Molnar et al. (2016); Molnár et al. (2016); Huang et al. (2010, 2011); Gedalin (1991); Gedalin and Oiberman (1995) 333It is customary to use the following tensor decomposition of the energy-momentum tensor

$T^{\mu\nu}=(\epsilon+P_{0})\,u^{\mu}u^{\nu}\,+\,g^{\mu\nu}\,P_{0}\,+\,\pi^{\mu\nu}\,,$

where $P_{0}$ is the equilibrium pressure and $\pi^{\mu\nu}$ is the shear viscous tensor. Nonetheless, for highly anisotropic systems, Mólnar et. al Molnar et al. (2016); Molnár et al. (2016) showed that using $T^{\mu\nu}$ (5) is convenient. It should be noted that both formulations are equivalent Molnar et al. (2016) and explicit examples for the Bjorken Molnár et al. (2016) and Gubser flows Martinez et al. (2017); Behtash et al. (2018b) have already been discussed in the literature.

[TABLE]

with the projector operator $\Xi^{\mu\nu}=g^{\mu\nu}+u^{\mu}u^{\nu}-l^{\mu}l^{\nu}$ which is orthogonal to both $u^{\mu}$ and $l^{\mu}$. The energy density $\epsilon$, the transverse and longitudinal pressures ($P_{L}$ and $P_{T}$ respectively) can be written in terms of the moments $c_{nl}$ Behtash et al. (2018a)

[TABLE]

It follows that $\epsilon=2P_{T}+P_{L}$ from these expressions. We also mention the nonnegativity of pressure components, $P_{T},P_{L}\geq 0$, sets the physical range for $c_{01}$ as

[TABLE]

These bounds are satisfied by the exact solution of the RTA-BE (1) but it is not expected to be satisfied by a particular truncation scheme for the distribution function.

The energy-momentum conservation together with the Landau matching condition for the energy density imply $c_{00}\equiv 1$. The only independent (normalized) shear viscous component $\bar{\pi}\equiv\tau^{2}\,\pi^{\varsigma\varsigma}/\epsilon$ for the Bjorken flow is proportional to the moment $c_{01}$ Behtash et al. (2018a)

[TABLE]

In a previous work of some of us Behtash et al. (2018a), it was shown that the time evolution of the temperature $T$ and the moments $c_{nl}$ with $n\geq 0$ and $l\geq 1$ is described by an infinite number of coupled nonlinear ODEs as follows

[TABLE]

where the coefficients are given by

[TABLE]

The hierarchy of equations in (9) constitute a dynamical system where moments of different degree $n$ and $l$ mix amongselves in a nontrivial way. The nonlinear nature of the RTA-BE is manifest in the set of ODEs (9b) by the mode-to-mode coupling term $\sim c_{nl}\,c_{01}$. Nonetheless, one observes that the low energy modes $c_{0l}$ decouple entirely from the high energy ones $c_{nl}$ with $n>0$. As a result, the time evolution of the energy momentum tensor is fully reconstructed from the solutions of the temperature and the modes $c_{0l}$. Notice that one cannot deduce the same thing about the high energy modes whose evolution receive major contribution from the lower energy modes. Furthermore, the high energy moments shall play a role in the stability of the system and its convergence to its asymptotic thermal state, which is the subject of discussion in Sect. V.

III Transseries solutions and Costin’s formula

In this section, we construct the transseries solution to the dynamical system in (9). Once one reduces to the Boltzmann equation to a dynamical system, next immediate step would be to construct an exact form of the transseries around each asymptotic fixed point and obtain all the coefficients recursively from the evolution equations. In what follows, we will attempt to build the general form of the exact transseries solution by slightly modifying the original set of ODEs. For technical reasons, we also start with a truncation of the dynamical system at $0\leq l\leq L$ and $0\leq n\leq N$, hence the original Boltzmann distribution will be reproduced by taking the limit $N,L\rightarrow\infty$. Since we are interested in building the solutions starting at late times, we will assume the convergence at infinity, i.e., $c(\tau)\rightarrow 0$ for $\tau\rightarrow\infty$.

Let us prepare the set of ODEs in the dynamical system (9) in the following asymptotically linearized form:

[TABLE]

where $\hat{\Lambda}$ and ${\mathfrak{B}}$ are constant matrices, ${\bf c}$ is an $((N+1)\times(L+1)-1)$-dimensional vector given by

[TABLE]

and ${\bf A}$ is a constant vector. Here, we have defined a new time coordinate as $w=\tau T(\tau)$ which behaves asymptotically like $w\sim\tau^{2/3}$ at late times. We remind that ${\rm rank}(\hat{\Lambda})={\rm rank}({\mathfrak{B}})=(N+1)\times(L+1)-1=:I$. Furthermore, we assume $\hat{\Lambda}$ is a diagonal matrix proportial to the unit matrix. To construct the transseries solution, we suitably diagonalize ${\mathfrak{B}}$ by defining an invertible matrix $U$ such that

[TABLE]

These transformations cast Eqs. (11) in the following form

[TABLE]

We hereby call the components of the vectors $\tilde{\bf c}$ pseudomodes due to being generally complex-valued, thus not physical. This makes the components of the matrix $U$ complex-valued as well. But the inverse transformations of $\tilde{\bf c}$ achieved by the action of $U^{-1}$ always yield real vectors that eventually contribute to the observables in our theory.

The classical asymptotics beyond naive asymptotic power series expansion can be carried out with the help of a transseries ansatz first put forward by O. Costin in Ref. Costin (1998). In that seminal work, the author proves that the set of ODEs written in the prepared form (14) has the exact transseries solution

[TABLE]

where $I={\rm dim}\,({\bf c})$, ${\bf m}\in{\mathbb{N}}_{0}^{I}$ is an integer vector, and the dot denotes the inner product between any two vectors. The real numbers $\sigma_{i}$ are going be referred to as “integration constants” throughout this work as they would really symbolize the constants to be obtained if we were to integrate the ODEs in Eqs. (14). Here, for simplicity we have defined $E^{({\bf m})}_{k}(w)$ that stand for the basis of transmonomials (i.e., the exponential factors and the fractional powers $w^{-1}$).

The transseries data such as the coefficients $\tilde{u}^{({\bf n})}_{i,k}$, the Lyapunov exponents $S_{i}$ and the anomalous dimensions $\tilde{b}_{i}$ can be recursively determined by the evolution equations up to normalization of $\sigma_{i}$. Without loss of generality, we pick the normalization fixed by $u_{j,0}^{({\bf m})}=\delta_{ij}$ for $m_{j}=\delta_{ij}$. This leaves no room for ambiguity in determining the coefficients $\tilde{u}_{i,k}^{({\bf m})}$. It is noteworthy that the transseries solution of $c_{i}(w)$ can be reproduced by the inverse transformations of (13), and one can find that $c_{i}(w)$ has essentially the same transseries as $\tilde{c}_{i}(w)$ due to the fact that the matrix $U$ acts on the index $i$ (mode number) in $\tilde{u}^{({\bf m})}_{i,k}$ only.

Although the transseries of $\tilde{c}_{i}(w)$ given by (15) is generally complex-valued due to $\tilde{b}_{i}$ and $\tilde{u}^{({\bf 0})}_{ik}$ taking values in complex numbers, we can still recover a real transseries solution for the moments $c_{i}(w)$. The important fact is that if $\tilde{b}_{i}$ for some $i$ is complex, it is always accompanied by its complex conjugate counterpart, i.e., $\tilde{b}_{i}=\tilde{b}^{*}_{j}$ for some $j$, where the associated coefficients satisfy $\tilde{u}^{({\bf m})}_{i,k}=\tilde{u}^{({\bf m})*}_{j,k}$ under a certain normalization of the eigenvectors of $\hat{\mathfrak{B}}$ in such a way that $U^{-1}_{Ii}=1$ for every $i$ where again $I$ is the rank of matrix $\hat{\mathfrak{B}}$. Therefore, the reality of $c_{i}(w)$ is guaranteed if the following conditions are satisfied:

[TABLE]

III.1 Evolution equations: $N=0$ case

In terms of $w$, Eqs. (9) can be reduced to the following non-autonomous dynamical system:

[TABLE]

The second equation is valid for any $l>0$. Since temperature $T$ is now sort of washed away from Eq. (22b), we shall solve these ODEs for $c_{nl}$ by the mathematical techniques discussed in Refs. Iwano (1982); Costin (1998). These tools are well known in the context of resurgence theory. A curious reader is invited to check e.g., Ref. Costin (2008) for technical details, Ref. Dorigoni (2014) for a nice introduction to the subject, and Ref. Dunne and Ünsal (2016) for the summary of its applications to quantum mechanics and quantum field theories). We will afterwards substitute $c_{01}$ in Eq. (22a) to solve for $T$.

It is convenient to cast Eqs. (22a)-(22b) in the following familiar form

[TABLE]

where

[TABLE]

Here, the index of vectors and matrices runs over $i=1,\dots,L$. To directly apply Costin’s formula, one has to diagonalize ${\mathfrak{B}}$ in Eq. (23) using the transformations in (13) to get the equation

[TABLE]

We can easily see that $c_{0i}$ can be reproduced by $\tilde{\bf c}$ and $U$ as

[TABLE]

As is seen in Eq. (31), from now on whenever $N=0$, we may drop the index [math] in $c_{0l}$ for ease of notation. Notice that the asymptotic expansion of $\tilde{c}_{i},c_{i}$ both start at order $\mathcal{O}(w^{-1})$ asymptotically, and due to this we have $\tilde{u}^{(0)}_{i,0}=0$.

Before trying to solve the dynamical system in (23), it would be instructive to first linearize the r.h.s. of it to obtain the so-called IR data $S_{i}$, $\tilde{b}_{i}$ in Eq. (18) of the transseries ansatz. These are going to be the input data obtained uniquely from the profile of the solution of the linarized equation around the asymptotic fixed point of the dynamical system at $w\rightarrow\infty$. To do so, we expand the factor $\frac{1}{1-c_{1}/20}$ to get

[TABLE]

where ${\bf c}_{1}(w)=\left(c_{1}(w),0,\dots,0\right)^{\top}$. The coefficients of the lowest power of $w^{-1}$, i.e., $u^{({\bf 0})}_{i,k}$, can be found from this equation, and one can readily see that $c_{i=l}=\mathcal{O}(w^{-l})$. In particular,

[TABLE]

Now, we may proceed to linearize the equation (III.1) by substituting $c_{i}\rightarrow\bar{c}_{i}+\delta c_{i}$ with $\bar{c}_{i}=\sum_{k}u^{({\bf 0})}_{i,k}w^{-k}$, and expanding to first order in $\delta c_{i}$

[TABLE]

Using the matrix $U$, the solution to the linearized equation is found to be given by

[TABLE]

where $\sigma_{i}$ is the integration constant and $\delta\tilde{\bf c}=U\delta{\bf c}$. Therefore, we have the IR data (Lyapunov exponent and anomalous dimension) as

[TABLE]

Finally, we substitute the ansatz (15) in Eq. (32) and then insert the resulting transseries for $\bf c$ in Eqs.(30)-(31) to find the recursive relation for the transseries coefficients as

[TABLE]

where $u^{({\bf m}^{\prime})}_{1,k^{\prime}}=\sum_{i=1}^{L}U^{-1}_{1i}\tilde{u}^{({\bf n}^{\prime})}_{i,k^{\prime}}.$ Fixing $\tilde{u}^{({\bf m})}_{i,0}=1$ for $m_{i^{\prime}}=\delta_{i,i^{\prime}}$, we can solve this equation order by order to get the transseries data $u^{({\bf n})}_{1,k}$, $S_{i}$ and $\tilde{b}_{i}$ without any (imaginary) ambiguity. Note that the IR data should match the numbers found using linearization around the asymptotic fixed point in Eq. (37).

III.2 Evolution equation: $N\neq 0$ case

The technique advocated in the $N=0$ case above can be extended to include higher order moments $c_{nl}$. The energy dependence of hydrodynamization is only materialized by $c_{nl}~{}(n=1,\dots,N)$ even if $c_{0l}$ have been for the most part sidelined in the literature. To get a complete and correct picture, however, all the moments ought to be taken into consideration. In principle, this means that one has to eventually solve an infinite-dimensional dynamical system for getting the full hydrodynamization process in both IR and UV regimes, which is obviously ambitious. As a result, a practical approach is to consider a truncated dynamical system whose proposed transseries solution has of course the ability to be generalized to the exact result. So our starting point in this section will be Eq. (22b),

[TABLE]

The equation (39) can be written in a concise matrix form,

[TABLE]

where the quantities ${\bf c}$, $\hat{\Lambda}$, ${\bf A}$ are explicitly given by

[TABLE]

and we have defined the block matrices ${\mathfrak{B}}$ and ${\mathfrak{D}}$ as

[TABLE]

with

[TABLE]

Here, the blocks $\bar{\mathfrak{B}}_{00}$, $\bar{\mathfrak{D}}_{00}$ are both $L\times L$ matrices, $\bar{\mathfrak{B}}_{10}$,$\bar{\mathfrak{D}}_{10}$ are $(L+1)\times L$ matrices, and the remaining blocks are $(L+1)\times(L+1)$ matrices. As in the $N=0$ case, we will be defining the matrix $U$ to diagonalize ${\mathfrak{B}}$,

[TABLE]

We finally employ the transseries ansatz in the diagonalized form of Eqs. (40)-(41) as before to obtain the recursive relations involving the transseries data, namely

[TABLE]

where $\tilde{\mathfrak{D}}=U{\mathfrak{D}}U^{-1}$, $\tilde{c}_{i}(w)=\sum_{|{\bf m}|\geq 0}^{\infty}\sum_{k=0}^{\infty}\tilde{u}^{({\bf m})}_{i,k}E^{({\bf m})}_{k}(w)$, and $\bf m$ is the index of each non-perturbative sector of the pseudomodes, $\tilde{c}_{i}(w)$. Eq. (51) can be solved recursively. The higher non-perturbative sectors are connected to the lower ones through mode-to-mode couplings of the form $c_{01}c_{nl}$, whereas operations in charge of promoting the asymptotic order are differentiation and multiplication by a $1/w$ factor. Once the perturbative order is fixed, the 1st-order non-perturbative sector can be constructed explicitly.

For pedagogic purposes, we show how to compute the leading order contribution in the perturbative sector of the transseries corresponding to $\bm{m}=\bm{0}$ (also known as the leading order in the asymptotic expansion). As an example, first, we determine which moments behave like $\mathcal{O}(1/w)$. This is effectively done by taking $k=0$ in (51). Since the asymptotic fixed point has all the moments other than $c_{00}$ go to [math] at $w\rightarrow\infty$, the coefficients $\tilde{u}^{\bm{0}}_{i,0}$ have to vanish. By default, all the other coefficients promoting the asymptotic order also vanish. For the lowest asymptotic behavior for each $n$ and $l$, Eq. (40) gives

[TABLE]

Here, we relabeled the index $i$ by $n,l$ for convenience. Since $c_{nl}\sim{\cal O}(1/w)$ and $A_{nl}\neq 0$ only when $(n,l)\in\{(0,1),(1,1)\}$, we find that both $c_{01}$ and $c_{11}$ have a non-zero coefficient for the 1st-order asymptotics given by

[TABLE]

Since $\bm{A}$ is constructed by $\gamma_{01}$ and $\phi_{1}$, only these two moments survive at the 1st asymptotic order, that is, $c_{01}$ (cf. Eq. (34) aka Navier-Stokes) and $c_{11}\sim-\phi_{1}\theta_{0}/w$. For the other moments, the leading asymptotic order is dictated by their neighbors $c_{n^{\prime}l^{\prime}}$ alongside mode-to-mode coupling. The equation responsible for identifying the leading order of $c_{nl}\sim\mathcal{O}(1/w^{k})\,\,(k>1)$ reads

[TABLE]

All the other operations promoting the asymptotic order vanish, such as differentiation. For example, the leading asymptotic order of $c_{2l}~{}(l=0,1,2)$ is identified by

[TABLE]

As can be seen, they both have the lowest asymptotic order $1/w^{2}$ because of $c_{11}$. By virtue of Eq. (51), the leading order asymptotics of higher moments is determined to be $c_{nl}\sim\mathcal{O}(1/w^{k})$ where $k=\max(n,l)$ for $n,l\geq 2$. The higher modes, therefore, decay faster in general as all the Lyapunov exponents are equal due to the RTA-BE having a single scale which in this case is $\theta_{0}$, and the suppression in the IR is just determined by the leading asymptotic order. All the modes satisfying $n,l\leq 1$, however, are exempt from this rule as $c_{10}\sim 1/w^{2}$, and as mentioned earlier, the two remaining lower modes $c_{01}$ and $c_{11}$ are the slowest modes of all in the IR.

Once the coefficients of the bare asymptotic series for each $c_{nl}$ are obtained, the transseries coefficients can be achieved by solving the equation including higher ${\bf m}>0$ non-perturbative sectors. As an example, we show the transseries and leading order bare asymptotics of five moments $c_{01},\dots,c_{21}$ of the truncated dynamical system $N=2,~{}L=1$ in Fig. 1  444In this work the initial condition for the distribution function is given by the RS ansatz of the distribution function Romatschke and Strickland (2003), i.e.,

$f_{0}(\tau,p_{T},p_{\varsigma})=\exp\left[-\frac{\sqrt{p_{T}^{2}+(1+\xi_{0})\left(\frac{p_{\varsigma}}{\tau}\right)^{2}}}{\Lambda_{0}}\right]\,.$

where $\Lambda_{0}$ and $-1<\xi_{0}<\infty$ are the initial effective temperature and momentum anisotropy along the $\varsigma$ direction. Thus the set of initial conditions for the moments $c_{nl}$ depend on the initial time $\tau_{0}$ and initial anisotropy parameter $\xi_{0}$. See App. E for further details. . Each moment is renormalized by the inclusion of non-perturbative/perturbative sectors available using the transasymptotics as $C_{nl,k}/w^{k}+\mathcal{O}(1/w^{k})$, where the transasymptotic matching is approximated up to $15$th-order transmonomials. On the numerical front, the integration constants $\bm{\sigma}$ are determined by employing simultaneous least-squares method, which aims to minimize the difference of the exact and transseries solutions of all the moments involved in a given truncated dynamical system simultaneously. So the overall average error will be reduced across all the moments. In this figure, five integration constants $\sigma_{1},\dots,\sigma_{5}$ are given by their optimized values along with individual standard deviations reported in the plots. Also, the light blue dashed lines stand for the transseries solution with the best optimized $\sigma_{i}$ while the blue shaded areas highlighting the possible variation of transseries due to standard deviation. Furthermore, bare asymptotic expansion at $1/w,~{}1/w^{2},$ and $1/w^{3}$ orders for each moment is plotted as a comparison. These plots show that continuing to a higher asymptotic order will not necessarily lead to a better result, as opposed to adding more transmonomials at each order, which results in a better approximation of the transasymptotic matching overall (see next section).

III.3 Transasymptotic matching

In this section, we construct the transasymptotic matching condition responsible for the full-blown form of the time evolution of the transport coefficients including all the effects of non-perturbative sectors, which is a 1st-order PDE. In doing so, we first redefine $\tilde{c}_{i,k}(w)$ as

[TABLE]

and sum over ${\bf m}$ in Eq. (51). Then the transasymptotic matching condition turns out to have the following formal form

[TABLE]

where $\hat{\zeta}_{i}:=\partial/\partial\log(\sigma_{i}\zeta_{i})$. Here, $\tilde{C}_{i,k}$ is supposed to be only a function of $\sigma_{i}\zeta_{i}$. We again note that Eq. (60) is a 1st-order PDE, whose solution yields the time evolution of the transport coefficients with all the $\sigma_{i}\zeta_{i}(w)$ included. Moreover, the integration constants may be determined from the transseries coefficients $\tilde{u}^{({\bf m})}_{i,k}$.

Due to its partial differential nature, it is technically not easy to solve (60) for the general case even though for $L=1,N=0$, it can be exactly solved. For instance, the first three coefficients are given by solving (60) in this case as Behtash et al. (2018a)

[TABLE]

where $W_{\zeta}:=W(-\sigma\zeta/20)$ is Lambert $W$-function.

For a system consisting of a larger number of moments, the above equation becomes a PDE, which renders the matching condition hard to solve. However, the transasymptotic matching, in this case, is at least effectively consolidated by summing over a large number of transmonomials just in the same manner as in the renormalization group equation for a running coupling constant, where loop diagram contributions are added to the r.h.s. in order to capture a more realistic RG flow. As a special example, we guide the readers to take a look at App. D.1 in which the transasymptotic matching for $N=0,~{}L=1,~{}\Delta=1$ system yields polynomials of finite degree, which are the exact solutions to the dynamical system. In Fig. 1, the transmonomials are taken up to 15th order to get as close as possible to the exact solution of the truncated system. We again mention that the renormalized moment is meant to be a solution compatible with the full-fledged transasymptotics at its leading asymptotic order.

III.4 Comment on resurgence and cancellation of the imaginary ambiguities

Once one has an ODE of the form (11), we can realize imaginary ambiguity cancellation in a systematic way. In general cases, it is hard to prove the cancellation of imaginary ambiguities. However, for ODEs of the type (11) there are mathematically rigorous proofs for the ambiguity cancellations due to O. Costin (cf. Costin (1998) and references within for more technical details). In this section, we consider an easy example involving only one mode $i=1$, and therefore we omit the index for simplicity.

Let us begin with a formal transseries ansatz given by

[TABLE]

Suppose that $F(w)$ solves (11) and let us only focus on the [math]th non-perturbative sector $a^{(0)}_{k}$. By calculating the radius of convergence one can see that the power series (64) is divergent. The approximate form of the growing upper bound is

[TABLE]

with some positive real constants, $M$ and $S$. This means that $\hat{\Phi}^{(0)}(w)$ is asymptotically of Gevrey-1 class 555The Gevrey-$n$ class means that there exists a smooth function $f$ on $\mathbb{R}^{d}$ such that on every compact subset $C$, there are constants $p,q$ such that $|D^{\alpha}f(x)|\leq pq^{k}(k!)^{n}$. Here, $D^{\alpha}$ is a differential operator of some multi-indices $\alpha$ such that $|\alpha|=k$., hence the Borel transform is just enough for convergence purposes, defined by

[TABLE]

Note that $\Phi^{(0)}(\xi)$ has a finite convergence radius and $r_{c}=S$. Since $\hat{\Phi}^{(0)}(w)$ is a divergent series, there should exist a branch-cut on the positive real axis with a branch-point $S$ on the Borel plane. The position of the first singularity corresponds to the Lyapunov exponent of the flow lines approaching the asymptotically stable fixed point in the dynamical system. Notice that the appearance of a branch-cut is a reminder that the series expansions are coupled to higher order transmonomials such as $\exp(-Sw)/w$ which independently satisfy the same ODE, i.e. Eq. (11).

The original divergent series (64) can be reproduced through the Laplace integral,

[TABLE]

and subsequently taking the asymptotic limit

[TABLE]

It is noteworthy that the imaginary ambiguity depending on the contour of integration is nothing but an artifact of going to the Borel plane and it is exponentially suppressed in the large-$w$ limit. However, when one observes a singularity on the Borel plane, one can build a relationship among non-perturbative sectors by taking the Hankel contour going around the singularity in the Laplace integral,

[TABLE]

The relations obtained in this way are called resurgence relations.

We now are in position to introduce the Stokes automorphism $\underline{{\mathfrak{S}}}^{\theta}$. Formally speaking, the integration constants in the transseries may jump once the singularity (Stokes) rays on the Borel plane are crossed. This is associated with the existence of an automorphism that takes one integration constant to another plus some new contribution, defined by

[TABLE]

where $\{\rho_{\theta}\}$ denotes the set of singular points along the Stokes ray with the contour enclosing the singularities being controlled by an angle parameter $\theta$. Here, $\Delta_{\rho}$ represents an abstract derivative operator known as alien derivative which can be understood as an infinitesimal change in the asymptotic behavior due to a singularity in the direction parametrized by $\theta$. Note that

[TABLE]

if $\hat{\Phi}^{(0)}(w)$ is not a divergent series along $\theta$ direction. However, when a singularity appears for a particular angle $\theta$, we can make the resurgence relations among different sectors. In our case, the set of singularities on the real positive axis can be given by $\{\rho_{0}\}=\{nS\,|\,n\in{\mathbb{N}}\}$, where again $n=1$ shows the Lyapunov exponent of the flow lines flowing to the asymptotic stable fixed point of the dynamical system. By taking the Hankel contour with a singularity at $nS$, one can obtain the so-called bridge equation

[TABLE]

and the relationships among different sectors are given by

[TABLE]

where $A\in i\,{\mathbb{R}}$ the famous Stokes constant. This equation means that information of some sector is carried over to higher sectors once Stokes rays are crossed.

Although the construction of resurgence relations is generally an independent issue from the imaginary ambiguity cancellation, we can systematically find such a cancellation mechanism for the transseries solution of the ODE in (11) without the need to resort to the discussion of this section. The significantly important fact is that a nonzero r.h.s. Eq. (76), satisfies the same ODE. Therefore, the Borel-summed transseries along the $\theta=0$ ray may well be expressed by the arbitrariness in the integration constant, and the fact that nothing prevents one from extending $\sigma\in{\mathbb{R}}$ to $\sigma\in{\mathbb{C}}$. The relationship between two integration constants across a Stokes ray is then given by

[TABLE]

Therefore, the imaginary ambiguity can be cancelled by shifting the integration constant in such a way that

[TABLE]

Let us demonstrate the imaginary ambiguity cancellation for simple cases. For $L=1$ and $N=0$, the leading large-order growth of the asymptotic expansion is given by Basar and Dunne (2015)

[TABLE]

where $S_{1}=\frac{3}{2\theta_{0}},\tilde{b}_{1}=-18/35$, and $C(\theta_{0})$ is an overall factor depending on $\theta_{0}$. Here, we have used $\tilde{u}_{1,0}^{(1)}=1$ in the second line. Since $w$ always appears with a constant factor $1/\theta_{0}$ in $\tilde{c}_{l}(w)$, one finds that

[TABLE]

where $C_{0}$ is a constant. We measured $C_{0}$ by numerical fitting with the a prepared form (79) and obtained

[TABLE]

Therefore, the Stokes constant $A_{1}$ is related to $C(\theta)$ in the following way

[TABLE]

leading to the imaginary ambiguity cancellation once we set

[TABLE]

For general $L$ and $N$, $\tilde{u}^{(0)}_{l,k}$ are complex-valued so are the integration constants $\sigma_{l}$ for the reality condition to hold. For $L=2$, the overall factor in (79) is found to be

[TABLE]

in which $\tilde{b}_{1}=(75+i\sqrt{10655})/110=\tilde{b}_{2}^{*}$. By fitting numerically we can estimate the complex constant $C_{0}$ as

[TABLE]

Therefore, the Stokes constant as a function of $\theta_{0}$ can be fully evaluated from Eq. (82). Finally, the imaginary ambiguity cancellation follows by just shifting the integration constant as

[TABLE]

IV Global dynamics and asymptotic stability analysis

In this section we will explore the generic properties of the Bjorken flow from the perspective of dynamical systems with some useful definitions provided in App. B for readability.

IV.1 Phase space and asymptotic UV/IR fixed points

The explicit time-dependence of (9) means that the system is non-autonomous and the phase space is larger than that of time-independent (autonomous) system with the same number of variables. Let us assume that the total number of Legendre and Laguerre terms in the expansion of the distribution function is $L+1$ and $N+1$, respectively. Then one can write down the phase space of this truncated dynamical system as

[TABLE]

where each subspace $M_{i}$ consists of $L$ dimensions, and we excluded the constant component $c_{00}=1$ from the phase space. So (87) is just the hyperspace $c_{00}=1$ of the full phase space. Eventually, the limit $N,L\rightarrow\infty$ has to be taken, meaning that the true physical dynamical system at its full glory is infinite-dimensional 666This is indeed expected by noting that the distribution function $f$ solving the Boltzmann equation has an exponential kernel which admits infinitely many terms in its series expansion, leading to an infinite number of moments.. The tuple $(M_{0},\dots)$ is just short for the product of all the entries, and the time manifold $\mathbb{t}$ is taken to be the set of positive real numbers $\mathbb{R}^{+}$. This choice is made due to the fact that Eq. (9a) is not regular at $\tau=0$ so to keep the continuity, $\tau<0$ is not allowed. Since temperature $T$ has its own equation in (9), we also add its manifold $\mathbb{T}$ to the crowd and let it take values over $\mathbb{T}=\mathbb{R}^{+}$. In total we are then left with $(L+1)(N+1)+1+1-1=(L+1)(N+1)+1$ dimensions for $\mathcal{M}$.

By looking at the case where $N=0$, the dynamical system reduces to $L+2$ dimensional subspace $(M_{0},\mathbb{T},\mathbb{t})$ which is denoted by $\mathcal{M}_{0}\subset\mathcal{M}$. The special facet of this subspace is that it indeed encompasses the invariant manifold 777This manifold is basically an embedding in an $L+2$ dimensional phase space. Because the exact solution of RTA-BE has only the parameters $\xi_{0}$, $T_{0}$, and $\tau_{0}$ controlling the shape of flow lines, the invariant manifold has to be a three dimensional object. For the argument in the $w$ coordinate, cf. Sect. IV.2. because its equations are completely decoupled from the rest $(n>0)$ in case we started with a larger system involving more moments, i.e., $N>0$. Therefore, any solution of this sub-system entirely lies inside $\mathcal{M}_{0}$ at all times $\tau>0$. The second observation is that the mode-to-mode coupling terms in $\eqref{eq:ODEsBjor}$ for $n>0$ all depend on the elements of the subspace $\mathcal{M}_{0}$. So it is crucial to understand the dynamics happening in this subspace by first locating the asymptotic fixed points and analyzing their stability, and then obtaining the flow lines connecting these fixed points to find the qualitative shape of the invariant manifold.

In terms of the variable $w=\tau T(\tau)$, the non-autonomous dynamical system is described by

[TABLE]

where $\mathbb{w}=\mathbb{R}_{+}:=\mathbb{R}^{+}\cup\{0\}$. Note that $w=0$ has been allowed since the combination of $\tau T(\tau)$ does essentially have a well-defined zero limit. The invariant manifold now is inside $W_{0}$ parametrized by ${c}_{0l}$. The structure of this invariant manifold will be the subject of next subsection.

There are however three major differences between the two parametrizations of the phase space of the Bjorken flow. Namely,

• 1-

In terms of $\tau$, the early-time $\tau<1$ or UV limit of the theory does have a pair of fixed points for the truncation order $L,N=\infty$; one being the maximally oblate point at which the longitudinal pressure $p_{L}$ vanishes and transverse pressure $p_{T}$ is maximized, and one more fixed point, at which $p_{L}$ is maximum and transverse pressure is vanishing, to be called maximally prolate point. The early-time shape of flow lines in $\mathcal{W}$ looks completely different by allowing a continuous flow from the maximally prolate point, while in $\mathcal{M}$ there is no such flow in general since the flows cannot reach $T=0$ line on which both these points are located. Either one of the two parametrizations, nonetheless, leads to the same IR structure for flow lines in both phase spaces $\mathcal{W}$ and $\mathcal{M}$.

• 2-

A nonphysical singularity hypersurface $c_{01}=20$ emerges in $\mathcal{W}$ that again affects the UV physics of the flow lines initiated in the basin of attraction of the invariant manifold admitting the range given in Eq. (7).

• 3-

In the original version of the dynamical system, Eqs. (9), the temperature does not converge in the UV so practically speaking it is impossible to connect continuously to the maximally oblate fixed point at $T=0$ by running the evolution equations backward in time, say, from somewhere in the IR. It is also not possible to find a complete flow connecting to another UV fixed point present in the Bjorken flow phase space $\mathcal{M}$ Behtash et al. (2018a). In terms of $w$, the UV structure of the theory is altered through a change in the stability of fixed points at $w\rightarrow 0$ even though the IR will remain intact as mentioned above. It is now possible to search for flow lines starting at the maximally prolate point (see Fig 3 for instance), but still it will not be feasible to have a flow that begins its journey at exactly the maximally oblate fixed point. This latter solution, if existed, would be a critical line due to the extreme fine-tuning required around a saddle point to initiate a structurally unstable flow on the boundary of the invariant manifold that connects to the IR fixed point $c_{nl}=\delta_{n0}\delta_{l0}$ at $w\rightarrow\infty$. This is explained below.

If we consider the full theory with $L=\infty$ and $N=0$ in the phase space $\mathcal{M}$ of the Bjorken flow, there are two UV fixed points in general whose coordinates are given by

[TABLE]

where we have adopted the notation $(\tau,T,\{c_{0l},c_{02},\dots\})$. In Fig. 2 the stability analysis of these two fixed points in the truncated system $N=0,L=31$ as well as $N=30,L=31$ is shown. At a general truncation order, the stability analysis is summarized as follows for any relaxation time $\tau_{r}=\frac{\theta_{0}}{T^{1-\Delta}}$ ($0\leq\Delta\leq 1$):

1. (a)

$N=0,~{}L=\text{odd}$: The system has 2 fixed points in the UV. The maximally prolate fixed point is a source while the maximally oblate fixed point is a saddle point with one repelling direction, i.e., index $1$ fixed point); 2. (b)

$N=0,~{}L=\text{even}$: The system has only an index $1$ maximally oblate point. The total index is still conserved since the other fixed point becomes a complex saddle; 3. (c)

$N>0,~{}L=\text{odd}$: The maximally prolate fixed point becomes an $\text{index}(I-N+1)$ saddle point whereas the maximally oblate fixed point is an $\text{index}(N+1)$ saddle point; 4. (d)

$N>0,~{}L=\text{even}$: The maximally oblate fixed point is a saddle of $\text{index}(N+1)$.

Here, we have again used the notation $I=(L+1)(N+1)-1$. In $L={\rm odd}$ case the total UV index is always $I+2$. Considering that an asymptotically stable hydrodynamic fixed point has zero index, then the total Morse index is $I+2$. A comparison between the numerical solutions of the flow lines in $L={\rm odd}$ system and the exact RTA-BE shows that $I+2$ is the correct topological invariant of the real phase portrait of the RTA-BE. Note that the temperature is always an attractive direction in both UV and IR, and therefore it does not contribute to the index. Furthermore, due to the lack of one real fixed point (maximally prolate point) in the UV when $L=\rm{even}$, this kind of truncation will not be physically sensible and will be hereby omitted from our discussions.888Again, the past of the dynamical system in $w$ coordinate is altered from the original case in a way that $w=0$ does not capture the physics of original theory, which is an argument related to the stability of individual fixed points. Nonetheless, the conclusion about the index is true across the board. The phase portraits of different 3d sectors of the dynamical system in (39) at odd truncation orders are plotted in Figs. 3 and 4.

An interesting question that comes to mind concerns the existence of a flow line that connects the maximally oblate fixed point to the asymptotic hydrodynamic fixed point in the phase portrait of the Bjorken dynamical system using $w$ coordinate. There are two things that one might keep in mind here. First, the invariant manifold is in the subspace $W_{0}\subset\mathcal{W}$ or simply the subspace of the phase space with coordinates $c_{0l}$. Since the maximally oblate point for $N>0$ is a saddle point of index $N+1$ (see (c) above), it must be located on the boundary of the invariant manifold. So at least in the $N>0,~{}L=\text{odd}$ system, there cannot be a critical line since it will always flow outside of the basin of attraction of the invariant manifold due to the repelling directions being in the subspace $\bigcup_{n=0}^{N}W_{n}$. Second, for $N=0$, the maximally oblate point is essentially a sink (time is the only repelling direction). As a result, the $c_{0l}$ set initially to be the coordinates of this fixed point at $w=0$ will change its position with time, while being a fixed point at any moment $w>0$ until eventually it becomes the hydrodynamic fixed point at $w\rightarrow\infty$. So there will not be a flow line (process) to begin with, proving that the critical line does not exist in the $w$ parametrization of the Bjorken flow just like it could not exist in the $\tau$ coordinate but for a completely different reason.999Another way to phrase this is that the r.h.s. of the dynamical system in Eq. (23) is obviously zero at the maximally oblate fixed point, but it does remain zero at all times $w>0$ even if the position of this fixed point changes. This means that $dc_{0l}(w)/dw=0$ so a flow will not happen to exist at any later time if ${\bf c}_{0}$ are equal to the values given in (89).

IV.2 Initial value problem for transseries and invariant manifold

Although the transseries contains a large amount of information about a given dynamical system, fixing the integration constants is an important, yet challenging, problem. On the one hand, since our (truncated) dynamical system has $I$ variables-thus $I$ first-order ODEs, it consequently has $I$ integration constants needed to be tuned in order for the basin of attraction (of the invariant space) to be determined. On the other hand, the initial condition of the exact integral solution of the Boltzmann equation has only three parameters, namely $\tau_{0},T_{0}$, and $\xi_{0}$.

By redefining the ODEs using $w$ as a time coordinate, $T_{0}$ can be completely separated from the rest of the dynamical system which now has the form (22a). The new initial time $w_{0}$ determines a point on a flow (more precisely a process), hence $\sigma_{i}$ do only depend on the anisotropy parameter $\xi_{0}$. We could behold from the numerical results of the integral solution of the Boltzmann equation that every flow is closed on an open subset of the $\mathcal{W}$ as long as $\xi_{0}$ is chosen from the domain $[-1,+\infty)$. Moreover, we recall that if a flow is initiated in the subspace $W_{0}\subset\mathcal{W}$ parametrized by $c_{0l}$, it always stays in that subspace, meaning that the invariant manifold has to be an embedding in the space $(W_{0},\mathbb{w})$. But how does transseries know about the invariant manifold? Since $\xi_{0}$ determines the basin of attraction of the isolated invariant space of the Boltzmann equation, there has to be a map from the $\xi_{0}$-space to $W_{0}$ such that

[TABLE]

where ${\mathbb{w}}=[0,+\infty)$ is the time manifold and ${\cal F}$ is the space of integral curves given by the $n=0$ ODEs in the dynamical system. Because of the uniqueness theorem, the map in (91) is a bijection only when $L,N=\infty$, and consequently one can find that the invariant manifold is a two dimensional surface embedded in the space $\cal F$.

It is a quite challenging problem to figure out the correspondence defined by (92), but we can narrow ${\mathbb{R}}^{I}$ down to a subspace formed by those $\bm{\sigma}(\xi_{0})$ which yield flow lines starting in the basin of attraction of the invariant manifold. This means that $\bm{\sigma}({\xi_{0})}$ is able to give the invariant manifold once the stability analysis of each fixed point in the UV is performed and the approximate boundary of this space is determined by finding possible critical lines. To facilitate the search for the invariant manifold, one can construct the transseries around every fixed point individually at $w\rightarrow 0$, which would be something for the form

[TABLE]

where the anomalous dimensions $\tilde{b}_{i}$ of the pseudomodes are now the eigenvalues of the linearization matrix about the corresponding fixed point.101010Since pseudomodes are not physical, $\tilde{b}_{i}$ can take complex values in general. Choosing $\sigma_{i}=0$ for every $i$ for which $\Re(\tilde{b}_{i})<0$, we get a maximally prolate point that is a source (that is otherwise a spiral) fixed point. We will discuss the stability of the UV fixed points in Sect. IV.1.

An immediate question that comes to mind is whether one can construct a complete flow line and/or a critical line from the transseries solutions around individual fixed points. One part of the answer involves topological arguments to be discussed in an upcoming work from the perspective of dynamical systems and equivariant cohomology Conley index Behtash et al. (work in progress). Another piece comes from transferring the information from one transseries to another by tuning $\sigma_{i}$ to a particular set of numbers from both UV and IR sides, provided that the flow line achieved this way remains on the invariant manifold at all times.

As a final remark, let us mention that one can analytically continue the transseries by rotating in the time axis $w\rightarrow e^{i\alpha}w$ such that transseries solution around the asymptotic hydrodynamic fixed point is now exactly a Fourier decomposition at $\alpha=\pi/2$. Neglecting an overall imaginary factor at each non-perturbative order, it is easy to see that $\sigma_{i}$ serve as amplitudes of the individual wave components in the Fourier decomposition. The significant advantage of analytic continuation is that the attractor entails oscillatory components which were exponentially suppressed in real time, thus fitting the transseries to the complex numerical solutions of the analytically continued Boltzmann equation is much more tractable.

IV.3 Dynamical renormalization of 2nd order transport coefficients

In this section, we aim to construct an RG equation for the transport coefficients in analogy with the gradient descent approach to the RG flows in the context of quantum field theories using the language of dynamical systems (cf. Ref. Gukov (2017)).

The variable $w$ so far has been playing the role of flow ‘time’ in our system of ODEs, i.e., Eq. (22b). But in this section, we want to interpret it as playing the role of ‘energy scale’ for the renormalization scheme arising from resuming all the non-perturbative fluctuations around the asymptotic expansion of non-hydrodynamic modes ${\bf c}:=\{c_{nl}\}$. The reason is simple: the Knudsen number $Kn=\tau_{r}(\tau)|D_{\mu}u^{\mu}|\sim w^{-1}$ in Bjorken flow provides the only parameter by which the system evolves from the UV regime all the way to IR. Since the transseries coefficients of individual moments include transport coefficients with exponentially suppressed factors in $w$ following Behtash et al. (2018a), the far-from-equilibrium dynamics of the moments can then be associated to the running of these coefficients as $w$ changes. Hence, preparing a time-dependent (non-autonomous) dynamical system for $\bf c$ automatically amounts to having a renormalization group equation on the phase space of moments $\bf c$ and $w$, together denoted by $\mathcal{W}$, the so-called phase space of the Bjorken flow in $w$ parametrization as in (88).

A flow (process) on the phase space is described by a continuous map

[TABLE]

such that ${\phi}^{w,w_{0}}(\mathbf{c}_{0})=\mathbf{c}(w)$, and ${\phi}^{w,w_{0}}({\phi}^{s,w_{0}}(\mathbf{c}))={\phi}^{w+s,w_{0}}(\mathbf{c})$ where $w\in\mathbb{R}^{+}$ is the RG time. The IR regime of the theory is captured by the behavior of the flow lines approaching to an asymptotic stable fixed point aka hydrodynamic fixed point satisfying $c_{l}^{eq}=\delta_{0l}\in\mathcal{M}$ for which ${\phi}^{w,w_{0}}(c_{l}^{eq})=c_{l}^{eq}$ for all $w\in\mathbb{R}^{+}$. We alternatively refer to $c_{l}^{eq}$ as the asymptotic IR fixed point since it is reached at $w\rightarrow\infty$.

Let us formulate our dynamical system as

[TABLE]

where $\mu:=w/\theta$ is the RG scale and $\log(\mu)\in\mathbb{R}$. As $\mu$ varies, the moments ${\bf c}$ are mapped to themselves in a self-similar way, and any initial RG scale $\mu_{0}$ can initiate a process to a later $\mu$ by the renormalization group action. Since the r.h.s. only depends on ${\bf c}$ and $\mu$, one can write

[TABLE]

This is nothing but an RG equation in the space of moments $\bf c$ and the vector $\bm{\beta}(\mu,{\bf c})$ consists of $\beta$-functions, each encoding the dependence of every non-hydrodynamic mode on the renormalization scale, $\mu$, in the process of equilibration. In this expression, it is convenient to redefine our notations as

[TABLE]

Note that the ordinary derivative along an RG flow can be expressed by two partial derivatives as

[TABLE]

In the formalism we have been seeking to build so far, one prominent assumption is that an observable ${\cal O}$ is able to be expressed in terms of $c_{i}(w)$ (or $\tilde{c}_{i}(w)$), namely ${\cal O}={\cal O}({\bf c}(w))$. Therefore by solving the RG equation via transseries, we can achieve a renormalized form of the observable, say a transport coefficient, up to an arbitrary order of our choosing.

It is straightforward to define the RG equation by starting with the ODE in (40). To make things more tractable, we again go back to the RG time $w$ and compute the derivative of an observable ${\cal O}({\bf c}(w))$ w.r.t. $\log w$

[TABLE]

where $\tilde{\bf f}=U{\bf f}$ and $\tilde{\bm{\beta}}:=w\tilde{\bf f}$. Here, use was made of $d/d\log w=(\tilde{\bf b}-{\bf S}w)\cdot\hat{\zeta}+\partial/\partial\log w$ and the fact that $\tilde{c}_{i}(w)=\sum_{k=1}^{\infty}\tilde{C}_{i,k}(\bm{\sigma\zeta}(w))w^{-k}$. We should also point out that $\hat{\zeta}_{i}:=\partial/\partial\log(\sigma_{i}\zeta_{i})$. The first line describes the scaling behavior of ${\cal O}$ in terms of $w$, and the second line contains information of the dynamical system through $\bm{\beta}$. Hence, solving the RG equation determines a renormalization of the observable ${\cal O}$. Since we want to consider the RG equation for the transport coefficients, we are slightly better off with the definition of ${\cal O}$ changed to ${\cal O}={\cal O}({\bf C}_{k}(w))$. Roughly speaking, this implies that ${\cal O}$ is a function of the transport coefficients and depends only on $\bm{\zeta}(w)$. Consequently, the scaling behavior could be derived in a fashion similar to what we did to obtain (99), so

[TABLE]

However, the link between the scaling behavior and the dynamical information hidden in $\bm{\beta}$ is very non-trivial. Yet, one can proceed to calculate the RG equation for the $L=1,\,N=0$ system.

To do so, we employ the scaling behavior in (100) which gives

[TABLE]

Suppose now that $\zeta_{1}$ and $w$ are independent of each other, i.e., $\tilde{c}_{1}=\tilde{c}_{1}(\sigma_{1}\zeta_{1},w)$ and $\tilde{C}_{1,k}=\tilde{C}_{1,k}(\sigma_{1}\zeta_{1})$. In this case we have

[TABLE]

To remove $w$ from the l.h.s. of Eq. (102), we take a contour integration around the origin after multiplying by $w^{k-1}$, which then gives

[TABLE]

Hence, Eq. (100) for $L=1$ and $N=0$ reads

[TABLE]

where $\hat{\zeta}_{1}\tilde{C}_{1,k}(\sigma_{1}\zeta_{1})$ is given by Eq. (103).

The definition in (104) has a direct connection with the dynamical information encoded in $\beta_{1}$, so that (104) is regarded as the RG equation of the transport coefficients by setting ${\cal O}(C_{1,k}(w))=C_{1,k}(w)$. In the $L=1,N=0$ system, choosing $k=1$ in Eq. (104) and solving for the renormalized quantity $C_{1,1}$ gives that it is proportional to the shear-to-entropy ratio at the equilibrium Behtash et al. (2018a). Plugging $C_{1,1}$ back in the expansion of the non-hydrodynamic moment $c_{1}$ obtained using the linear response theory (see App. A) automatically promotes the 1st-order transport coefficient to the non-equilibrium case compatible with the transasymptotic matching along an RG flow approaching the IR fixed point asymptotically. This essentially gives a dynamically renormalized $\eta/s$ Behtash et al. (2018a) as follows

[TABLE]

where $-\frac{3}{40}C_{1,1}(w\rightarrow\infty)\equiv(\eta/s)_{0}$ Behtash et al. (2018a). Likewise, the 2nd order renormalized transport coefficients are given by

[TABLE]

where again $\frac{9}{40}C_{1,2}(w\rightarrow\infty)\equiv T\left[\tau_{\pi}\left(\frac{\eta}{s}\right)_{0}-\left(\frac{\lambda_{1}}{s}\right)_{0}\right]$.

We conclude this section by discussing whether there is some sort of a Lyapunov function that would identify the RG flows from a global dynamical standpoint. It is well-known that one can formally define an object analogous to a dynamical potential from the Newtonian mechanics. For simplicity, we assume that $\beta$ is independent of $w$ (a scenario that would hold once the limit $\theta_{0}\rightarrow\infty$ is taken).111111When we keep $w$ in the $\beta$-function, we have to promote the non-autonomous system to an autonomous system of one dimension higher by introducing an ODE for $w$ in terms of a new flow time $\rho$ as follows

$\displaystyle\frac{d{\bf c}(\rho)}{d\log\rho}=\bm{\beta}({\bf c}(\rho),w(\rho)),\quad\frac{dw(\rho)}{d\log\rho}=\beta_{w}(w(\rho)).$

(109)

Now, we define a positive-definite differentiable function $V$

[TABLE]

We can easily show that $V$ is a monotonically decreasing function in terms of $w$ as

[TABLE]

thus a candidate for a global Lyapunov function satisfying the conditions in Eq. (140).121212In case $V$ is not positive definite, one can always add a positive constant to it such that it remains positive at all times $w<\infty$.

V Hydrodynamization of soft and hard modes

In previous sections, we demonstrated that the moments $c_{nl}$ solving the ODEs in the dynamical system of (22b) turn out to be of the multi-parameter transseries form. In this section, we continue our quest for what type of information one can get from these solutions on a more physical ground by analyzing the late-time behavior of different momentum and energy sectors of the distribution function. In doing so, we stumble upon something interesting: a flow line in the space of moments initiated at an arbitrary state in the basin of attraction of the asymptotically stable fixed point, is controlled at late times by not only the known non-hydrodynamic moment $c_{01}$ but also the mode $c_{11}$. The latter also happens to possess the same perturbative decay as $c_{01}$. We recall that this immediately leads to the statement that the distribution function and any observable projecting onto/or involving at least the $l=1$ sector of the distribution function would receive two major IR contributions. In this section, we shall analyze in great details the impact of the slowest ($=\mathcal{O}(1/(\tau T))$) non-hydrodynamic modes in the IR.

Let us start rewrite our ansatz as

[TABLE]

where the soft (s), semi-hard (sh), and hard (h) sectors of the distribution function are respectively defined by

[TABLE]

In order to assess which sector of the distribution function hydrodynamizes faster, we focus on the decay of the following normalized moments Denicol and Noronha (2016); Strickland (2018) 131313In this work we follow the notation of Strickland Strickland (2018).

[TABLE]

$\bar{M}^{nm}$ are then certainly sensitive to the energy and momentum tails of the distribution function. To this explicitly, first note that $\bar{M}^{nm}\geq 0$ if and only if $f(\tau,p_{T},p_{\varsigma})\geq 0$. This condition holds for the exact solution of the RTA-BE, but not necessarily for the approximate distribution function. Second, the Landau matching condition for the energy density implies $\bar{M}^{20}\equiv 1$. Last, some of the moments in (114) such as $\bar{M}^{10}=n/n_{eq.}$ (, where $n$ and $n_{eq.}$ are the non-equilibrium and thermodynamic particle densities, respectively) and $\bar{M}^{01}=P_{L}/P_{0}$ have a direct interpretation in terms of the usual macroscopic variables.

After taking the time $\bar{M}^{nm}$ in Eq. (114), we end up with the following equation

[TABLE]

The solutions to this system indicate that different momentum tails of the distribution function reach values at equilibrium state asymptotically, that is the asymptotic hydrodynamic fixed point. In order to put the problem into the framework of dynamical systems, we write $\bar{M}^{nm}$ (114) as a linear combination of the moments $c_{nl}$ using the ansatz (3) of the distribution function, i.e.141414In order to get this expression, we utilize the following identities Gradshteyn and Ryzhik (2007); Abramowitz and Stegun (1964)

$\displaystyle\int_{0}^{\infty}\,dx\,e^{-x}\,x^{\gamma-1}\,\mathcal{L}^{(\mu)}_{n}(t)=\frac{\Gamma\left(\gamma\right)\,\Gamma\left(1+\mu+n-\gamma\right)}{n!\,\Gamma\left(1+\mu-\gamma\right)}\,\,,\hskip 28.45274pt\text{Re}[\gamma]>0\,,$

$\displaystyle\int_{0}^{1}\,dx\,x^{\rho}\,\mathcal{P}_{m}(x)=\frac{\pi^{1/2}}{2^{\rho+1}}\,\frac{\Gamma\left(\rho+1\right)}{\Gamma\left(1+(\rho-m)/2\right)\,\Gamma\left(m/2+3/2+\rho/2\right)}\,,\hskip 28.45274pt\text{Re}[\rho]>-1\,.$

[TABLE]

where $\Gamma(x)$ is the Gamma function and $\binom{a}{b}$ is the binomial coefficient. Therefore the solutions to the dynamical system in (115) are also written in the form of multi-parameter transseries.151515Alternatively one can also find a compact form of the moment $\bar{M}^{nm}$ from the exact RTA-BE solution of the distribution function (see Eq. (129) in App. F). These solutions flow to the stable fixed point at late times once the initial state is chosen in the basin of attraction of the invariant manifold (of the truncated system), which indicates that each momentum sector is guaranteed to equilibrate eventually. In other words, $\bar{M}^{nm}\to 1$ since $c_{nl}^{eq.}\to\delta_{n0}\delta_{l0}$ asymptotically.

In Sect. III.2, (cf. Eq. (53)), we showed that the slowest non-hydrodynamic moments decay perturbatively like $c_{01},c_{11}\sim 1/w$ and $c_{10},c_{20}\sim 1/w^{2}$ at large $w$. Thus, the asymptotic behavior of $\bar{M}^{nm}$ is dominated by the slowest moments appearing in each sector of the distribution function.

We study first the normalized moments $\bar{M}^{nm}$ for $n,l>0$ where there are contributions coming only from both soft and semi-hard sectors the distribution function whose slowest non-hydrodynamic modes are $c_{01}$ and $c_{11}$, respectively. In this case, the impact of these non-hydrodynamic modes is determined by taking the following asymptotic limits in Eq. (117):

• •

Soft (s) regime in which there is only the leading order asymptotic contribution of $c_{01}\sim 1/w$, i.e.,161616In Eq. (4.6) of Ref. Strickland (2018), the author assumes explicitly the 14-moment approximation when truncating the distribution function, that is

\begin{split}f&=f_{eq.}\left(1+\frac{p_{\mu}p_{\nu}\pi^{\mu\nu}}{2(\epsilon+p)T^{2}}\right)\,,\\ &=f_{eq.}\left(1+\frac{3}{16}\left[3\left(\frac{p_{\varsigma}}{\tau}\right)^{2}-\left(p^{\tau}\right)^{2}\right]\frac{\bar{\pi}}{T^{2}}\right)\,,\\ &\approx\,f_{eq.}\left(1+\frac{1}{3T^{2}}\left[\left(p^{\tau}\right)^{2}-3\left(\frac{p_{\varsigma}}{\tau}\right)^{2}\right]\frac{\eta}{s}\,\frac{1}{T\,\tau}\,\right)\,.\end{split}

By substituting this expression in our definition (114), one obtains the result derived by Strickland, Eq.(4.8) Strickland (2018), which obviously differs from ours (118). This discrepancy arises from the truncation procedure. Here, we do not truncate the momentum and energy dependence of the distribution function and rather keep these in their exact form. Instead, we truncate the number of moments entering the distribution function at our discretion. We remind the reader of the work of Denicol et al. Denicol et al. (2012b), where it was shown that the 14-moment approximation does not provide a unique set of equations of motion for the dissipative currents, an ambiguity which was resolved in Ref. Denicol et al. (2012a).

[TABLE]

• •

Soft + semi-hard ($s+sh$) regime which incorporates the leading order asymptotic contribution of both $c_{01}$ and $c_{11}$ up to $\mathcal{O}(w^{-1})$, namely

[TABLE]

In Fig. 5, a number of $\bar{M}^{nm}$ are plotted versus the variable $w$, obtained from the exact RTA-BE solution (200) (black line), the $s$ (118) (green dot-dashed line), and the $s+sh$ (119) (blue dashed line) regimes. We verify numerically that the IR behavior is not affected by the choice of initial conditions as long as they agree with the bounds set by the invariant manifold. In order to get a hold of how fast the normalized moments equilibrate, we set an arbitrary saturation bound at which a given moment falls within 5% of its equilibrium value, that is to say $M^{nm}=1\pm\delta$ with $\delta=0.05$, and the $+$ or $-$ is taken depending on whether the moment monotonically increases or decreases, respectively. In this setting, as seen in Fig. 5, it can be verified that independently of the truncation scheme, the normalized moments will asymptote to one, as expected. In the UV, however, the $s$ and $s+sh$ regimes disagree with the exact RTA-BE result. This should not come as a surprise since both limits are only valid when the distance from the hydrodynamic fixed point is small. But in the IR region, the moments $\bar{M}^{nm}_{s}$ first reach the saturation bound while $\bar{M}^{nm}_{s+sh}$ and $\bar{M}^{nm}_{exact}$ approach this value much later on, and above all else, almost at the same time. In the same limit, we behold a remarkable agreement between the exact and $s+sh$ solutions. Thus, the asymptotic behavior of the normalized moments $\bar{M}^{nm}$ ($n,l>0$) can only be given entirely by a linear combination of both non-hydrodynamic modes $c_{01}$ and $c_{11}$.

In Fig. 6, the value of $w=w_{c}$ is numerically evaluated. As is shown, if the indices for energy, $n$, and longitudinal momentum, $m$, are increased, $w_{s}$ for the exact and $s+sh$ moments marks a later time. The convergence of $\bar{M}^{nm}_{s}$ depends only on $m$, which is evident from Eq. (118). The numerical results in this figure attests that in fact the exact and $s+sh$ moments are in good agreement, and therefore $s+sh$ is the correct description of the IR regime of different energy and momentum sectors of the distribution function solely based on a linear combination of both $c_{01}$ and $c_{11}$ as written in Eq. (118) 171717A careful reader would notice a minor disagreement (less than 5$\%$) between the values of $w_{c}^{exact}$ and $w_{c}^{s+sh}$. This small mismatch is due to the difficulty in extracting numerically the large-time behavior of the exact RTA-BE solver used here. While concluding the draft of this paper, we were informed of an optimized algorithm that uses logarithmically-spaced grid in proper time and reduces the numerical error at late times. This new code is available in the ancillary files included in the arXiv version of Ref. Strickland (2018).. Given this agreement, one concludes that the larger the index $n>0$, the more relevant the inclusion of the new non-hydrodynamic mode, $c_{11}$. We want to mention that this is solely due to the power of dynamical systems since the $1/w$ decay of both $c_{01}$ and $c_{11}$ in the perturbation theory around the hydrodynamic fixed point is directly confirmed from the form of ODEs in (115).

We are now in position to analyze the normalized moments $\bar{M}^{nm}$ for $n>0,l=0$. In this case, only the hard sector of the distribution function contributes and the slowest modes playing a role in Eq. (114) are $c_{10}$ and $c_{20}$, both asymptotically decaying as $1/w^{2}$. As a double check, consider for instance the mode $\bar{M}^{10}=n/n_{eq.}$. For the Bjorken flow, the Chapman-Enskog expansion in the RTA approximation Jaiswal et al. (2015) shows that the NS corrections to the particle density $n$ vanishes as chemical potential $\mu$ is turned off, and thus the leading order dissipative corrections are of the second order in the Knudsen number $\sim 1/w^{2}$. In general, it is known that second-order fluid dynamics theories, e.g. Israel-Stewart theory, do not reproduce properly the behavior of the heat flow and particle density as well as the dominant corrections arising from couplings between the particle density and the shear stress tensor Jaiswal et al. (2015); Denicol et al. (2014c).

Now, as done before, we study the impact of the non-hydrodynamic modes entering Eq. (114) by taking two truncation limits of the hard ($h$) regime

[TABLE]

In Fig. 7 the evolution of several normalized moments $\bar{M}^{n0}$ are shown. In the IR $\bar{M}^{n0}_{c_{10}+c_{20}~{}(h)}$ and $\bar{M}^{n0}_{exact}$ reach the asymptotic hydrodynamic fixed point faster than the moments $\bar{M}^{n0}_{c_{10}~{}(h)}$, supporting the compatibility of the exact normalized moments with the former. Note that the asymptotic behavior of the normalized moments $M^{n0}$ is entirely determined by the linear combination of the non-hydrodynamic modes $c_{10}$ and $c_{20}$. This conclusion does not hold when $n\geq 3$, where only the mode $c_{10}$ contributes.

In Fig. 8 the values of $w_{c}$ vs. $n$ for the moments $M^{n0}$ are plotted. It is seen that the saturation bound is approached as $n$ increases if $n\geq 3$, as opposed to when $n=0,1$ for which $w_{c}$ decreases. Now, when comparing with the results plotted in Fig. 6, we observe that $\bar{M}^{nm}$ ($n>0,m\geq 1$) asymptote to their values at the equilibrium state later than the moments $\bar{M}^{n0}$ do. In terms of our orthogonal polynomial basis expansion, it is easy to understand this as the slowest non-hydrodynamic modes in the soft and semi-hard sectors decay like $\sim 1/w$, whereas in the hard sector, the slowest modes die away like $1/w^{2}$. On the other hand, for the given bound $\delta=0.05$, there is a disagreement between all the hard asymptotic limits (120) and the exact results due to the truncation limitations. As we emphasized previously, close to the IR fixed point, the inclusion of both moodes $c_{10}$ and $c_{20}$ (Eq. (120b)) is a must, which is also seen in Fig. 8.

The rest of this section is dedicated to showing how the presence of both $c_{01}$ and $c_{11}$ is a crucial factor in unveiling the asymptotic behavior of the exact distribution function in the IR. We calculate the deviations from the thermal equilibrium distribution function defined as

[TABLE]

following the procedure outlined in Ref. Heinz and Martinez (2015). In the $s$ and $s+sh$ regimes, $\delta f$ reads

[TABLE]

where $T_{NS}$ denotes the NS solution of the conservation law (9a). The Fig. 9 outlines the time evolution of $\delta f$ for the exact and asymptotic limits. As a double check, other valid initial conditions are considered by varying the anisotropy parameter $\xi_{0}$. Two remarks are due here. First, no momentum and energy sector of the distribution function does in fact thermalize homogeneously, and second, the low energy particles ($p^{\tau}/T<1$) thermalize faster than the highly energetic ($p^{\tau}/T>1$) particles. This is not surprising since it is known that for weakly coupled systems, and in the absence of external fields, soft particles equilibrate faster Bazow et al. (2016b, 2015); Ernst (1981); Krook and Wu (1976). On the one hand, the IR limit of $s$ truncation (122a) does not capture the relaxation of the energy and momentum tails when compared against the exact RTA-BE solution. On the other hand, the exact RTA-BE result is perfectly compatible with the $s+sh$ limit (122), which in turn confirms the necessity of including the new mode $c_{11}$. Therefore, the proper description of the IR in the Bjorken flow has to incorporate the dynamics of this non-hydrodynamic mode as well.

In sum, the numerical results presented in this section together with the previous analysis of the multi-parameter transseries solutions for the moments $c_{nl}$ lead to the true mechanism behind the nonlinear relaxation processes of different momentum sectors of the distribution function, which is completely governed by various mode-to-mode couplings among the moments in the framework of dynamical systems. Furthermore, the IR regime of the distribution function is determined uniquely by two non-hydrodynamic modes $c_{01}$ and $c_{11}$. The conclusions derived in this section would still hold if one used other relaxation time models (cf. App. F).

VI Conclusions

In this work, we have proposed a new dynamical renormalization scheme which allows us to study the nonlinear transient relaxation processes of a weakly-coupled plasma undergoing Bjorken expansion. The distribution function was expanded in terms of orthogonal polynomials. The nonlinear dynamics of the RTA-BE is then investigated through the kinetic equations of the coefficients entering this expansion, i.e. the average momentum moments $c_{nl}$ of the distribution function. We not only developed further our earlier findings Behtash et al. (2018a) to include higher modes, but we also found new interesting results summarized in the following

• •

Based on the seminal works of Costin Costin (1998); Costin and Costin (2001) we show that the coupled system of nonlinear ODEs for the moments $c_{nl}$ admit analytic multi-parameter transseries solutions. At every given order in the time-dependent perturbative asymptotic expansion of each mode, the summation over all the transient non-perturbative sectors appearing in the trans-series leads to renormalized transport coefficients. This presents a new description of the transport coefficients in the regimes far from equilibrium with an associated renormalization group equation, going beyond the usual linear response theory.

• •

As long as $0\leq\Delta<1$, the $T=0$ limit of the Boltzmann equation is explicitly time-independent, meaning that the distribution function $f$ does not evolve with time. In the language of our dynamical system, this limit is equivalent to saying that the system is autonomous. The stability of the maximally oblate point shows in this case that it is a sink and consequently, there is no flow line that could connect it to the asymptotic hydrodynamic fixed point at $\tau\rightarrow\infty$. As soon as $T>0$, we have that the dynamical system explicitly depends on $\tau>0$ (i.e., it is non-autonomous) and continuously connected to the IR theory (hydrodynamics) at $\tau\rightarrow\infty$. But now the maximally oblate point is out of reach, hence a flow line (process) cannot be found with $\tau_{0}=0,T_{0}=0$ which leads to hydrodynamical behavior at late times. In other words, the actual phase space $\mathcal{M}$ of the Bjorken flow is a disjoint union of hypersurfaces $T=0$ and $T\neq 0$.

• •

The above conclusion does not hold true in the $w=\tau T$ parametrization. The phase space $\mathcal{W}$ is different in the UV in the sense that the system stays always non-autonomous for all $w\geq 0$. Also, notice that the temperature has been washed away from $\cal W$ at the cost of introducing a new singularity hypersurface $c_{01}=20$. However, the maximally prolate point is the only fixed point in the UV that can connect to the IR as shown on the l.h.s. of Fig. 3. The maximally oblate point that in principle should correspond to the free-streaming limit, is a saddle point of index $N+1$. On the one hand, if $N=0$, it is basically a sink from the perspective of an observer sitting in the subspace $W_{0}\subset\mathcal{W}$ parametrized by $c_{0l}$. Therefore, if the initial conditions are set to be the coordinates of this fixed point, the observer only sees the maximally oblate fixed point moving to the IR as time passes. On the other hand, there is no bifurcation or change of stability along the way. Rather, the same fixed point at $w=0$ is just dislocated to become the fixed point at $w\rightarrow\infty$. As a result, there is no complete flow line connecting this fixed point to the hydrodynamic fixed point.

Also, if $N>0$, there appear to be repelling directions in the subspace $\bigcup_{n=1}^{N}W_{n}$. However, the invariant manifold is in $W_{0}$, and hence the search for a solution hits a snag again. If there was such a line, then we would have an example of a critical line. This critical flow line is often mistaken in the literature as a global attractor as if the system was autonomous. In non-autonomous systems such as the Bjorken flow, the attractor is a concept encoding information of the past or future of the system of equations governing its dynamics independently of each other, as shown in Fig. 3. The relevant attractor is then either a forward or pullback attractor. In sum, given our study of the dynamical system for the Bjorken flow, we conclude that there is no critical line connecting the maximally oblate point to the hydrodynamic fixed point, which is otherwise known as “attractor solution” in the hydrodynamics literature.

• •

For the Bjorken flow the asymptotic behavior of the kinetic equations unveils the existence of a new non-hydrodynamic mode whose decay rate is exactly the same as the usual NS shear viscous component, i.e., $\sim(\tau T)^{-1}$. The origin of the new mode’s decay is due to the nonlinear mode-to-mode coupling $~{}c_{11}c_{01}$ which dominates close to the asymptotic hydrodynamic fixed point. This non-hydrodynamic mode is the slowest high energy mode and it determines quantitatively the late-time behavior of the transient high energy tails of the distribution function.

In the present paper, we left out many interesting questions which certainly require appropriate answers. For instance, within our approach it is important to understand the analytic behavior of the retarded correlators together with the subtle interplay between branch cuts and poles as a function of the coupling Romatschke (2016); Kurkela and Wiedemann (2017). It will be also important to investigate the possible phenomenological consequences of the new non-hydrodynamic mode $c_{11}$. For instance, one might wonder about the impact of this high energy mode in the modeling of a jet crossing an expanding QGP (quark-gluon plasma) using kinetic theory. Another possibility is to understand if the origin of azimuthal anisotropies at large momenta is connected to $c_{11}$. It has been argued by different authors that the azimuthal anisotropies at intermediate $p_{T}$ are related to non-hydrodynamic transport  Romatschke (2018b); Borghini and Gombeaud (2011); Kurkela et al. (2018). On a more theoretical ground, it would be very interesting to come up with extending our analysis to the challenging case of nonlinear PDEs such as Israel-Stewart hydrodynamic equations for either $2+1$ or $3+1$ dimensions. Finally, our work can be generalized to studying nonlinear aspects of time-dependent systems of ODEs whose description is cast in the form of the properties of a dynamical system. We have left these exciting topics for future research projects.

Acknowledgements

We thank O. Costin for his patience, useful discussions, explanations and clarifications of some of the technical and conceptual matters in his groundbreaking work Costin (1998). We also thank G. Dunne, M. Spalinski, M. Heller, J. Casalderrey-Solana, T. Schäfer and M. Ünsal for useful discussions on the subject of resurgence and non-equilibrium physics. We are grateful to T. Schäfer for bringing to our attention the possible relation between our calculations and non-Newtonian fluids. H. S. thanks C. N. Cruz-Camacho for his help and assistance with the numerical solver of the RTA Boltzmann equation. M. M. thanks L. Yan for clarifying the details of Teaney and Yan (2014) and M. Strickland for his insistence on trying to understand the asymptotic behavior of the normalized moments $M^{n0}$. A. B., M. M., and S. K. are all supported in part by the US Department of Energy Grant No. DE-FG02-03ER41260. M. M. is partially supported by the BEST (Beam Energy Scan Theory) DOE Topical Collaboration. A. B. was partially supported by the DOE grant DE-SC0013036 and the National Science Foundation under Grant No. NSF PHY-1125915.

Appendix A Asymptotic Chapman-Enskog expansion of the distribution function

In this section, we briefly describe the asymptotic Chapman-Enskog expansion of the distribution function for the RTA Boltzmann equation. The Chapman-Enskog expansion has been calculated up to second order by different authors Teaney and Yan (2014); Chattopadhyay et al. (2015); Bhalerao et al. (2014); Jaiswal (2013a, b); Romatschke (2012); Yan (2012); Denicol et al. (2012b). We follow closely the derivation of Teaney and Yan Teaney and Yan (2014) so any interested reader can take a look into their work for further details.

It is convenient to rewrite the general RTA-BE in the following suitable form Teaney and Yan (2014)

[TABLE]

where $C_{\bm{p}}=-T^{2}\tau_{r}/(u\cdot p)$ with $\tau_{r}=\theta_{0}/T^{1-\Delta}$. Next we expand the distribution function around $f_{eq.}$

[TABLE]

where $\alpha$ is an arbitrary parameter which keeps track of the order of the gradients. One obtains the following set of coupled equations for $\delta f_{1}$ and $\delta f_{2}$ at first and second order in $\alpha$

[TABLE]

by substituting (124) in the RTA-BE (1) and using the fact that the collisional kernel is $\mathcal{O}(\alpha)$. In order to write the general form of $\delta f_{1}$ and $\delta f_{2}$, it is necessary to classify all the possible irreducible tensors invariant under rotations Teaney and Yan (2014); Denicol et al. (2011); de Groot et al. (1980b). For $d=4$ one gets after some algebra Teaney and Yan (2014); Chattopadhyay et al. (2015); Bhalerao et al. (2014); Jaiswal (2013a, b); Romatschke (2012) 181818Here, we have used the Cauchy-Stokes decomposition of the fluid velocity

$\partial_{\mu}u_{\nu}=-u_{\mu}\left(u\cdot\partial u_{\nu}\right)+\Delta_{\mu\nu}(\partial_{\lambda}u^{\lambda})/3+\frac{\sigma_{\mu\nu}}{2}+\omega_{\mu\nu}\,,$

together with the (conformal) conservation laws to first order in gradients

$\displaystyle D\log T$ $\displaystyle=-\frac{\partial_{\mu}u^{\mu}}{3}\,,$

$\displaystyle Du_{\nu}$ $\displaystyle=-\Delta_{\mu\nu}\partial^{\nu}\log T\,.$

The use of conservation laws at first order in the gradients ensures that the value of $\delta f$ does not change the equilibrium energy density York and Moore (2009).

[TABLE]

where $A_{\langle\mu_{1}}\cdots B_{\mu_{l}\rangle}$ is a projector being traceless, symmetric and orthogonal to the fluid velocity $u^{\mu}$ de Groot et al. (1980b). In these expressions we have defined

[TABLE]

with $f_{eq.}^{\prime}=df_{eq.}(x)/dx$, $\theta=\nabla\cdot u$, $\tilde{p}:=p/T$ and tilde indicates that the quantity is dimensionless. The transport coefficients in $a_{E^{\ast}}$ (128e) will be given below.

For the Bjorken flow the gradient expansion of the distribution function reduces to the expression Blaizot and Yan (2017a, b)

[TABLE]

where $\tilde{\chi}_{p}=-C_{p}f^{\prime}_{eq}$. One immediately recognizes that the Chapman-Enskog expansion of the distribution function is an asymptotic series in terms of $1/(\tau T)$. Close to the thermal equilibrium $1/(\tau T)\approx\tau^{-2/3}$, which in the case of Bjorken flow, it becomes proportional to the Knudsen and inverse Reynolds numbers.

For the Bjorken flow using the second order Chapman-Enskog expansion for $\delta f$ (129) in Eq. (4) results in the asymptotic series expansion of the moments $c_{01}$ and $c_{02}$ as follows

[TABLE]

The transport coefficients appearing in Eq. (130) are given by Teaney and Yan (2014)

[TABLE]

Appendix B A self-contained dictionary for the dynamical systems

In this section, we list a dictionary of words from the terminology of the mathematical field of dynamical systems used in Sect. IV. This is a self-contained and compact dictionary that could help the reader with some terminology of time-dependent dynamical systems.

• 1-

Dynamical system. A set of differential equations in terms of a state vector whose components are real numbers determined by a set of points in some suitable state space. Any small variation in the state of the system leads in turn to a change in the numbers. In short, the evolution of a dynamical system is based on a fixed deterministic rule that describes what future states follow from the current state. When explicitly time-dependent, however, equations of the system also evolve in form, and thus the rule changes as time passes.

• 2-

Non-autonomous dynamical system. A differential equation of the form

[TABLE]

is explicitly time-dependent through the function $g$. This defines a non-autonomous dynamical system if the state function $x(\tau)$ is dynamically evolving with a flow time $\tau$. $T(\tau)$, and equivalent of $g$ in (132) is the r.h.s. of (9).

• 3-

(Asymptotic) fixed point. The time-dependent solutions $x(\tau)$ of $g(x(\tau),\tau)=0$ are the moving fixed points of the system (132). For $\tau\rightarrow\infty$ these define a point(s) $x_{*}$ at which the state vector is in a steady equilibrium state for all times $t>0$ if $\tau=\infty+t$, referred to as fixed points of the non-autonomous system in (132). Note that here the time $\tau$ takes values over $\mathbb{t}=\mathbb{R}^{+}$ 191919In math literature, a non-autonomous system is known to not admit any fixed point by promoting it to an autonomous system of one dimension higher and observing that the r.h.s. of this secondary system

$\displaystyle dx/d\tau=g(x),\quad d\tau/d\tau=1,$

(133) never vanishes. Nonetheless, our definition of a fixed point for a non-autonomous system is not for all times as in the autonomous case, rather over just enough period of time beyond $\tau=\infty$ in order to distinguish it as an equilibrium state. Studying topology of the flow lines around such fixed points is thus more complicated and requires a lot of additional mathematical input Behtash et al. (work in progress).. If we reparametrize flow time to $\rho$ such that $\mathbb{t}=\mathbb{R}$, the fixed points can now be equally defined at $\rho\rightarrow-\infty$. Hence $\tau(\rho)$ turns into a new independent variable in the dynamical system such that

[TABLE]

is autonomous. The fixed points are then defined as usual by finding those $(x_{*},\tau_{*})$ for which $g(x_{*},\tau_{*})=0=h(\tau_{*})$.

• 4-

Flow line (process). A flow line is a map defined by ${\phi}^{\tau,\tau_{0}}({x}_{0}):X\times\mathbb{t}\times\mathbb{t}\rightarrow X$ such that ${\phi}^{\tau,\tau_{0}}({x}_{0})={x}(\tau)$ solves (132) with some initial values $({x}_{0}\equiv{x}(\tau_{0}),\tau_{0})$.

• 5-

Phase portrait. It is a geometric representation of the flow lines for a dynamical system in its phase space given a set of initial conditions represented by $(x_{i},\tau_{i})$ with $i\in I$ where $I$ is simply an index set.

• 6-

Invariant manifold. As defined above, it is a subset of the phase portrait which contains all the flow lines initiated inside it (or at its boundary) which keep being restricted in there at all times.

• 7-

Complete flow line. Choosing a different parametrization of flow time $\tau$, say $w$, may allow us to connect two fixed points at $(x_{0},0)$ and $(x_{\infty},\infty)$. Note again that the time manifold $\mathbb{w}$ now allows 0 to be a valid input. This flow line if existed, would be called complete. The set of fixed points of the dynamical system all belong to the boundary of the invariant manifold iff there is one complete flow line between every two fixed points.

• 8-

Stability of a fixed point. Linearizing the system in (134) about a fixed point $(x_{*},\infty)$ leads to the understanding of how flow lines behave around that fixed point. Eigenvalues of the linearization matrix $J$ with components

[TABLE]

evaluated at $(x_{*},\infty)$, determine the stability of $(x_{*},\infty)$. Let us denote the eigenvalues of (135) by ${\lambda}_{\infty}$. Throughout this section we only consider hyperbolic fixed points for which $\Re(\lambda_{\infty})\neq 0$.

• 9-

The Hartman-Grobman (HG) theorem. It states that near a hyperbolic fixed point $(x_{*},\infty)$, the nonlinear system (132) has the same qualitative structure as the linear system

[TABLE]

Therefore, HG implies that two hyperbolic systems are locally topologically flow equivalent iff their unstable manifolds have equal dimensions.

• 10-

Index of $(x_{*},\infty)$. The number of components of ${\lambda}_{\infty}$ whose real parts are positive is called the (Morse) index of that fixed point, which is technically the dimension of the unstable manifold around every hyperbolic fixed point. So any nonzero index means that fixed point is basically unstable. This instability is severe for a larger index. For example, in a $d$ dimensional system $(x_{1},\dots,x_{d-1},\tau)$, an index $n$ fixed point is an unstable saddle whose stable manifold is $d-n$-dimensional. The instability of an index $n$ fixed point known as a source is far more pronounced than a saddle. Finally, a stable fixed point (so-called sink) has always a vanishing index 202020The imaginary part of ${\lambda}_{\infty}$ if present, means that the flow line is going to spiral its way in or out; thus the name stable or unstable spiral, respectively. .

• •

Invariant manifold. Suppose $\phi$ is a flow on the phase space $(X,\mathbb{t})$ equipped with a spatial metric $\rm dist$. A family ${\cal A}:=\bigcup_{\tau\in\mathbb{t}}A_{\tau}$ of nonempty subsets of $X$ is invariant with respect to $\phi$ if for all $\tau\geq\tau_{0}$

[TABLE]

${\cal A}$ is called the invariant manifold.

• 12(a)-

Forward attraction. Let $\phi$ be a flow (process). A nonempty, compact and invariant manifold ${\cal A}:=\bigcup_{\tau}A_{\tau}$ for compact $I$ is said to forward attract if

[TABLE]

for all $x_{0}\in X$ and $\tau_{0}\in\mathbb{t}$. The distance function between two sets $A,B$ is defined as usual on the metric space $X$ as ${\rm dist}(A,B)=\inf_{x\in A,y\in B}{\rm dist}(x,y)$.

• 12(b)-

Pullback attraction. A nonempty, compact and invariant non-autonomous set ${\cal A}:=\bigcup_{\tau}A_{\tau}$ is said to pullback attract if

[TABLE]

for all $x_{0}\in X$ and $\tau\in\mathbb{t}$.

• 13-

Basin of attraction. A neighborhood of an attractor is called the basin of attraction $B$ if it consists of all points that evolve to the attractor in the limit $\tau\rightarrow\infty$. In other words, $B$ is the set of all points $b\in B$ in the phase space provided that for any open neighborhood $N$ of the attractor, there is a positive constant $s$ such that $\phi^{\tau,\tau_{0}}(b)\in N$ for all real $\tau>s$.

• 14-

Lyapunov function. Let us assume that the point $(x,\tau)$ at $\tau\rightarrow\infty$ is $(0,\infty)$ which defines an asymptotically stable fixed point of a non-autonomous dynamical system. Then if there is a function $V(x,\tau)$ such that the Lyapunov stability conditions

[TABLE]

are satisfied, $V(x,\tau)$ is called a Lyapunov function.

It may be appropriate here to mention that in general, 11(a)-(b) are very different and independent concepts. Pullback attractor (139) contains information about the past (early-time) of a non-autonomous dynamical system, whereas forward attractor (138) makes use of information about the future. It well may be the case that one exists while the other one does not in a non-autonomous system. However, in the case of an autonomous system, they are equivalent and the resulting attraction is global with the set $\mathcal{A}$ referred to as global attractor.

Appendix C Constructing the transseries for $T(\tau)$

C.1 Algebraic properties of transseries

Let us suppose the following transseries is given

[TABLE]

where $\phi^{({\bf n})}_{k}$ are complex-valued coefficients. We define the basis of the transseries aka $E^{({\bf n})}_{k}(z)$ as

[TABLE]

The transseries can be regarded as an element of a vector space defined as

[TABLE]

The operations on this vector space are briefly explained below. Note that the transseries is always closed under these operations.

Product

The product operation $\times$ between two transseries is an additive operation on both transmonomial and asymptotic orders given by

[TABLE]

Exponential

The exponential operation ${\rm Exp}$ is an extension of the product to the exponential terms defined by

[TABLE]

where $\Phi^{({\bf n})}_{k}\in{\cal H}^{({\bf n})}_{k}$. Thus, the exponential is a linear mapping of the form

[TABLE]

Integration

We define an integration operation of the form

[TABLE]

where $\Phi^{({\bf n})}_{k}\in{\cal H}^{({\bf n})}_{k}$ and $s\in{\mathbb{N}}_{0}$. For realizing the closedness of transseries under the integration, we have to assume that

[TABLE]

Then after some algebra, Eq. (155) casts ${\rm Int}_{s}$ in the form of the following linear mapping

[TABLE]

where $R_{l}(a)$ is written in terms of the Pochhammer symbol $(a)_{l}$,

[TABLE]

We recall that to get the expression in (158), use has been made of the asymptotic form of the upper incomplete Gamma function for $z\rightarrow\infty$, namely

[TABLE]

C.2 Transseries solution for the temperature $T(\tau)$

As we have seen in the previous section, $\hat{T}(\hat{\tau})$ can be expressed by the same transseries ansatz used for $c_{l}(\hat{\tau})$ due to the fact that $c_{l}(\hat{\tau})=\mathcal{O}(\hat{\tau}^{-1})$. Using Exp and ${\rm Int}_{s}$ mappings defined in Eq. (154) and Eq. (158) respectively, we find that

[TABLE]

By assuming the ansatz for $c_{01}(\hat{\tau})$ as

[TABLE]

we can get from Eq. (163)

[TABLE]

Hence, it is easy to read off the coefficients

[TABLE]

Therefore, inserting Eq. (167) in (162) results in

[TABLE]

where $\Psi(\hat{\tau})$ converges to zero at late times, that is $\lim_{\hat{\tau}\rightarrow\infty}\Psi(\hat{\tau})=0$.

C.3 The ${\bf c}$-$T$ dynamical system

In the analysis of the previous section, we used $w$ as a time variable, but something else that may be done is to consider the same analysis for the dynamical system including explicitly the temperature, which will require constructing transseries for $T(\tau)$ as well. In this subsection, we seek to prove that this dynamical system entails the same structure as (11), which is compatible with Costin’s prepared form.

The dynamical system in the $\tau$-coordinate takes the following form

[TABLE]

where ${\mathfrak{B}}$ is a constant matrix, ${\mathfrak{D}}$ is a diagonal constant matrix, and ${\bf A}$ is a constant vector. From Eq. (172), one can obtain the formal solution of $T(\tau)$ as

[TABLE]

where $C_{0}$ is an integration constant, and $C_{T}(\tau)$ is defined as

[TABLE]

In the IR, or the late-time limit $\tau\rightarrow\infty$, $T(\tau)$ is observed to behave like

[TABLE]

with a constant $T_{0}$. Hence, by redefining this constant, the temperature $T(\tau)$ takes the form

[TABLE]

Now, by working (176) into (171) and changing the parameters $\hat{\tau}$ and $\hat{\theta}_{0}$ as

[TABLE]

the dynamical system becomes

[TABLE]

where

[TABLE]

Note that the temperature $T(\hat{\tau})$ can be reproduced by

[TABLE]

As discussed in App. C, we find that $\hat{T}(\hat{\tau})$ will be expressed by

[TABLE]

where $\Psi(\hat{\tau})$ is a continuous function with a transseries formula identical to that of $\tilde{\bf c}$, which converges to zero asymptotically. Therefore, Eq. (178) bears the familiar prepared form as in Eq. (11) with only one significant difference from a physical standpoint: the temperature now is a dimension in the phase space of the dynamical system. In the UV limit, the divergence of $T(\tau)$ as $\tau\rightarrow 0$ in general (see Eq. (172)) suggests that $T=0$ and $T\neq 0$ theories are disconnected from each other, and there cannot be any flow line (process) that connects a UV fixed point to the hydrodynamic fixed point achieved at $\tau\rightarrow\infty$ in the original kinetic model of the Bjorken flow.

Appendix D Bjorken dynamical system and its trasseries solutions for general $\tau_{r}=\theta_{0}/T^{1-\Delta}$

In this appendix, we extend the dynamical system in Eq. (9) to include general relaxation time, $\tau_{r}=\theta_{0}/T^{1-\Delta}$ with $0\leq\Delta\leq 1$, which is very straightforward. The generalized system has the same number of asymptotic fixed points with their stability being identical to the case $\Delta=0$ as discussed in Sect. IV.1. So except for some small quantitative alterations in the flow lines around UV and IR, the qualitative picture of the phase portrait remains the same. In other words, there are no changes in global dynamics.

The generalized dynamical system is described by

[TABLE]

As in the main bulk of the paper, we change to the coordinate $w=\tau T^{1-\Delta}$ and consider the dynamical system to be of the form

[TABLE]

where the matrices and vectors involved here have the same definitions as in (41). Being again compatible with Costin’s prepared form, we conclude that the transseries can be explained with the same form used to explore the solutions of the $\Delta=0$ model defined by (15). By defining an invertible matrix $U$ for the diagonalization purposes again, one finds the recursive relation involving the transseries data as

[TABLE]

with the IR data $S_{i}$ and $\tilde{b}_{i}$ in the transseries being

[TABLE]

where $b_{i}$ are the eigenvalues of the matrix ${\mathfrak{B}}$.

In addition, the transaymptotic mathing is obtained in the following, where $\hat{\zeta}_{i}=\partial/\partial\log{\zeta}_{i}$. Summing over ${\bf m}$ yields the transasymptotic matching for the general $\Delta$ as

[TABLE]

The transseries and leading order asymptotics of five moments $c_{01},\dots,c_{21}$ for $\Delta=1$ system are depicted in Fig. 10. To avoid repeating ourselves here, we refer the interested reader to the explanations given in the text prior to Fig. 1.

D.1 Transasymptotic matching for $c_{01}$ of the $N=0,L=1,\Delta=1$ dynamical system

As a result of Eq. (187), the mode-to-mode coupling term $c_{01}\frac{dc_{nl}}{dw}$ vanishes, and therefore the transasymptotic matching effectively becomes a finite expansion of the transmonomials packaged into $\zeta$ at each order. The first five transseries coefficients are then found to be exactly given by

[TABLE]

Appendix E Review of the exact solution of the RTA-BE

In this section, we briefly touch down on some aspects of the exact RTA-BE that have been taken advantage of throughout this paper. For a more complete and detailed explanation of this solution, we refer the reader to Refs. Baym (1984); Florkowski et al. (2013b, a).

The exact solution of the RTA-BE (1) is

[TABLE]

where $D(\tau_{2},\tau_{1})=\exp\left[-\int_{\tau_{1}}^{\tau_{2}}d\tau^{\prime}/\tau_{r}(\tau^{\prime})\,\right]$ and $f_{0}\left(\tau_{0},p_{T},p_{\varsigma}\right)$ is the initial distribution function at $\tau=\tau_{0}$. Here, we consider the RS distribution function Romatschke and Strickland (2003)

[TABLE]

where $\Lambda_{0}$ and $\xi_{0}$ are the initial temperature and initial momentum anisotropy along the $\varsigma$ direction, respectively. With (194) at hand, the energy-momentum conservation law leads to the following integral equation Baym (1984); Florkowski et al. (2013a, b):

[TABLE]

The functions $\mathcal{R}_{nqr}$ appearing in the previous expression are given by Molnár et al. (2016)

[TABLE]

Eq. (195) is solved with the help of the iterative methods discussed in Refs. Banerjee et al. (1989); Florkowski et al. (2013a, b); Denicol et al. (2014b). By plugging the exact solution (193) in the definition of the moments $c_{nl}$ in (4), one finds

[TABLE]

where

[TABLE]

In Eq. (197) the reader must bear in mind that the Landau matching condition for the energy density at $\tau=\tau_{0}$ when using the RS distribution function (194) implies $T_{0}=\mathcal{R}^{1/4}(\xi_{0})\Lambda_{0}$, where the function $\mathcal{R}(\xi)$ is defined by Martinez and Strickland (2010)

[TABLE]

The moments $\bar{M}^{nl}$ (114) can now be evaluated from the exact solution (193) as follows

[TABLE]

where

[TABLE]

The exact solutions for both $c_{nl}$ and $\bar{M}_{nl}$ are determined numerically by having the expression for the temperature from Eq. (195).

Appendix F Hydrodynamization of soft and hard modes for $\Delta=1$ system

In Sect. V we analyzed the hydrodynamization process of the normalized moments $\bar{M}^{nm}$ for the conformal system, i.e., $\tau_{r}=\theta_{0}/T$. Here, we give the results of hydrodynamization when $\tau_{r}=\theta_{0}/T^{1-\Delta}$ by fixing $\Delta=1$. For this choice, the variable $w=\tau(fm/c)$ becomes dimensionful. The soft and semi-hard limits of the normalized moments are given by the expressions in Eqs. (118) and (119), respectively. Also, the exact moments $\bar{M}^{nm}$ are given by Eq. (200).

The results of the hydrodynamization of the normalized moments are presented in Figs. 11-14. From these plots, we arrive at the conclusions already discussed in Sect V: the larger $n$ and $m$, the longer it takes for $\bar{M}^{nm}$ ($n>0,l\geq 1$) to thermalize while $\bar{M}^{n0}$ ($n\geq 0$). Furthermore, the IR behavior is uniquely determined by both non-hydrodynamic modes $c_{01}$ and $c_{11}$ for $\bar{M}^{nm}$ ($n>0,l\geq 1$) while the non-hydrodynamic modes $c_{10}$ and $c_{20}$ for $\bar{M}^{nm}$ ($n\geq 0$).

When calculating the deviation from equilibrium by following the same procedure discussed in Sect. V we reach the same conclusions presented in there: the soft regime does not capture the late-time asymptotics correctly while there is an extremely good match between the soft+semi-hard regime and exact results. This non-trivial double check gives a solid evidence of the importance of considering not only the effective shear mode $c_{01}$ but also the new non-hydrodynamic mode $c_{11}$.

Appendix G Exact solution vs. truncation of the dynamical system

In the numerical analysis presented in this work, we used various truncations in the number of moments in the original dynamical system (9b). Both the exact solution to RTA-BE and the solution to the truncated dynamical system share the same IR fixed point but as discussed in Sect. IV, the solution of the truncated dynamical system cannot fully agree with the exact solution at early times due to the missing corrections due to the truncation. For simplicity we take $\Delta=0$ in Eq. (2).

In this section, we test how good the truncation schemes would match the exact RTA-BE result. We compare the moments $c_{nl}$ obtained from Eq. (200) by solving RTA-BE exactly against solutions to the truncated dynamical system in (22b), where the truncation is controlled by the number of involved moments. For this reason, we solve numerically the dynamical system in (9b) by varying its dimension from $1$ to $L$ and compare against exact result (4) by using the same set of initial values of the moments $c_{l}$. The initial conditions for the exact $c_{nl}$ (198) are chosen to be $\tau_{0}=0.25~{}\mathrm{fm/c},T_{0}=0.6~{}\mathrm{GeV}$ and $\xi_{0}=1000$. In Fig. 16 we present the numerical results for our findings by considering the initial conditions $c_{0l}(w_{0})=c_{0l}(w_{0};\xi_{0})$. When the dimension of the truncated dynamical system increases, they get closer to the exact solutions of the full RTA-BE. We notice that the larger the index $l$, the more moments are needed in order to match the exact solution of RTA-BE.

In Fig. 17, the exact moments $c_{nl}\quad(n=1,3,5)$ are compared to the solutions of the truncated system. This figure encodes the information of the non-linear structure of the original dynamical system: the evolution equation of $c_{nl}$ couples only to the moments of the same or lower order in $n$, and at the same time it couples to the next moment $c_{n,l+1}$. So $c_{nl}$ in Fig. 17 become closer to the exact RTA result if $L$ is increased for a fixed value of $n$. We comment that the numerical results presented in Fig. 16 and 17 do not depend on the particular choice of initial conditions as long as we are bound to the basin of attraction of the invariant manifold of the truncated system.