Exponential Approach to the Hydrodynamic Attractor in Yang-Mills Kinetic Theory

We use principal component analysis to study the hydrodynamic attractor in Yang-Mills kinetic theory undergoing the Bjorken expansion with Color Glass Condensate initial conditions. The late time hydrodynamic attractor is characterized by a single principal component determining the overall energy scale. How it is reached is governed by the disappearance of single subleading principal component characterizing deviations of the pressure anisotropy, the screening mass and the scattering rate. We find that for wide range of couplings the approach to the hydrodynamic attractor at late times is well described by an exponential. Its decay rate dependence on the coupling turns out to translate into a simple dependence on the shear viscosity to entropy density ratio.


I. INTRODUCTION
Ultrarelativistic heavy-ion collisions at RHIC and LHC produce a collective state of matter comprised of the fundamental quark and gluon constituents of strong interactions. Describing ab initio formation and evolution of this quark-gluon plasma (QGP) presents an important theoretical problem. Since a first principles description directly using the theory of strong force, quantum chromodynamics (QCD), remains an outstanding challenge, the standard model of the space-time evolution of heavyion collisions is based on effective descriptions of QCD [1]. The initial non-equilibrium dynamics after the collision has been addressed in a variety of different microscopic models [2,3]. Eventually, on a time scale of ∼ 1 fm/c, the QGP can be well described by relativistic hydrodynamics. This is to some extent surprising, as the system can feature large spatial and temporal gradients and remains significantly out of equilibrium for a much longer period of time, as indicated e.g. by large values of the pressure anisotropy.
Due to these concerns the applicability of relativistic hydrodynamics in high-energy collisions of heavy and light nuclei has repeatedly been questioned [4], and different proposals have been put forward to extend the applicability of hydrodynamics to earlier times and more anisotropic systems [5]. One proposal to potentially extend the validity of hydrodynamics goes by the name of hydrodynamic attractors [6], which represent emergent constitutive relations away from equilibrium. Over the course of the last few years, the emergence of such hydrodynamics attractors has been firmly established in a variety of effectively 0+1 dimensional microscopic models which feature a rather high degree of symmetry [7][8][9][10][11][12][13][14], and different theoretical [15][16][17][18][19][20][21][22][23] and phenomenological aspects of the emergence of hydrodynamic attractors [24][25][26][27] have been explored. See [28] for a recent review.
Despite efforts in this direction [15,29], the generalization of the concept of hydrodynamic attractors to higher dimensional systems represents an outstanding challenge. With this as a motivation, it was recently pointed out that the hydrodynamic attractor can be viewed as a dimensionality reduction in a space of natural observables associated with nuclear collision problems [13]. From this perspective, the number of relevant degrees of freedom required to describe the evolution of a system is effectively lowered due to a rapid memory loss of initial conditions. Viewing thermalization in this way enables to approach the study of thermalization in a data-driven way and borrow methods from data science.
In the present work we adopt the approach of [13] to study the evolution of a 0+1D Bjorken flow in an effective kinetic theory (EKT) of pure glue QCD [24,25] with highly anisotropic color glass condensate [30,31] initial conditions. Building on the methodology developed in [13], we use principal component analysis (PCA) to study the dimensionality reduction for a set of observables that describe the thermalization process at weak coupling. The mechanism behind the information loss can come from the expansion at early time or interactions at late times [19,32]. We also look more closely at the late time behaviour to gain information on transient nonhydrodynamic contributions to the pressure anisotropy, and compare our results in pure glue QCD kinetic theory to previous studies in conformal Relaxation Time Approximation (RTA) [33,34].
The paper is organized in the following way: In Sec. II, we introduce our model of the initial state, along with the effective kinetic theory of pure glue QCD describing the non-equilibrium dynamics of the QGP. Subsequently, in Sec. III, we perform a PCA of the evolution for different initial conditions to study information loss and dimensionality reduction over the course of the thermalization process of the pre-equilibrium QGP. Sec. IV is devoted to a more detailed study of the evolution towards equilibrium, focusing on the exponential decay of variance of the pressure anisotropy in the QGP. We conclude this paper in Sec. V with a short summary of our most important findings and an outlook on potential future works.

II. SETUP
We use the Color Glass Condensate (CGC) effective theory of high-energy QCD [30,31], to evaluate the phase-space distributions of gluons produced in the collision. Due to physical differences and uncertainties in modeling the initial state, this gives rise to a multidimensional parameter space of early time phase-space distributions of gluons. They serve as initial conditions for EKT, which we subsequently utilize to evolve expanding nuclear matter into the hydrodynamic phase.

A. Color Glass Condensate initial conditions
We will compute the initial spectrum of gluons based on the k T -factorization formula [35,36] where dNg d 2 bd 2 Pdy describes the transverse momentum (P) spectrum of gluons produced per unit rapidity (y) and transverse area (b). By N c = 3 we denote the number of colors, g is the Yang-Mills coupling, b 0 denotes the impact parameter of the nucleus-nucleus collision and φ A/B (b, k) is the un-integrated gluon distribution from each nucleus (A or B). We employ a particularly simple parametrization of the nuclear gluon distributions, due to Golec-Biernat and Wusthoff (GBW) [37], for which where Q 2 A,B = Q 2 s (b) denotes the (adjoint) saturation scale [38] for nucleus.
Based on Eqns. (1) and (2), the initial gluon spectrum in the GBW model can then be expressed analytically as denote the root mean square average and the ratio of the saturation scales. Since the saturation scales Q 2 A/B are generically proportional to local density of nuclear matter, we can anticipate that for different collision geometries the ratio x = Q A /Q B will vary, and we will consider variations 1 ≤ x ≤ 6 in the following.
Since the transverse spectrum in Eq. (3) features a ∝ 1/P 2 infrared behavior, soft observables, such as screening mass m D and scattering rate g 2 T * in Eq. (14), are infrared divergent. However, this divergence can be regulated by non-linear effects, going beyond the factorization formula, see e.g. the discussion in [35]. Within our study, we will simply model this by introducing an additional infrared regulator, by virtue of the replacement 1 P 2 → P 2 (P 2 +m 2 ) 2 in the first factor of Eq. (3), which ensures that m D and g 2 T * in Eq. (14) remain finite, while the initial energy density is not affected by the regulator as long as max(Q 2 A , Q 2 B ) m 2 . Since corrections to the factorization formula become important for P min(Q A , Q B ), in the following we choose the infrared regulator as µ ≡ m . We further note that due to the particular simplicity of the GBW model, the initial spectrum in Eq. (3) exhibits and exponential decay at high momentum P Q, whereas a more realistic parametrization should give rise to a Q 4 /P 4 power-law tail [39]. However, since the energy density of the system is dominated by momenta P ∼ Q, we believe that the model is adequate to study the thermalization of the bulk QGP, while neglecting the impact of high-energy degrees of freedom.
Based on the transverse momentum spectrum, the initial phase-space distribution f (x, p) can then be obtained as [40] where η denotes the space-time rapidity and we employ these as initial conditions for the subsequent kinetic description at an initial proper time τ 0 = 1/Q, where a quasi-particle description first becomes applicable [2,3]. While in the high-energy boost-invariant limit, the initial distribution is proportional to δ(y − η), any interactions will immediately broaden the longitudinal momentum distribution and we therefore consider a smearing form of the delta function in the initial conditions. We treat the width of longitudinal rapidity distribution σ as the second free parameter in the initial conditions and consider variations in the range 0.05 < σ < 0.50 in the following. Starting from the phase-space distribution, the energymomentum tensor T µν is defined as a its second moment where ν g = 2(N 2 c − 1) denotes the degeneracy factor for gluons. Specifically, one may evaluate the energy density , the transverse and the longitudinal pressure P T /L in Milne coordinates as = T τ τ = d 2 P (2π) 2 dp (2π) p ν g f g (τ, b, P, p ), (7a) where we have re-expressed the phase-space distribution f g in terms of the momentum variables p = p τ = P 2 + p 2 and p = τ p η = |P|sinh(y − η). They denote the total and the longitudinal momentum in the co-moving frame. When considering variations of the parameters x and σ, we adjust the value of the average saturation scale Q, to keep the initial transverse pressure per unit rapidity g 2 τ 0 P T (which is independent of the coupling g) constant in physical units. Hence the parameter space of the initial conditions is two dimensional and spanned by the ratio of the nuclear saturation scales x and the longitudinal smearing width σ.

B. Effective kinetic theory of pure glue QCD
We model the evolution of excited nuclear matter using kinetic description. We follow previous works [12,41,42] and treat the system as longitudinally boost-invariant and locally homogeneous in the transverse plane. The dynamics reduces to a (0+1)-dimensional problem and is governed by an effective kinetic for pure glue QCD [43] denotes the leading order elastic collision integral for gluons and C 1↔2 g [f ] is the inelastic collision integral that describes emission/absorption of gluon radiation We note that the matrix element |M gg→gg (p, p 2 |p 3 , p 4 )| 2 for elastic scattering is self-consistently screened following the prescription of [44], and that the effective inelastic rates dΓ g gg dz (p, z) account for the Landau-Pomeranchuk-Migdal effect [45][46][47] via an effective vertex resummation [43]. Details of the algorithms and numerical implementations can be found in [48]. Early initialization requires a fine discretization of the momentum space variables for the EKT simulations, and if not stated otherwise we employ N p = 1024 and N cos(θ) = 512 to discretize p and cos θ = p p in the range between p min /(g 2 τ 0 P T ) 1/3 = 0.01 and p max /(g 2 τ 0 P T ) 1/3 = 20.
Finally, note that in pure glue QCD kinetic theory the interaction strength g and the number of colors N c enter the Boltzmann together as the 't Hooft coupling λ [49] λ ≡ g 2 N c .
Therefore, in the following we will express the interaction strength in terms of λ.

A. Hydrodynamization & Emergence of attractors
One prominent indicator of transition to hydrodynamics is the longitudinal pressure over what would be the equilibrium pressure ratio P L /( /3) or any of its closely related variants. We present it in Fig. 1 as a function of the universal time scale [14,[50][51][52] where η is the shear viscosity and P = /3 is the thermodynamic pressure for a conformal system. It can be thought of as measuring the physical proper time in units of an effective equilibrium relaxation time where T is an effective temperature of non-equilibrium system determined by = ν g π 2 T 4 30 and s is the entropy density. While the overall scaling with T follows on dimensional grounds for a conformal system, the normalization factor of 1/4π is conventional and inspired by the strong coupling result η/s = 1/4π of [53,54].
Different colored curves in Fig. 1 show to the results obtained for different ratios of the saturation scales x and different longitudinal smearing width σ. Different dash styles correspond to the evolution for three different coupling strengths λ = 5, 10, 20, for which the relevant viscosities η/s ∼ O(1) as summarized in Tab. I in the Appendix.
Starting from the initial conditions, one observes a significant variation of P L /( /3), where small smearing parameters σ feature a highly anisotropic initial distribution, with almost vanishing longitudinal pressure. Since in all cases the initial time τ 0 = 1/Q is kept fixed, the curves for different simulation parameters x, σ, λ start at different values ofw; in this way larger saturation scale ratios x also tends to provide relevantly more anisotropic initial distribution. During the initial stages the evolution is dominated by the longitudinal expansion of the system [19,32], such that the ratio of P L /( /3) undergoes a quick memory loss and at intermediate times approaches a small value irrespective of the initial conditions. Eventually, on a time scalew 1, all the different curves merge towards the hydrodynamic limit, which to first order in gradients can be universally expressed as [52] .
While the emergence of a pressure attractor at 0.1 <w < 1 clearly indicates an early reduction of the variance due to the initial free-streaming expansion dynamics, it is equally important to point out that the full convergence to a common attractor only occurs at later timesw > 1, when the system relaxes towards hydrodynamics.

B. Principal component analysis
Next, in order to better understand the emergence and properties of the attractor in pure glue QCD kinetic theory, we will scrutinize the evolution further by including additional observables, which go beyond probing the bulk anisotropy of the system. Beyond the anisotropic pressure P L and P T in Eqns. (15), we will study the screening mass m 2 D and the collision rate Since m 2 D and g 2 T * govern the strength of elastic and inelastic interactions in pure glue QCD kinetic theory, they present the most natural quantities to include in a weak coupling analysis.
Since all of of the above quantities P L , P T , m 2 D , T * are dimensionful, it is natural to consider dimensionless ratios along with one single quantity that is sensitive to the overall energy scale. Specifically, for our analysis, we choose to consider P L /( /3) and normalize m 2 D , T * by the transverse pressure P T . This gives us the following dimensionless quantities All of them by construction approach unity in thermal equilibrium. While initially, the transverse pressure P T is kept constant for different initial conditions, differences in the kinetic evolution gives rise to variations of the overall energy scale, which we monitor in terms of the dimensionless quantitȳ where we employ C ∞ = 0.98 [25]. By following the arguments of [14,25], the denominator in Eq. (15d) provides an estimate for the late time asymptotic value of τ 4/3 P T , such that for typical initial conditions this quantity can again be expected to be close to unity at late times. We can directly verify these expectations in Fig. 2, where we present the evolution of the observables P T ,P L ,m 2 D ,T * as a function of the universal time scalẽ w for a variety of different initial conditions at three different coupling strengths λ = 5, 10, 20. While the selfnormalized quantitiesP L ,m 2 D ,T * all converge to a common attractor behavior at late times, the overall energy scaleP T develops a ∼ 15% variation over the course of the non-equilibirum evolution of the system.
With the set of observables in Eq. (15), we can further perform a PCA to extract the most significant contributions in the emergence of the attractor phenomenon, namely the principal components of the attractor. For this purpose, rather than studying these quantities separately, we can also consider these quantities as parametrizing a four dimensional space. Each solution is represented as a vector evolving in time. Evidently, to perform a meaningful comparison of the different directions in X(w) space, the units and overall scales of the different components of X(w) must be comparable, which is precisely the reason that the observables have been normalized as in Eq. (15). From this point of view, thermalization is achieved through dimensionality reduction. With sufficient variation of the initial phase-space distribution function, the set of initial conditions {X 0 } collectively spans some volume in the four dimensional space. However, since the initial conditions used here only depend on two parameters, we only expect to span (at most) a two dimensional subspace. Over the course of the thermalization process, the dimensionless variablesP L ,m 2 D ,T * all approach unity, and onlyP T which is sensitive to the overall energy scale should have a significant variation at late times, implying that ultimately the solutions will span a one dimensional subspace.
PCA is a simple method to study this process of dimensionality reduction. Given a set of points {X}, the PCA calculates the eigenvectors and eigenvalues of the covariance matrix C mn = Cov(X m , X n ). (17) to produce a set of orthonormal vectors such that the first vector points along the direction of the largest variance of the data set, and the subleading vectors are similarly optimized in the space orthogonal to the leading vector. Each vector is associated with an explained variance, quantified by the variance of the set of points in its direction. The number of non-negligible explained variances gives a measure of the dimensionality of the region occupied by the states of interest at a given value of time variable. We illustrate this behavior in Fig. 3, where we show the distribution of values in theP L ,P T plane at three different timesw = 0.2, 0.9, 1.3 of the evolution. Each point in Fig. 3 corresponds to the evolution for a particular initial condition, while the blue and orange arrows indicate the first and second principle components. While at early times (w = 0.2) the different initial conditions cover a two-dimensional subspace, with the largest variations in theP L direction, the effective reduction of the dimensionality of the distribution is clearly visible, as at late times w = 1.3 all points converge towards a one dimensional manifold oriented along theP T direction.
By applying PCA to the set of solutions at each point in the universal timew, we can study the dimensionality of the dataset through the explained variances, and identify which directions in the space of observables that show the greatest variance. Our results are compactly summarized in Figs The decomposition of the first and second principal components are shown in Fig. 5. The abrupt shift in behaviour atw ≈ 0.5 is correlated with the crossing of explained variances in Fig. 4, and should be thought of as the two vectors switching identity. At late times the first component is dominated byP T , which was chosen to be sensitive to the energy scale of the state. The second component is dominated byP L andm 2 D . While Fig. 2 already shows that the observablesm 2 D ,T * also features a similar attractor behaviour as the well-knownP L , it does not tell us whether these attractors are actually the same. In contrast, the principal component analysis in Figs. 4 and 5 reveals that the evolution of the different observables is highly correlated and can be captured with a single principal component.

IV. EVOLUTION TOWARDS EQUILIBRIUM
So far we have statistically analyzed the emergence of hydrodynamics attractors starting from CGC initial conditions at very early timesw 1. While in this case the dynamics is initially dominated by the longitudinal expansion, we performed additional simulations that focus on the late time approach to hydrodynamic behavior to further analyze the memory loss and convergence towards hydrodynamic behavior in pure glue QCD kinetic theory. By initializing the system according to a Romatschke-Strickland type distribution [55], for different initial anisotropies ξ 0 = 1.25, 2.5, 5, 10 at initial timew 0 = 0.1, 0.3 and ξ 0 = 1.1, 1.2, 1.3 at initial timẽ w 0 = 1, 3, 10, we are then able to analyze the effective memory loss at late times for a large range of couplings λ = 0.1−20. We note that in all cases, the normalization factor c(ξ 0 ) = 2 is chosen such that the initial energy density (w 0 ) = π 2 30 ν g T 4 0 remains the same irrespective of the initial anisotropy ξ 0 ; while in order to initialize the simulations at the samew 0 , the initial proper time is adjusted as τ 0 = 4πη/s T0w 0 in Eq. (18) for the respective coupling strength. For completeness, our numerical implementation employs the following discretization: N p = 64 and N cos(θ) = 64 with p min /T 0 = 0.01 and p max /T 0 = 8.

A. Exponential approach to hydrodynamics at late times
In order to investigate direct approach to viscous hydrodynamics at late times, we will focus on the evolution of P L /( /3) inw. As we discussed in Sec. III B, this quantity has a significant contribution to the leading principal component characterising deviations from the attractor, and exhibits a universal late time hydrodynamic behavior determined by Eq. (13).
The evolution of the P L /( /3) for the initial conditions (18) is depicted in Fig. 6. Different colored curves correspond to the results for different coupling strength λ = 0.1 − 20, while curves of the same color correspond to different initial conditions ξ 0 at different initialization timesw 0 . Irrespective of the coupling strength, one observes that very quickly after the initialization, different initial conditions appear to converge towards a common attractor curve. By careful inspection, one notes that for varying coupling strength slight differences in the attractors persist at intermediate timesw 3, before eventually all curves converge towards the same (by construction) hydrodynamic late time behavior (13). Note that contributions from the terms second and higher order in derivatives are expected to exhibit residual coupling dependence, which might at least partially explain the slight differences between attractors at earlier times.
In order to further characterize the approach towards an attractor, we compute the variances of P L /( /3) for the different initial conditions ξ 0 at each value of the initialization timew 0 and coupling strength λ. Before we discuss our results in pure glue QCD kinetic theory, shown in the middle panel of Fig. 6, it proves insightful to recall the behavior previously observed in different microscopic models of early-time dynamics of the QGP.
Starting with [6] (see also [56,57] for a more complete discussion), it was understood that in a class of models employing the Müller-Israel-Stewart (MIS) approach to embed hydrodynamics in a framework compatible with relativistic causality [58][59][60], the late time behavior of P L /( /3) can be elevated into a transseries [61] of the form Without dwelling into details of this formal expression, the key aspect for us are exponentially suppressed in times effects with decay rates associated with Ω j . In MIS different Ω j are just integer multiples (due to nonlinear effects) of a single relaxation scale present in this class of theories. The same structure was verified to appear in RTA kinetic theory [33,34,62] and in holography [62][63][64][65].
Most importantly, the contributions Ω i describe the decay of transient non-hydrodynamic contributions to the pressure anisotropy, which in the aforementioned examples can be related to the analytic structure of retarded correlation functions of the energy-momentum tensor in equilibrium [5,66]. In the conformal RTA kinetic theory the shear viscosity and the relaxation time are related by .
Following [67], the relaxation time determines the decay rate of non-hydrodynamic contributions in the boostinvariant background as where we have used the fact that T (τ ) scales as τ −1/3 at sufficiently late times in Bjorken flow and dropped subleading contributions at late time. While Eq. (22) provides the rate of decay of non-hydrodynamic contributions in each individual realization, the variance measures the square of these contributions and thus decays twice as fast.
While the structure of non-hydrodynamic excitations in QCD kinetic theory is generally expected to be rather complicated [68,69], it is nevertheless interesting to investigate to what extent a simple parametrization of the type in Eq. (20) can describe the convergence towards an attractor.
To this end, inspection of the evolution of the variances in the middle plot in Fig. 6 shows a clear exponential decay of the variance for most values of λ andw 0 . Clear deviations from an exponential decay are only seen for data initialized at very early times, particularly for small coupling, where subleading contributions to the pressure anisotropy can be expected to be more important.
We extract the effective decay rate Ω eff of the variance for each choice of λ andw 0 by a linear fit in a variety of windows. Extracted values of the decay rate are presented in the bottom plot of Fig. 6, along with the mean value of those fits and for comparison we also show the decay rate 12 5 π in conformal RTA kinetic theory. While not perfectly exponential, the data initialized at early times have a smaller effective decay rate at smaller couplings. For the data initialized at late times, which show a clear exponential decay, the decay rate is approximately independent of the coupling strength λ (with variations within 5% as the coupling varies by over two orders of magnitude) and turns out to be rather close to the value in RTA. While our finding suggests that, just like in RTA, the decay rate in QCD kinetic theory appears to be closely connected to the value of viscosity, we certainly do not have a good explanation of this behavior. We note however, that similar observations were made also at the level of second order transport coefficients [70]. In particular, this implies that the contribution atw −2 in Eq. (13) ex-hibits very weak residual dependence on the coupling λ.

V. SUMMARY & OUTLOOK
We have studied the transition to hydrodynamics of Yang-Mills kinetic theory undergoing Bjorken expansion. For initial conditions inspired by the CGC effective theory, we analyzed the subsequent evolution of a natural set of four observables given by Eq. (15) using PCA. Such an analysis may be sensitive to the precise initial conditions used, which is why we studied realistic initial conditions from the CGC effective theory.
Our studies showed that at late times there is a single dominant principal component that describes the variation of the overall energy scaleP T . Furthermore, the late time evolution ofP L ,m 2 D andT * is highly correlated and can be captured by a second subleading principal component. While at this point it remains a logical possibility that this correlation is an artifact of the considered initial conditions, this could be further explored within a higher dimensional parameter space.
The second principal component represents a transient contribution, which we studied in more detail by consid-eringP L using other initial conditions initialized at different times. A clear exponential decay can be seen at late times by comparing different profiles ofP L , signaling the presence of an exponential approach to the hydrodynamic attractor at late times.
Quite surprisingly, for the whole range of the couplings considered, i.e. for λ between 0.1 and 20, the decay rate is close to the value predicted by the conformal RTA kinetic theory. What this means in practical terms is that the effective decay rate is simply expressible in terms of the shear viscosity-to-entropy density ratio, see Eq. (21).
It is an interesting question to ask if the similarities between the EKT and conformal RTA go beyond this crude characteristic. To this end, an in-depth analysis of the RTA kinetic theory in [34] reveals that the exponentially decaying contribution there is not a single excitation, but actually a sum of infinitely many contributions whose occupation numbers are initial condition dependent. While they all are characterized by the same exponential decay rate, their subleading behaviour is distinct yet impossible to disentangle over the short time scales probed in the present project. This indicates that in reality the exponential decay we are reporting here is likely to be understood as an effective description (appropriate for the conformal Bjorken flow with no transverse dynamics) of an underlying more complicated structure.
Beyond further explorations of the surprising similarity in the highly symmetric Bjorken flow, it would also be interesting to compare the evolution of the energymomentum tensor in a variety of other settings. While first steps in this direction have been reported in [42,71] for linearized perturbations around Bjorken flow and in [29,72,73] for systems undergoing both longitudinal and transverse expansion, a particularly clean probe concerns  Table I for numerical values.
the evolution of energy-momentum (or metric) perturbations in equilibrium. While for conformal RTA kinetic theory the corresponding retarded correlators of the energy-momentum tensor were calculated analytically in [74], a dedicated study in QCD kinetic theory would certainly shed further light on the structure and importance of non-hydrodynamic excitations. One can view the results of [68] obtained in RTA kinetic theory with momentum-dependent relaxation time (see, however, [75] for a subtlety in this model that emerged after [68] was completed) as a first step in this direction. By inverting Eq. (A.1) for η/s, we can use the ratio 9/16(1/3−P L / )T τ to extract the ratio of shear viscosity to entropy density η/s from the asymptotic behavior of each numerical solution, and in subsequent analysis we use the mean value for each λ. We illustrate this extraction in the left panel of Fig. 7, while the right panel of Fig. 7 shows the extracted values in comparison to NLL parametrization of η/s in [76]. We note that for small couplings λ, the extracted values of η/s are in excellent agreement with the NLL parametrization, while for larger values of λ the NLL parametrization breaks down, and as discussed in [44] the transport coefficients also become more sensitive to the precise implementation of the screening of the elastic matrix elements. We also provide the extracted values in Table I.