# Designing Pu Intermetallics with First Principle Calculations

Matthew S. Cook† David C. Arellano† Derek V. Prada† Sven P. Rudin†† Eric D.
Bauer††† W. Adam Phelan† †MST-16, Los Alamos National Laboratory, Los Alamos,
New Mexico 87545, USA  ††T-1, Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, USA  †††MPA-Q, Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, USA  *corresponding author

###### Abstract

We present the ab initio supported discovery of two new Pu based intermetallic
compounds, PuNiSn and PuPtSn. Using density functional theory, the formation
energies within the relevant ternary phase diagrams were calculated to predict
the stability of both compounds. Simultaneously, Pu-Ni-Sn and Pu-Pt-Sn
materials were arc-melted and subsequently characterized with magnetization,
specific heat, and resistivity measurements from 2-300 K. Magnetization
measurements show that PuNiSn and PuPtSn order antiferromagnetically at
$T_{N}=11\ K$ and $T_{N}=15\ K$, respectively. Specific heat measurements show
an enhanced residual electronic specific heat that is indicative of strong
electron correlations. Resistivity measurements are indicative of Kondo
behavior for PuNiSn while crystal field effects may play a role in the
observed temperature dependence for PuPtSn.

††preprint: APS/123-QED

##  I Introduction

Plutonium (Pu) is a unique element with a complex phase diagram that hosts
seven allotropic phases [1, 2], an electronically driven 25% volume
contraction from the $\delta$ phase to the low temperature $\alpha$-phase [3],
and an anomalously large residual electronic heat capacity has been observed
for its meta-stable alloy phases (e.g., Ga-stabilized $\delta$-Pu) [4]. Many
of theses extraordinary properties of Pu are generally associated with the
unstable nature of its 5f electron configuration, which resides at the border
between itinerant and localized behavior in the periodic table.

This complexity is also clearly expressed in Pu compounds, where the 5f-states
can vary from 5f2 to 5f6 [5], or non-integer filling [6, 7, 8]. In this
regard, Pu-based intermetallic compounds have attracted interest for their
complex physical behavior which results from strong hybridization of 5f-states
with the conduction electron states of the surrounding ligands. The strong 5f
correlations in Pu intermetallics often result in compounds with an enhanced
electronic specific heat coefficient, which corresponds to a large effective
mass of charge carriers [9]. Various magnetic ground states have also been
observed in Pu compounds, such as ferromagnetism [10], antiferromagnetism
[11], and paramagnetism [12]. The four known examples of Pu-based
superconductors, PuCoGa5, PuRhGa5, PuCoIn5, and PuRhIn5 are characterized by
enhanced electronic specific heat coefficients and relatively large
superconducting transition temperatures for this superconducting class [13,
14, 15, 11]. For heavy fermion superconductors, PuCoGa5 was found to have a
record high transition temperature of 18.5 K, which is several times larger
than isostructural Ce based counterparts, CeCoIn5 and CeIrIn5 [16, 17].

From the preceding statements, the discovery of new Pu-containing materials is
a worthy pursuit. Yet, the experimental search for new Pu-based intermetallics
is heavily constrained by the acquisition of starting material, its associated
safety hazards concerning radioactivity, and the limited laboratory spaces
equipped to handled Pu. On the other hand, the theoretical search for such
intermetallics is limited by the challenges that the strong 5f correlations
pose on methods. Thus, it is vital that experimentalists and theorists work
together to make informed predictions on the stability of potential new
phases. In this work, by using density functional theory (DFT) to calculate
thermodynamics parameters of various ternary and binary Pu intermetallic
compounds in the Pu:T:Sn (T = transition metals) phase space, we attempt to
marry both experiment and theory in the hopes of laying ground for future
theory driven Pu-materials discovery.

In addition to the aforementioned DFT component, we report the synthesis and
physical properties characterization of two new Pu-based intermetallics,
namely, PuNiSn and PuPtSn. These materials expand the phase space of Pu 1:1:1
intermetallic compounds, as similar PuTX (T = transition metals and X = Al,
Ga) compounds were first synthesized nearly 5 decades ago [18]. More recently,
the physical characterization of PuPdSn was reported by Gofryk et al. [19].
Notably, PuNiSn crystallizes in an orthorhombic TiNiSi-type structure (space
group $Pnma$) while PuPtSn shares the same structure as the Pd analog.
Further, using magnetization we find that PuNiSn orders antiferromagnetically
at $T_{N}=11\ $K while PuPtSn also shows an antiferromagnetic transition at
$T_{N}=15\ $K with a potential secondary transition at 3 K. Curie-Weiss fits
to the high temperature magnetic susceptibility show that the effective moment
for PuNiSn ($\mu_{eff}=0.93\ \mu_{B}$) is close to that expected for an
intermediate coupling scheme, while PuPtSn ($\mu_{eff}=0.77\ \mu_{B}$) is
slightly lower than that expected in a Russell-Saunders coupling scheme. Low
temperature fits to $C/T$ vs $T^{2}$ taken above the magnetic transitions find
a Sommerfeld coefficient $\gamma$ of $243(3)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$ for PuNiSn and $216(3)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$ for PuPtSn. Resistivity measurements show a weak temperature
dependence in both materials with clear Kondo behavior present in PuNiSn.
Finally, using DFT, we calculate the phonon dispersion and corresponding
partial density of states for both compounds and compare to the lattice modes
determined from the specific heat. We find that two Debye modes can
effectively capture the vibrational frequencies coming from the consituent
atoms due to their large difference in mass.

##  II Methods

DFT calculations were performed to determine total energies and phonon
dispersions/density of states for a variety of materials within the Pu-Ni-Sn
and Pu-Pt-Sn phase diagrams. The DFT calculations allowed (collinear) spin
polarization (orbital polarization is not included). The computational cells
contained multiple Pu atoms, and so all possible spin structures were tested.
Here we allow spin polarization because the ground state of Pu compounds often
have a net magnetic moment. The inclusion of spin polarization has also been
shown to result in a more accurate description of elemental Pu in
approximation of the strong correlations present in the 5f electronic
structure [20]. For the two ternary compounds reported here, the spins on all
Pu atoms aligned in the same direction (ferromagnet). Spin-orbit coupling is
neglected, as its effects on the results for representative structures show
negligible changes. The DFT calculations were performed using the VASP package
[21, 22]. The electronic states were treated with the projector augmented wave
method [23], in the generalized gradient approximation of Perdew, Burke, and
Ernzerhof [24], with first-order Methfessel-Paxton smearing (width 27 meV)
[25], and a cutoff energy for the plane wave basis set of 500 eV. The k-point
meshes were chosen for each simulation cell such that their density is at
least 60 per Å-1. The convergence criteria are 10-7 eV for the electronic
self-consistency loop and 10-6 eV for the ionic relaxation. The phonon
calculations employ 2x2x2 supercells. Ternary phase diagrams were explored by
calculating formation energies.[26] Since the focus is on the central part of
the phase diagrams, the complete convex hull is not calculated. Instead,
formation energies are calculated relative to known, nearby compounds.

Figure 1:  (a) Powder x-ray diffraction data (red circles) for PuNiSn with
LeBail fit overlain in black. (b) Powder x-ray diffraction data (red circles)
for PuPtSn with LeBail fit overlain in black. The difference between the
calculated and experimental data are shown in blue in both panels.

Polycrystalline samples of PuNiSn and PuPtSn were arc-melted on stoichiometry
in an Zr-gettered argon atmosphere. Each button was flipped and remelted twice
to help ensure homogenization. Following arc-melting, the buttons were wrapped
in Ta-foil and sealed into separate evacuated silica tubes. The tubes were
placed into a muffle furnace at 773 K for two weeks.

Powder X-ray diffraction data were collected at room temperature using a
Bruker D8 Discover with DAVINCI design powder diffractometer equipped with a
LynxEye detector and a Cu K$\alpha$ radiation source. Subsequent LeBail fits
were generated using Bruker’s TOPAS v. 6 sofware. PuNiSn crystallizes in an
orthorhombic TiNiSi-type structure (space group $Pnma$), as shown in Fig.
1(a). PuPtSn crystallizes in a hexagonal ZrNiAl-type structure (space group
$P\overline{6}2m$), as shown in Fig. 1(b).

Specific heat measurements were collected using a commercial physical
properties measurement system (PPMS) cryostat via a semi-adiabatic thermal
relaxation technique. As Pu is a radioactive element, decay in this material
causes self-heating. This limits the temperature of which we were able to cool
the sample, such that larger samples were more difficult to cool. For PuNiSn,
the mass of the measured sample is 9 mg, thus we were only able to measure to
approximately $T=5.5\ \mathrm{K}$. For a sample mass of only 1.7 mg for
PuPtSn, this allowed us to reach $T=2.5\ \mathrm{K}$. Magnetization
measurements were also performed in the same commercial PPMS cryostat using a
vibrating sample magnetometer. In order to safely measure Pu samples and to
prevent possible contamination hazards, the magnetization samples were mounted
using layers of encapsulation. This includes small amounts of plastic wrap and
GE varnish glued inside a plastic straw. The plastic straw was then sealed
with custom caps that contain micron sized platinum frits which were glued to
the straw with stycast epoxy. Resistivity measurements were done with a
standard 4-wire configuration in the PPMS Dynacool cryostat using native
electronics.

##  III Results and Discussion

###  III.1 Convex Hull Calculations

To identify potential candidates for Pu-materials discovery, we determine the
theoretical thermodynamic stability of various compounds in the Pu-Ni-Sn and
Pu-Pt-Sn phase spaces using formation energy calculations and a convex hull
construction [26]. We calculate the formation energy of a given compound of
composition PulTmSnn relative to its constituent elements as:

$\displaystyle\Delta{E_{f}}={E_{\mathrm{Pu}_{l}\mathrm{T}_{m}\mathrm{Sn}_{n}}-lE_{\mathrm{Pu}}-mE_{\mathrm{T}}-nE_{\mathrm{Sn}}},$
(1)

where $E_{\mathrm{Pu}_{l}\mathrm{T}_{m}\mathrm{Sn}_{n}}$ is the total energy
per formula unit of $\mathrm{Pu}_{l}\mathrm{T}_{m}\mathrm{Sn}_{n}$, while
$E_{\mathrm{Pu}}$, $E_{\mathrm{T}}$, and $E_{\mathrm{Sn}}$ are the total
energies per atom of the constituent elements.

Figure 2:  DFT predicted phase diagrams for the (a) Pu-Ni-Sn and (b) Pu-Pt-Sn
systems. The phases which fall on a vertex of the facets are expected to be
thermodynamically stable compounds.

In Fig. 2, we plot each phase for the (a) Pu-Ni-Sn and (b) Pu-Pt-Sn series at
its particular composition in a two-dimensional diagram. We then form both
convex hulls using the calculated formation energies of each composition
relative to their constituent elements. This allows for identification of
competing phases in a relative compositional space. We also calculate the
formation energy ($E_{comp}$) with respect to these competing phases to
determine if decomposition is energetically favorable and ultimately determine
a compound’s theoretical thermodynamic stability at zero temperature. The
predicted stable compounds from this process lie on the vertex of each facet
in Fig. 2. Our calculations show that both PuNiSn and PuPtSn are predicted to
be thermodynamically stable relative to nearby competing phases. Table 1
presents the complete results of this analysis for both compounds.

| $E_{int}$ | $\Delta{E_{f}}$ | $\Delta{E_{comp}}$ |  |
---|---|---|---|---|---
Compound | [eV/f.u.] | [meV/atom] | [meV/atom] | Competing Phases | Space Group
PuNiSn | -24.377 | -343 | -67 | Pu2Ni2Sn, Pu3Ni3Sn4 |  $Pnma$ (62)
Pu2Ni2Sn | -44.339 | 296 | -4 | PuNiSn, PuNi2, Pu |  $P4/mbm$ (127)
Pu3Ni3Sn4 | -76.542 | -266 | 40 | PuNiSn, Sn |  $I\bar{4}3D$ (220)
Pu3Sn | -46.708 | -185 | -49 | Pu, PuSn2 |  $Pm\bar{3}m$ (221)
PuSn2 | -22.801 | -363 | -185 | Pu3Sn, PuSn3 |  $I4_{1}/amd$ (141)
PuSn3 | -26.946 | -350 | -77 | PuSn2, Sn |  $Pm\bar{3}m$ (221)
PuNi2 | -25.837 | -286 | -217 | Ni, PuNi |  $Fd\bar{3}m$ (227)
PuNi | -19.719 | -103 | 111 | Pu, PuNi2 |  $Cmcm$ (63)
Ni3Sn | -21.004 | -192 | -22 | Ni3Sn2, Ni |  $Pm\bar{3}m$ (221)
Ni2Sn | -15.269 | -167 | 69 | Ni3Sn, Ni3Sn2 |
Ni3Sn2 | -25.427 | -271 | -38 | Ni3Sn, NiSn |  $Pnma$ (62)
NiSn | -9.822 | -260 | -11 | Ni3Sn2, Ni3Sn4 |  $Pbam$ (55)
Ni3Sn4 | -33.374 | -234 | -11 | NiSn, Sn |  $C2/m$ (12)
|  |  |  |  |
PuPtSn | -26.187 | -738 | -169 | Pu2Pt2Sn, PtSn, PuSn2, PuPt |  $P\bar{6}2m$ (189)
Pu2Pt2Sn | -47.245 | -627 | 126 | PuPtSn, PuPt |  $P4/mbm$ (127)
Pu2Pt | -35.629 | -480 | 34 | Pu, PuPt |  $Pnma$ (62)
Pu5Pt3 | -93.018 | -562 | 16.2 | Pu,PuPt |  $P63/mcm$ (193)
PuPt | -21.686 | -771 | -217 | Pu,PuPt2 |  $Cmcm$ (63)
PuPt2 | -28.459 | -738 | -26 | PuPt, PuPt3 |  $Fd\bar{3}m$ (227)
PuPt3 | -35.074 | -683 | -129 | PuPt2, Pt |  $Pm\bar{3}m$ (221)
Pu3Sn | -46.708 | -185 | -49 | Pu, PuSn2 |  $Pm\bar{3}m$ (221)
PuSn2 | -22.801 | -363 | -185 | Pu3Sn,PuSn3,Sn |  $I4_{1}/amd$ (141)
PuSn3 | -26.946 | -350 | -77 | PuSn2, Sn |  $Pm\bar{3}m$ (221)
Pt3Sn | -23.662 | -383 | -90 | Pt, PtSn |  $Pm\bar{3}m$ (221)
PtSn | -11.102 | -585 | -129 | Pt3Sn, Pt2Sn3, PtSn2 |  $P6_{3}/mmc$ (194)
Pt2Sn3 | -26.268 | -514 | -6 | PtSn, PtSn2 |  $P6_{3}/mmc$ (194)
PtSn2 | -15.133 | -456 | -18 | Pt2Sn3, PtSn4 |  $Fm\bar{3}m$ (221)
Table 1: Relevant phases with associated space groups and calculated energies
for Pu-Ni-Sn and Pu-Pt-Sn phase diagrams. Determined competing phases for
which the formation energy $\Delta{E_{comp}}$ is calculated are also listed.

###  III.2 Magnetization

The magnetic susceptibility $\chi=M/H$ and its inverse $\chi^{-1}$ are shown
in Figs 3(a) and 3(d) for PuNiSn and PuPtSn, respectively. Both compounds were
fit with the Curie-Weiss law given by:

$\displaystyle\chi(T)=\chi_{0}+\frac{C}{T-\theta_{C}},$ (2)

where $\chi_{0}$ is the temperature independent contribution that results from
a combination of Pauli paramagnetism, Van Vleck paramagnetism, and the sample
holder background, $\theta_{C}$ is the Curie-Weiss temperature, and $C$ is
related to the effective moment $\mu_{eff}=\sqrt{8C}\ \mu_{B}$. The results of
the fit give $\theta_{C}=-37\ K$, $\mu_{eff}=0.93\ \mu_{B}$, and
$\chi_{0}=161\times 10^{-6}$ emu mol-1 for PuNiSn, and $\theta_{C}=-26\ K$,
$\mu_{eff}=0.77\ \mu_{B}$, and $\chi_{0}=-1.80\times 10^{-5}$ emu mol-1 for
PuPtSn.

The effective moment for PuNiSn is close to the value expected for
intermediate coupling (1.1 $\mu_{B}$) and similar to that found for Pu2Ni3Si5
[27]. An intermediate coupling scheme has been previously shown to be
appropriate in the description of 5f orbitals for f5 Pu compounds [28]. The
negative $\theta_{C}=-37\ K$ is indicative of antiferromagnetic correlations
that lead to an antiferromagnetic transition that can been seen in magnetic
susceptibility at $T_{N}=11\ K$ in Fig. 3(b). Here, $T_{N}$ is defined as the
maximum in the derivative $d\chi/dT$ (not shown). A magnetic field of 14 T has
a very small effect on the transition temperature.

For PuPtSn, the effective moment $\mu_{eff}=0.77\ \mu_{B}$ is slightly below
the value expected in an LS coupling scheme for a $5f^{5}$ manifold
($0.84\mu_{B}$). Although the intermediate coupling scheme is generally
accepted for Pu3+ compounds due to strong spin orbit coupling, lower values
for the effective moment have also been observed in other Pu intermetallics
[29]. The magnetic transition at $T_{N}=15\ K$ is in agreement with the
negative $\theta_{C}$ and also has a negligible dependence on magnetic field,
as shown in Fig. 3(e). The AFM transition in PuPtSn is extremely broad
compared to PuNiSn and this character is also observed in the specific heat.
There is a secondary transition in $\chi$ that is seen as an upturn in $\chi$
vs $T$ at approximately 4 K. This transition is suppressed by the application
of a field of 9 T.

The magnetism and strong electron correlations of related uranium ternary
compounds with the orthorhombic TiNiSi-type and hexagonal ZrNiAl-type
structures have been previously discussed in the context of electron
hybridization that is a direct consequence of the large spacial extent of the
5f wavefunction [18]. The hybridization of 5f states with valence states of
the surrounding atoms in these materials typically results in the
delocalization of 5f electrons and the formation of a narrow 5f band that
intersects the Fermi energy ($E_{f}$). This leads to magnetic moments that are
much smaller than expected for a free U3+ ion as well as a large density of
states at $E_{f}$ that is reflected in an enhanced electronic specific heat
coefficient $\gamma$. Here, the effective moment estimated in both compounds
is not greatly reduced from the expected moment for a free Pu3+ ion, thus we
can deduce that the 5f electrons in these materials have a localized
character.

Figure 3:  The magnetic susceptibility $\chi$ (left axis) and corresponding
inverse magnetic susceptiblity $\chi^{-1}$ (right axis) for (a) PuNiSn and (d)
PuPtSn. The solid red lines correspond to Curie-Weiss fits to the data.
Antiferromagnetic transitions in the low temperature magnetic susceptibility
$\chi$ at various magnetic fields for (b) PuNiSn and (e) PuPtSn. There is a
secondary transition near 4 K for PuPtSn in a field of 3 T. The field
dependence of the magnetization $M$ at various temperatures for (c) PuNiSn and
(f) PuPtSn.

###  III.3 Specific Heat

Figure 4:  Specific heat of (a) PuNiSn and (d) PuPtSn with corresponding fits
to equation (3) in solid red. Insets: Peak in $C/T$ vs $T^{2}$ from AFM order
with linear fit in solid red. Magnetic specific heat divided by temperature is
plotted on the left axis of (b) for PuNiSn and (e) PuPtSn, while the magnetic
entropies are shown in red on the right axes. The magnetic field dependence of
the specific heat divided by temperature for (c) PuNiSn and (f) PuPtSn.

The zero-field specific heat of PuNiSn and PuPtSn for $5\ \mathrm{K}\leq 300\
\mathrm{K}$ is shown in Fig. 4(a) and Fig. 4(d), respectively. For PuNiSn,
there is an anomaly in the specific heat at $T_{N}=11\ \mathrm{K}$ that
corresponds to the onset of antiferromagentic order. For PuPtSn, the ordering
temperature occurs slightly higher, at $T_{N}=15\ \mathrm{K}$. This is in
contrast to the isostructural PuPdSn which shows two magnetic transitions in
specific heat and magnetization measurements at $T_{1}=21\ $K and $T_{2}=9.6\
$K [19]. To get an estimate of the Sommerfeld coefficient at low temperatures,
the specific heat of both compounds is plotted as $C/T$ vs $T^{2}$, as shown
in the inset of Figs 4(a) and (d). We fit the linear region above $T_{n}$ to
$C/T=\gamma+{\beta}T^{2}$, where $\gamma$ is the Sommerfeld coefficient and
$\beta=1944/\Theta_{D}^{3}$, where $\Theta_{D}$ represents the Debye
temperature. The fit for PuNiSn yields $\gamma=243(3)\ \mathrm{mJ}\
\mathrm{mol}^{-1}\ \mathrm{K}^{-2}$ and $\beta=0.653(8)\ \mathrm{mJ}\
\mathrm{mol}^{-1}\ \mathrm{K}^{-4}$. Similarly, the fit for PuPtSn yields
$\gamma=216(2)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\ \mathrm{K}^{-2}$ and
$\beta=0.540(5)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\ \mathrm{K}^{-4}$. The
estimated $\gamma$ for both compounds is fairly enhanced when compared with
PuPdSn, which is estimated as $8\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$ [19]. However, the value found for PuPdSn must be met with
caution as it was determined by fitting the specific heat below the magnetic
transition with an electronic and debye term. To make a more direct comparison
to PuPdSn, we digitize the specific heat data in Ref. [19] and determine
$\gamma$ via extrapolation to 0 K in the paramagnetic state (not shown). Using
this method, we find a much larger $\gamma=290$ mJ mol-1 K-2, which is closer
to the values found for PuNiSn and PuPtSn. The values estimated for $\gamma$
are also similar to other magnetically ordered plutonium based heavy-fermion
systems, such as PuIn3 ($307\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$) [30], PuAl2 ($260\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$) [31], and Pu2Co3Si5 ($95\ \mathrm{mJ}\ \mathrm{mol}^{-1}\
\mathrm{K}^{-2}$), indicative of strong electron correlations.

In order to determine an accurate contribution to the specific heat from the
lattice, we employed a fit from $25\ \mathrm{K}\leq 300\ \mathrm{K}$. The
magnetic transitions at low temperatures causes difficulty in getting a
reliable fit over the entire temperature range, thus the chosen temperatures
are well above the transition. We assume that the form of the specific heat is
given by $C=C_{e}+(n-r)C_{D1}+(r)C_{D2}$, where $n=3$ is the number of atoms
per formula unit, r is a fitting parameter used to keep the amount of
oscillators consistent, $C_{e}=\gamma{T}$ where $\gamma$ is the Sommerfeld
coefficient, and $C_{D}$ is a Debye term given by

$C_{D}=9Nk_{B}x^{-3}\int_{0}^{x}\frac{x^{4}e^{x}}{(e^{x}-1)^{2}}\,dx$ (3)

where $T_{D}$ is the Debye temperature and $x=T_{D}/T$. Attempts to fit the
specific heat using a single Debye mode were attempted, yet were unsuccessful.
The resulting fits that are plotted as the solid red lines in Fig. 4(a) and
(d) yield $T_{D1}=259(10)\ \mathrm{K}$, $T_{D2}=113(12)\ \mathrm{K}$,
$\gamma=16.2(5)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\ \mathrm{K}^{-2}$, $r=1.0(2)$
for PuNiSn and $T_{D1}=266(11)\ \mathrm{K}$, $T_{D2}=141(8)\ \mathrm{K}$,
$\gamma=5.2(5)\ \mathrm{mJ}\ \mathrm{mol}^{-1}\ \mathrm{K}^{-2}$, $r=1.3(2)$
for PuPtSn. We will further comment on the lattice specific heat contributions
in the next section in relation to the phonon density of states.

To isolate the magnetic contribution to the specific heat, we subtract the fit
to the high temperature specific heat such that
$C_{mag}=C-(C_{e}+(n-r)C_{D1}+(r)C_{D2})$. To get an estimate of the magnetic
entropy associated with the AFM order in this system, we tentatively
extrapolate $C_{mag}/T$ vs $T$ linearly to $T=0\ \mathrm{K}$, as shown in blue
on the left axes in Figs 4(b) and (e). $C_{mag}(T)/T$ vs $T$ is then
numerically integrated to calculate the magnetic entropy $S_{mag}(T)$, which
we plot on the right axes of Fig. 4(b) and (e).

For both compounds, it is possible that crystalline electric field (CEF)
effects give rise to a doublet ground state; the J = 5/2 multiplet of a Pu3+
ion in an orthorhombic or hexagonal site symmetry split into 3 Kramers
doublets due to crystal field effects. Here, the estimated magnetic entropy
for both PuNiSn and PuPtSn is reduced from the Rln(2) expected for a ground
state doublet in the CEF. This may be due to similar energy scales between the
Kondo effect and RKKY interation that both depend on the coupling parameter
$J$, which describes the hybridization between conduction electrons and local
moments. Magnetic order was observed in Pu2Ni3Si5 and Pu2Co3Si5 with a reduced
specific heat jump expected for a spin 1/2 system, reduced entropy expected
for a ground state doublet, and a large linear specific heat coefficient both
above and in the magnetically ordered state [27]. This is thought to be due to
Kondo and RKKY interactions of similar magnitude and is based on theoretical
work which describes magnetic Kondo compounds with a doublet crystal field
ground state using an $S=1/2$ resonant level model [32]. This work predicts a
reduced specific heat jump at the antiferromagnetic ordering temperature as a
function of $T_{K}$/$T_{N}$ and is in agreement with experimental evidence
from various Ce and Yb compounds. In estimating the jumps in specific heat in
Fig. 2, we find they are reduced from the expected 12.48 J mol-1 K-1 for both
PuNiSn ($\Delta{C(11\ \mathrm{K}})\approx{3}$ J mol-1 Pu K) and PuPtSn
($\Delta{C(15\ \mathrm{K}})\approx{3.4}$ J mol-1 Pu K). This fact along with
reduced magnetic entropy may indicate similar energy scales for Kondo and RKKY
interactions in these materials.

The application of magnetic field has an extremely small effect on the peak in
the specific heat of both compounds, as shown in Figs 4(c) and (f). For
PuNiSn, a field of $\mu_{0}H=14\ T$ pushes the transition to lower
temperatures by only 0.25 K and slightly in magnitude. The peak for PuPtSn
behaves very similarly in field and additionally has some field dependence at
lower temperatures. The field dependence at low temperatures appears to be
from a broad peak centered at approximately 3 K that is broadened with field
as well as enhanced in temperature and magnitude. This is likely related to
the small upturn seen in magnetic susceptibility at the same temperature.

###  III.4 Resistivity

Figure 5(a) shows the temperature dependence of the resistivity for PuNiSn.
PuNiSn shows typical Kondo behavior of a minimum in the resistivity followed
by an increase as the temperature is lowered. There is a sharp decrease in the
resistivity from a reduction in magnetic scattering as the moments order
antiferromagnetically. As shown in the inset of Fig. 5(a), the peak in
resistivity coincides with the peak in magnetic susceptibility at $T_{N}=10$
K. NpNiSn, which also orders antiferromagnetically at 37 K, shows a similar
overall temperature dependence [33]. However, the resistivity in the Np
analogue has a negative slope from high temperatures down to the magnetic
ordering temperature, which is thought to be from dynamic spin correlations
rather than Kondo scattering.

The resistivity for PuPtSn is plotted in Fig. 5(b). It shows a very different
temperature dependence than PuNiSn, where we find a monotonic decrease in
resistivity down to low temperatures and a shallow minimum around 20 K, likely
related to Kondo scattering or magnetic fluctuations in proximity to $T_{N}$.
The maximum in resistivity at 12 K coincides well with the maximum in magnetic
susceptibility. The temperature dependence of the resistivity is comparable to
CePtSn, a Ce3+ ($J=5/2$) analogue that crystallizes in a TiNiSi-type structure
[34]. CePtSn is a Kondo compound with $T_{K}/T_{N}\approx{1}$, where crystal
field splitting has an effect on the temperature dependence of the
resistivity, as it was determined from inelastic neutron scattering that the
first two excited states are $\Delta_{1}=273\ K$ and $\Delta_{2}=425\ K$ [35].
It is possible here that the thermal population of crystal field levels in
PuPtSn also plays a role in the observed temperature dependence.

Figure 5:  Resistivity versus temperature for (a) PuNiSn and (b) PuPtSn.
Insets: Manifestations of the AFM transitions show a sudden drop in
resistivity in both compounds.

###  III.5 Phonon Calculations

Figure 6:  Phonon dispersion (left) and corresponding density of states
(right) for (a) PuNiSn and (b) PuPtSn. The dotted lines represent the
experimental Debye temperatures from fits to the specific heat using equation
(3).

The calculated phonon dispersion and partial phonon density of states (PPDOS)
for PuNiSn and PuPtSn are presented in Fig. 6(a) and Fig. 6(b), respectively.
For PuNiSn, the high frequency region is characterized by a gap, likely due to
the large difference in mass of Pu compared to Ni and Sn, as relatively no
phonon density comes from Pu above 5 THz. It has been previously shown that
the phonon band gaps present in nitrides within the wurtzite structure can be
increased by changing the cation-anion mass mismatch [36]. Upon inspection of
the PPDOS, the different vibrational properties of the constituent elements
become apparent. The low frequency spectrum, especially below 2.5 THz, is
mainly due to the vibration of Pu atoms, while the high frequency spectrum
above 3.5 THz mainly comes from the lighter elements, Ni and Sn. There is
however some Pu phonon density that coincides with the lighter elements in
between these two regimes. This suggests the existence of phonon modes where
all three elements vibrate together. In comparing this to fits to the lattice
contribution of the specific heat we find that the two reported Debye
temperatures which are show as dotted lines in Fig. 6 correspond fairly well
to the separate densities of the constituent elements. Our specific heat fit
produces a Debye temperature of 113 K with $r=1$ which is consistent with Pu
being responsible for most of the low frequency phonons.

Similar arguments can be made for PuPtSn with some clear differences. We still
find a gap in the phonon spectrum which clearly separates the Sn contributions
in the higher frequencies from the Pu and Pt contributions in the lower
frequencies. In this material, the transition element Pt (195.08 u) is much
closer in mass to Pu (244.06 u) than Ni (58.69 u). Therefore we see in the
phonon PDOS that most of the density of the lower frequencies are from Pu and
Pt. This difference is also represented in the fit to the specific heat, where
now $r=1.3$ and the smaller Debye temperature $T_{D2}$ has changed to 141 K as
the Pu and Pt phonons occupy higher frequencies. Here, two Debye temperatures
seem to capture the expected phonon density of states more accurately than for
PuNiSn, as the dashed lines clearly show.

##  IV Conclusion

In conclusion, we have successfully used DFT calculations to predict the
thermodynamic stability of two new Pu intermetallics, PuNiSn and PuPtSn. We
have shown that exploratory synthesis of Pu compounds can be well informed
through theoretical calculations, which offers a path forward when considering
the difficulties associated with elemental Pu. We also show that our
calculated phonon PDOS agree fairly well with modeling of the lattice specific
heat. Using arc melting, we have grown these two compounds, confirmed the
structures using PXRD, and subsequently characterized both materials using
specific heat, magnetization, and resistivity measurements. Magnetization
measurements show Curie-Weiss behavior at high temperatures followed by
antiferromagnetic ordering of Pu3+ moments in both compounds at low
temperatures. The low temperature anomalies in specific heat also support
magnetic order, while a reduced magnetic entropy and jumps in the specific
heat suggest a competition between the on-site Kondo and intersite RKKY
interactions. Fits to the specific heat above the ordering temperature give an
enhanced linear specific heat coefficient which is indicative of strong
electron correlations in these materials. Resistivity measurements corroborate
magnetic ordering in both PuNiSn and PuPtSn, as there is a prototypical drop
in resistivity upon entering the ordered state from a reduction in magnetic
scattering.

###### Acknowledgements.

All authors are grateful for support from LDRD’s early career research
catorgory under proposal number 20220538ECR. Further all authors are grateful
to P.H. Tobash and L.M. Casados for useful discussions. E.D.B. was supported
by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of
Materials Science and Engineering, under the project ”Quantum Fluctuations in
Narrow Band Systems.”

## References

  * Tróc and Fournier [1985] R. Tróc and J. M. Fournier, Bulk Properties of the Actinides, in  _Handbook of the Physics and Chemistry of the Actinides_ , Vol. 2, edited by A. J. Freeman and G. H. Lander (North-Holland, Amsterdam, 1985) pp. 29–174.
  * Hecker [2004] S. S. Hecker, The magic of plutonium: 5f electrons and phase instability, Metallurgical and Materials Transactions A 35, 2207 (2004).
  * Harrison  _et al._ [2019] N. Harrison, J. B. Betts, M. R. Wartenbe, F. F. Balakirev, S. Richmond, M. Jaime, and P. H. Tobash, Phase stabilization by electronic entropy in plutonium, Nature Communications 10, 3159 (2019).
  * Lashley  _et al._ [2003] J. C. Lashley, J. Singleton, A. Migliori, J. B. Betts, R. A. Fisher, J. L. Smith, and R. J. McQueeney, Experimental Electronic Heat Capacities of $\alpha$ \- and $\delta$ -Plutonium: Heavy-Fermion Physics in an Element, Physical Review Letters 91, 205901 (2003).
  * N/A [2006] N/A,  _The chemistry of the actinide and transactinide elements_ , 3rd ed., edited by N/A (Springer, Dordrecht, 2006).
  * Joyce  _et al._ [2006] J. J. Joyce, J. M. Wills, T. Durakiewicz, M. T. Butterfield, E. Guziewicz, D. P. Moore, J. L. Sarrao, L. A. Morales, A. J. Arko, O. Eriksson, A. Delin, and K. S. Graham, Dual nature of the 5f electrons in plutonium materials, Physica B: Condensed Matter 378-380, 920 (2006).
  * Shick  _et al._ [2009] A. Shick, L. Havela, T. Gouder, and J. Rebizant, Electronic-structure theory of plutonium chalcogenides, Journal of Nuclear Materials 385, 21 (2009).
  * Ramshaw  _et al._ [2015] B. J. Ramshaw, A. Shekhter, R. D. McDonald, J. B. Betts, J. N. Mitchell, P. H. Tobash, C. H. Mielke, E. D. Bauer, and A. Migliori, Avoided valence transition in a plutonium superconductor, Proceedings of the National Academy of Sciences 112, 3285 (2015).
  * Bauer and Thompson [2015] E. Bauer and J. Thompson, Plutonium-Based Heavy-Fermion Systems, Annual Review of Condensed Matter Physics 6, 137 (2015).
  * Griveau  _et al._ [2013] J.-C. Griveau, E. Colineau, D. Bouëxiere, K. Gofryk, T. Klimczuk, and J. Rebizant, Magnetic properties of ferromagnetic Pu2Pt3Si5, Journal of Alloys and Compounds 576, 409 (2013).
  * Bauer  _et al._ [2012a] E. D. Bauer, P. H. Tobash, J. N. Mitchell, and J. L. Sarrao, Single crystal growth of plutonium compounds from molten metal fluxes, Philosophical Magazine 92, 2466 (2012a).
  * Wastin  _et al._ [1993] F. Wastin, J. Rebizant, J. Spirlet, J. Fuger, B. Kanellakopulos, and V. Sechovsky, Magnetic study of new Pu(Mt: 3d, 4d or 5d transition metal)2Si2 intermetallic compounds, Journal of Alloys and Compounds 193, 119 (1993).
  * Sarrao  _et al._ [2002] J. L. Sarrao, L. A. Morales, J. D. Thompson, B. L. Scott, G. R. Stewart, F. Wastin, J. Rebizant, P. Boulet, E. Colineau, and G. H. Lander, Plutonium-based superconductivity with a transition temperature above 18 K, Nature 420, 297 (2002).
  * Wastin  _et al._ [2003] F. Wastin, P. Boulet, J. Rebizant, E. Colineau, and G. H. Lander, Advances in the preparation and characterization of transuranium systems, Journal of Physics: Condensed Matter 15, S2279 (2003).
  * Bauer  _et al._ [2012b] E. D. Bauer, M. M. Altarawneh, P. H. Tobash, K. Gofryk, O. E. Ayala-Valenzuela, J. N. Mitchell, R. D. McDonald, C. H. Mielke, F. Ronning, J.-C. Griveau, E. Colineau, R. Eloirdi, R. Caciuffo, B. L. Scott, O. Janka, S. M. Kauzlarich, and J. D. Thompson, Localized 5f electrons in superconducting PuCoIn ${}_{\textrm{5}}$ : consequences for superconductivity in PuCoGa ${}_{\textrm{5}}$, Journal of Physics: Condensed Matter 24, 052206 (2012b).
  * Petrovic  _et al._ [2001] C. Petrovic, P. G. Pagliuso, M. F. Hundley, R. Movshovich, J. L. Sarrao, J. D. Thompson, Z. Fisk, and P. Monthoux, Heavy-fermion superconductivity in CeCoIn ${}_{\textrm{5}}$ at 2.3 K, Journal of Physics: Condensed Matter 13, L337 (2001).
  * D. Thompson and Fisk [2012] J. D. Thompson and Z. Fisk, Progress in Heavy-Fermion Superconductivity: Ce115 and Related Materials, Journal of the Physical Society of Japan 81, 011002 (2012).
  * Sechovsky and Havela [1988] V. Sechovsky and L. Havela, Chapter 4 Intermetallic compounds of actinides, in  _Handbook of Ferromagnetic Materials_, Vol. 4 (Elsevier, 1988) pp. 309–491.
  * Gofryk  _et al._ [2009] K. Gofryk, J.-C. Griveau, E. Colineau, R. Jardin, J. Rebizant, F. Wastin, and R. Caciuffo, Magnetic and related properties of PuPdSn, Journal of Nuclear Materials 385, 46 (2009).
  * Söderlind  _et al._ [2019] P. Söderlind, A. Landa, and B. Sadigh, Density-functional theory for plutonium, Advances in Physics 68, 1 (2019).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Physical Review Letters 77, 3865 (1996).
  * Methfessel and Paxton [1989] M. Methfessel and A. T. Paxton, High-precision sampling for Brillouin-zone integration in metals, Physical Review B 40, 3616 (1989).
  * Fredeman  _et al._ [2011] D. J. Fredeman, P. H. Tobash, M. A. Torrez, J. D. Thompson, E. D. Bauer, F. Ronning, W. W. Tipton, S. P. Rudin, and R. G. Hennig, Computationally driven experimental discovery of the CeIr 4 In compound, Physical Review B 83, 224102 (2011).
  * Bauer  _et al._ [2011] E. D. Bauer, P. H. Tobash, J. N. Mitchell, J. A. Kennison, F. Ronning, B. L. Scott, and J. D. Thompson, Magnetic order in Pu ${}_{\textrm{2}}$ M ${}_{\textrm{3}}$ Si ${}_{\textrm{5}}$ (M = Co, Ni), Journal of Physics: Condensed Matter 23, 094223 (2011).
  * Shorikov  _et al._ [2005] A. O. Shorikov, A. V. Lukoyanov, M. A. Korotin, and V. I. Anisimov, Magnetic state and electronic structure of the $\alpha$ and $\delta$ phases of metallic Pu and its compounds, Physical Review B 72, 024458 (2005).
  * Boulet  _et al._ [2005] P. Boulet, E. Colineau, F. Wastin, P. Javorský, J. C. Griveau, J. Rebizant, G. R. Stewart, and E. D. Bauer, Magnetic properties of the two allotropic phases of Pu Ga 3, Physical Review B 72, 064438 (2005).
  * Chudo  _et al._ [2014] H. Chudo, G. Koutroulakis, H. Yasuoka, E. D. Bauer, P. H. Tobash, J. N. Mitchell, and J. D. Thompson, Weak itinerant antiferromagnetism in PuIn ${}_{\textrm{3}}$ explored using ${}^{\textrm{115}}$ In nuclear quadrupole resonance, Journal of Physics: Condensed Matter 26, 036001 (2014).
  * Stewart and Elliott [1985] G. R. Stewart and R. O. Elliott, Anomalous behavior of Pu Al 2, Physical Review B 31, 4669 (1985).
  * Besnus  _et al._ [1992] M. Besnus, A. Braghta, N. Hamdaoui, and A. Meyer, A correlation between specific heat and the ratioTK/TN in magnetic Kondo lattices, Journal of Magnetism and Magnetic Materials 104-107, 1385 (1992).
  * Javorský  _et al._ [1999] P. Javorský, L. Havela, F. Wastin, J. Griveau, E. Bednarczyk, J. Rebizant, J. Sanchez, and P. Vulliet, Magnetic and transport properties of NpNiSn, Journal of Alloys and Compounds 283, 16 (1999).
  * Alsmadi  _et al._ [2011] A. M. Alsmadi, M. S. Torikachvili, M. K. Kothapalli, F. Nasreen, and H. Nakotte, Magnetic phase transitions in CePtSn under high pressure, Journal of Physics: Conference Series 273, 012025 (2011).
  * Divis  _et al._ [1994] M. Divis, H. Nakotte, F. R. D. Boer, P. F. D. Chatel, and V. Sechovsky, Determination of the crystal-field parameters for CePtSn, Journal of Physics: Condensed Matter 6, 6895 (1994).
  * Zhang  _et al._ [2013] J. J. Zhang, G. J. Zhao, and X. X. Liang, First-Principle Studies of Phonons III-N Compound Semiconductors in Wurtzite Structure, International Journal of Applied Physics and Mathematics , 227 (2013).
  * Shrestha  _et al._ [2018] K. Shrestha, D. Antonio, J.-C. Griveau, K. Prokeš, P. Gaczyński, E. Colineau, R. Caciuffo, and K. Gofryk, Magnetic and electronic properties of NpPdSn, Physical Review Materials 2, 074401 (2018).
  * Stewart  _et al._ [1984] G. R. Stewart, Z. Fisk, J. L. Smith, J. O. Willis, and M. S. Wire, New heavy-fermion system, NpBe 13 , with a comparison to UBe 13 and PuBe 13, Physical Review B 30, 1249 (1984).
  * Sarrao  _et al._ [2015] J. Sarrao, E. Bauer, J. Mitchell, P. Tobash, and J. Thompson, Superconductivity in plutonium compounds, Physica C: Superconductivity and its Applications 514, 184 (2015).
  * Bachmann  _et al._ [2019] M. D. Bachmann, G. M. Ferguson, F. Theuss, T. Meng, C. Putzke, T. Helm, K. R. Shirer, Y.-S. Li, K. A. Modic, M. Nicklas, M. König, D. Low, S. Ghosh, A. P. Mackenzie, F. Arnold, E. Hassinger, R. D. McDonald, L. E. Winter, E. D. Bauer, F. Ronning, B. J. Ramshaw, K. C. Nowack, and P. J. W. Moll, Spatial control of heavy-fermion superconductivity in CeIrIn ${}_{\textrm{5}}$, Science 366, 221 (2019).
  * Moore and Van Der Laan [2009] K. T. Moore and G. Van Der Laan, Nature of the 5 f states in actinide metals, Reviews of Modern Physics 81, 235 (2009).
  * Le  _et al._ [2008] M. Le, K. McEwen, F. Wastin, E. Colineau, R. Jardin, and J. Rebizant, Magnetisation, susceptibility and heat capacity studies of, Physica B: Condensed Matter 403, 1035 (2008).
  * Brodsky and Friddle [1975] M. B. Brodsky and R. J. Friddle, Magnetic susceptibility measurements on Actinide-Be13 intermetallic compounds, in  _AIP Conference Proceedings_, Vol. 24 (AIP, 1975) pp. 353–354, iSSN: 0094243X.
  * Nishigori  _et al._ [1996] S. Nishigori, H. Goshima, T. Suzuki, T. Fujita, G. Nakamoto, H. Tanaka, T. Takabatake, and H. Fujii, Pseudogap due to Coherence Kondo Effect in CeNiSn and CeRhSb, Journal of the Physical Society of Japan 65, 2614 (1996).
  * Palstra  _et al._ [1987] T. Palstra, G. Nieuwenhuys, R. Vlastuin, J. Van Den Berg, J. Mydosh, and K. Buschow, Magnetic and electrical properties of several equiatomic ternary U-compounds, Journal of Magnetism and Magnetic Materials 67, 331 (1987).
  * Sechovsky and Havela [1998] V. Sechovsky and L. Havela, Chapter 1 Magnetism of ternary intermetallic compounds of uranium, in  _Handbook of Magnetic Materials_, Handbook of Magnetic Materials, Vol. 11 (Elsevier, 1998) pp. 1–289, iSSN: 1567-2719.
  * Nellis and Brodsky [1974] W. J. Nellis and M. B. Brodsky, Magnetic Properties of Actinide Compounds, in  _The Actinides: Electronic Structure and Related Properties_, Vol. 1, edited by A. J. Freeman and J. B. Darby (Academic Press, Amsterdam, 1974) pp. 265–288.
  * Hecker [2008] S. S. Hecker, Plutonium—an Element Never at Equilibrium, Metallurgical and Materials Transactions A 39, 1585 (2008).


# QGP Physics from Attractor Perturbations

Xin An<EMAIL_ADDRESS>National Centre for Nuclear Research, 02-093
Warsaw, Poland  Michał Spaliński<EMAIL_ADDRESS>National Centre
for Nuclear Research, 02-093 Warsaw, Poland  Physics Department, University of
Białystok, 15-245 Białystok, Poland

###### Abstract

The strong longitudinal expansion characteristic of heavy-ion collisions leads
to universal attractor behaviour of the resulting drop of Quark-Gluon Plasma
already at very early times. Assuming approximate boost invariance, we
incorporate transverse dynamics of this system by linearizing the Mueller-
Israel-Stewart theory around the attractor. The result is a system of coupled
ordinary differential equations which describe the proper-time evolution of
Fourier modes encoding the transverse structure of the initial energy
deposition. The late time asymptotic behaviour of solutions is described by
transseries which make manifest the stability of the attractor against
transverse perturbations. In this framework, most of the physically relevant
information resides in the exponentially suppressed corrections to evolution
along the attractor, which are not yet negligible at freeze-out. These
findings also suggest a simple numerical approach to QGP dynamics which
accounts for the transverse dynamics using a finite number of Fourier modes.
We show that this approach is able to describe collectivity at the level of
the transverse anisotropy, as a surrogate for elliptic flow. Physical
observables can be expressed in terms of the asymptotic data evaluated at
freeze-out, which we illustrate by calculating the final multiplicity
distributions.

##  I Introduction

Quark-Gluon Plasma (QGP) is created in heavy-ion collision experiments in
highly anisotropic, nonequilibrium states. Many features of the subsequent
evolution are successfully described by models formulated in the language of
fluid dynamics, which are applied long before local equilibrium is
established. This implies a vast reduction in the number of degrees of freedom
at the earliest moments following the collision. A possible explanation of
this follows from a key kinematical feature of heavy-ion collisions: the
dominant longitudinal expansion at the prehydrodynamic stages of evolution. An
idealisation of this situation assumes boost-invariance along the collision
axis and neglects the transverse dynamics [1]. It has been shown in a number
of models that in such circumstances an early-time, far-from-equilibrium
attractor governs the dynamics until the QGP drop approaches a state amenable
to a hydrodynamic description with small gradients [2, 3, 4, 5]. Within such a
picture, the information about the initial state is contained in a single
scale which characterises the particular attractor background, up to
corrections which are exponentially suppressed — and thus effectively lost —
at asymptotically late times.

Of course, what actually happens in heavy-ion collision experiments is that
the system does not survive until such asymptotically late times, because as
the effective local temperature drops, the cooling QGP is converted into a
stream of hadrons which source the multitude of particles registered in
detectors. The result of this is a wealth of information which reflects the
structure of the initial state in the plane transverse to the collision axis.
Almost all of the information inferred from such measurements resides in the
exponentially suppressed terms which are not yet negligibly small at the time
of freeze-out. The precise form of these terms is a reflection of the
nonhydrodynamic modes present in the underlying microscopic theory.

In this paper we aim to make this description explicit and generalize it to
account for the transverse dynamics which is a crucial element of many
observed phenomena. The goal is to extend the Bjorken model in a way which
still allows for some analytic insights. Our analysis retains the assumption
of longitudinal boost-invariance, and accounts for the transverse dynamics at
the level of linearization around the boost-invariant and transversely-
homogeneous attractor. The resulting model consists of a system of coupled
ordinary differential equations for a set of modes parameterised by the
transverse wave vector. We establish the form of the exponential corrections
to the background and show that essentially all physical observables are
determined be these terms. Thus, the physics of QGP flow in heavy-ion
collisions provides an example of a very nontrivial dynamical system where the
exponentially damped corrections to asymptotic results are not just non-
negligible, but contain almost all of the physically relevant information. As
a bi-product, we arrive at a very simple framework which can easily be studied
numerically and captures effects of the transverse dynamics, such as elliptic
flow.

Focusing on a version of the hydrodynamic model originallty due to Müller [6]
and Israel and Stewart [7, 8] (MIS), we begin by studying this dynamical
system numerically, and then apply asymptotic methods to describe the
solutions analytically at late times. This model can be viewed as a “UV-
regularised” form of Navier-Stokes (NS) theory. It describes not just the
universal hydrodynamic regime, but also a nonhydrodynamic sector which is
essential for causality and stability and provides a very simple model of the
dynamics at earlier times (see e.g. [9, 10, 11]).

The basic degrees of freedom in the hydrodynamic description are the classical
fields $\mathcal{E}$ (the energy density), $u^{\mu}$ (the flow velocity) and
${\pi^{\mu\nu}}$ (the shear-stress tensor). They satisfy the following set of
partial differential equations

$\displaystyle u\cdot\nabla\mathcal{E}$
$\displaystyle=-(\mathcal{E}+p)\nabla\cdot u+u^{\nu}\nabla^{\mu}\pi_{\mu\nu},$
(1a) $\displaystyle(\mathcal{E}+p)u\cdot\nabla u_{\mu}$
$\displaystyle=-\Delta_{\mu\nu}\nabla^{\nu}p-\Delta_{\mu\nu}\nabla_{\lambda}\pi^{\nu\lambda},$
(1b)
$\displaystyle\Delta_{\mu\alpha}\Delta_{\nu\beta}u\cdot\nabla\pi^{\alpha\beta}$
$\displaystyle=-\left(1+\frac{4}{3}\tau_{\pi}\nabla\cdot
u\right)\pi_{\mu\nu}-2\eta\sigma_{\mu\nu},$ (1c)

where $\nabla_{\mu}$ is the covariant derivative, $\Delta_{\mu\nu}\equiv
g_{\mu\nu}+u_{\mu}u_{\nu}$ is the transverse projector,
$\sigma_{\mu\nu}=\frac{1}{2}\Delta_{\mu\alpha}\Delta_{\nu\beta}(\nabla^{\alpha}u^{\beta}+\nabla^{\beta}u^{\alpha}-\frac{2}{3}\Delta^{\alpha\beta}\nabla\cdot
u)$ is the shear tensor, $\eta$ is the shear viscosity and $\tau_{\pi}$ is the
relaxation time for $\pi_{\mu\nu}$. Throughout this paper we assume an
equation of state and transport coefficients dictated by conformal invariance:

$\displaystyle\mathcal{E}=\frac{1}{3}p=C_{e}T^{4},\quad\eta=\frac{4}{3}C_{e}C_{\eta}T^{3},\quad\tau_{\pi}=C_{\tau}T^{-1},$
(2)

where $C_{e}$, $C_{\eta}$, and $C_{\tau}$ are dimensionless, constant
transport coefficients (see e.g. [10]) and $T$ is the effective temperature.
For $\mathcal{N}=4$ supersymmetric Yang-Mills theory we have
$C_{e}=8\pi^{2}/15,C_{\eta}=1/4\pi,C_{\tau}=(2-\ln 2)/2\pi$. These values
often serve as a point of reference, and we have adopted them in our numerical
calculations.

##  II The attractor background

The basic physical picture we adopt is that of the Bjorken model [1]: at
sufficiently high energies, in the first approximation, the system exhibits
boost invariance in the longitudinal direction and homogeneity in the
transverse directions (perpendicular to the collision axis $z$). We will refer
to this approximate description as the background; perturbations dependent on
the transverse coordinates will subsequently be treated at the linearized
level. Under these assumptions, the energy-momentum tensor of the conformal
MIS theory can be parameterised in terms of only two functions of proper time
$\tau=\sqrt{t^{2}-z^{2}}$: the effective temperature $T(\tau)$, and the
pressure anisotropy $\mathcal{A}(\tau)\equiv 9\pi^{i}_{i}/2\mathcal{E}$, where
$i=1,2$ labels the transverse coordinates (for details see [11]). Eqs. (1)
then reduce to

$\displaystyle\tau\partial_{\tau}\ln T(\tau)$
$\displaystyle=-\frac{1}{3}+\frac{1}{18}\mathcal{A}(\tau),$ (3a)
$\displaystyle\tau\partial_{\tau}\mathcal{A}(\tau)$
$\displaystyle=8\alpha^{2}-\frac{\tau
T(\tau)}{C_{\tau}}\mathcal{A}(\tau)-\frac{2}{9}\mathcal{A}^{2}(\tau),$ (3b)

where $\alpha^{2}\equiv C_{\eta}/C_{\tau}$. In the NS limit ($C_{\tau}\to 0$),
Eq. (3b) becomes algebraic and one finds the well-known solutions

$\displaystyle T(\tau)$
$\displaystyle=\Lambda(\Lambda\tau)^{-1/3}\left(1-\frac{2C_{\eta}}{3}(\Lambda\tau)^{-2/3}\right),$
(4a) $\displaystyle\mathcal{A}(\tau)$ $\displaystyle=\frac{8C_{\eta}}{\tau
T(\tau)}\sim
8C_{\eta}(\Lambda\tau)^{-2/3}\left(1+\sum_{n=1}^{\infty}\left(\frac{2C_{\eta}}{3}\right)^{n}(\Lambda\tau)^{-2n/3}\right),$
(4b)

where $\Lambda$ is a integration constant with the dimension of energy. The
general solutions of Eqs. (3) are not available in closed form and can only be
found in certain regimes. However, there are some special (and exact)
solutions which are completely independent of any integration constants:
$T=0,\mathcal{A}=\pm 6\alpha$ and
$T=\frac{2(4-\alpha^{2})}{3C_{\tau}\tau},\mathcal{A}=-12$. We regard these
solutions as unphysical, since they requires fine-tuning of the initial
conditions. To understand physically interesting solutions we will resort to
approximate solutions valid at small or large proper times.

At early times Eqs. (3) admits two special families of solutions which can be
characterised by having a finite value of the pressure anisotropy at $\tau=0$:

$\displaystyle T_{\pm}(\tau)$
$\displaystyle\sim\mu(\mu\tau)^{-\frac{1}{3}(1\mp\alpha)}\left(1\pm\sum_{n=1}^{\infty}t_{n}^{(0)}(\mu\tau)^{\frac{n}{3}(2\pm\alpha)}\right),$
(5a) $\displaystyle\mathcal{A}_{\pm}(\tau)$ $\displaystyle\sim\pm
6\alpha\left(1\pm\sum_{n=1}^{\infty}a_{n}^{(0)}(\mu\tau)^{\frac{n}{3}(2\pm\alpha)}\right),$
(5b)

where $\mu$ is an integration constant defining the initial condition. The
early-time series coefficients $t_{n}^{(0)}$ and $a_{n}^{(0)}$ are rational
functions of the transport coefficients and the first few are given by

$\displaystyle
t^{(0)}_{1}=-\frac{3\alpha}{C_{\tau}(2+\alpha)(2+9\alpha)},\quad
t^{(0)}_{2}=\frac{9\alpha(2+9\alpha+12\alpha^{2})}{2C_{\tau}^{2}(2+\alpha)^{2}(2+5\alpha)(2+9\alpha)^{2}},\quad\cdots$
(6a) $\displaystyle a^{(0)}_{1}=\frac{18\alpha}{C_{\tau}(2+9\alpha)},\quad
a^{(0)}_{2}=\frac{54\alpha(2+7\alpha+7\alpha^{2})}{C_{\tau}^{2}(2+\alpha)(2+5\alpha)(2+9\alpha)^{2}},\quad\cdots$
(6b)

The power series appearing in Eqs. (5) have a finite radius of convergence.
The upper sign in Eqs. (5) defines the class of attractor solutions labelled
by $\mu$. It is easy to check that all these solutions are mapped to the
universal attractor introduced in Ref. [2] (see also the reviews [10, 12, 13,
11]). The solutions with the lower sign are mapped to the “repulsor” solution
noted in the original approach of Ref. [2] (see also the recent Ref. [14]).

There are also “generic” solutions, characterised by a pressure anisotropy
which diverges at early times. These solutions approach the attractor already
in the far from equilibrium regime. To understand this, it is important to
keep in mind that since Eqs. (3) are non-autonomous, the full phase space of
solutions is three-dimensional and can be naturally parameterised by
$(\tau,T,\tau T^{\prime})$ [15] (see also Refs. [13, 11]). The attractor is a
two-dimensional surface in this full phase space. On any constant proper-time
slice of the phase space, $\mu$ is a parameter which labels points along the
attractor curve on that slice. The generic solutions viewed on a sequence of
such phase-space slices at increasing values of $\tau$ gravitate toward this
locus, as a consequence of fast longitudinal expansion at early times. This is
illustrated in Fig. 1.

Figure 1: The left panel shows the early-time attractor in the 3D phase space
$(\tau,T,\tau T^{\prime})$ with uniformly distributed initial conditions
(points) at $\tau=0.1$ fm. The 2D attractor surface is obtained using Padé
approximation of the series in Eqs. (5). The right panel are the snapshots of
the 3D plot at different times, where the solid and dashed line represents the
series solution and its Padé approximant respectively.

At late times, we expect the system to approach equilibrium in a way
consistent with the original insights of Bjorken [1]. In that regime,
solutions to Eqs. (3) can be represented in the form of transseries:

$\displaystyle T(\tau)$
$\displaystyle\sim\Lambda\left[(\Lambda\tau)^{-\frac{1}{3}}\left(1+\sum_{n=1}^{\infty}t_{n}^{(\infty)}(\Lambda\tau)^{-\frac{2n}{3}}\right)+C_{\infty}(\Lambda\tau)^{-\frac{2}{3}(1-\alpha^{2})}e^{-\frac{3}{2C_{\tau}}(\Lambda\tau)^{2/3}}\left(1+\mathcal{O}((\Lambda\tau)^{-\frac{2}{3}})\right)+\cdots\right],$
(7a) $\displaystyle\mathcal{A}(\tau)$ $\displaystyle\sim
8C_{\eta}(\Lambda\tau)^{-\frac{2}{3}}\left(1+\sum_{n=1}^{\infty}a_{n}^{(\infty)}(\Lambda\tau)^{-\frac{2n}{3}}\right)+C_{\infty}^{\prime}(\Lambda\tau)^{-\frac{1}{3}(1-3\alpha^{2})}e^{-\frac{3}{2C_{\tau}}(\Lambda\tau)^{2/3}}\left(1+\mathcal{O}((\Lambda\tau)^{-\frac{2}{3}})\right)+\cdots$
(7b)

These are of the form of an asymptotic power series augmented by an infinite
set of exponential transseries contributions, of which only the leading one is
displayed above. The initial conditions are mapped to the dimensionful scale
$\Lambda$ and the dimensionless integration constants $C_{\infty}$ and
$C^{\prime}_{\infty}$. The first few coefficients of the late-time series
$t_{n}^{(\infty)}$ and $a_{n}^{(\infty)}$ are given by

$\displaystyle t^{(\infty)}_{1}=-\frac{2C_{\eta}}{3},\quad
t^{(\infty)}_{2}=-\frac{2C_{\eta}C_{\tau}}{9},\quad\cdots$ (8a) $\displaystyle
a^{(\infty)}_{1}=\frac{2C_{\tau}}{3}(1+\alpha^{2}),\quad
a^{(\infty)}_{2}=\frac{2C_{\tau}^{2}}{9}(4-\alpha^{2}+2\alpha^{4}),\quad\cdots$
(8b)

Note that in the limit $C_{\tau}\to 0$ the leading and next-to-leading terms
in Eqs. (7) reduce to their NS form, given in Eq. (4), upon substituting the
coefficients in Eqs. (8). The power series appearing in Eqs. (7) have a
vanishing radius of convergence and are best interpreted in the sense of
asymptotic analysis — through optimal truncation or by Borel summation (see
e.g. Refs. [16, 14]).

The description of the attractor presented above differs from the original
formulation of Ref. [2] in that it uses proper time as the evolution
parameter, rather than the dimensionless evolution variable $w\sim\tau T$.
This is natural when the dynamic system possesses scales in addition to
temperature, as will be discussed in Sec. III.

##  III The transverse perturbations

In the previous Section, we discussed an idealized description which was
homogeneous in the transverse plane. While this idealization provides a useful
first approximation for modeling the heavy-ion collisions at sufficiently high
energies, it cannot account for physical observables which depend on the
structure of the plasma in the transverse plane. In this Section, we relax
this condition by considering additional fields which aim to model the
transverse dynamics. These additional fields arise a perturbations of the
fully-nonlinear hydrodynamic equations Eq. (1) around the attractor
background.

###  III.1 The linearized equations

We will look for solutions that can be approximated by the boost-invariant and
translation-invariant background solution discussed in Sec. II and a
perturbation depending also on the transverse coordinates:

$\displaystyle T(\tau,\bm{x})=T(\tau)+\delta
T(\tau,\bm{x}),\quad{u^{\mu}}(\tau,\bm{x})=u^{\mu}+\delta{u^{\mu}}(\tau,\bm{x}),\quad{\pi^{\mu\nu}}(\tau,\bm{x})=\pi^{\mu\nu}(\tau)+\delta{\pi^{\mu\nu}}(\tau,\bm{x}),$
(9)

where $\bm{x}=(x_{1},x_{2})$ labels the coordinates of the transverse plane.
We shall always retain the argument $\bm{x}$ of the above quantities to
distinguish the full ones from the background, which depends only on $\tau$
(and, as the only argument, is often suppressed). The background fields $T$
and $\pi^{\mu\nu}$ are taken to be on the attractor locus defined in Sec. II,
and $u^{\mu}=(-1,\bm{0},0)$. Due to the assumed symmetries and the
transversality condition $u^{\mu}\pi_{\mu\nu}=0$, the shear-stress tensor has
only one independent background component (i.e., $\pi^{11}=\pi^{22}$) and
three independent perturbation components (i.e.,
$\delta\pi^{11},\delta\pi^{22},\delta\pi^{12}$) that couple to the
perturbation of hydrodynamic fields $\delta T$ and $\delta\bm{u}$. We also put
$\delta\pi^{i\eta}=0$, which is consistent due to $\delta u^{\eta}=0$.

The perturbation fields can also be normalized by the background energy scale
$T(\tau)$, i.e.,

$\displaystyle\delta\hat{T}(\tau,\bm{x})=\frac{\delta
T(\tau,\bm{x})}{T(\tau)},\quad\delta\hat{\pi}_{ij}(\tau,\bm{x})=\frac{\delta\pi_{ij}(\tau,\bm{x})}{C_{e}T(\tau)^{4}},$
(10)

such that all six perturbation fields, denoted collectively by
$\hat{\phi}(\tau,\bm{x})=(\hat{\delta}T,\delta u_{1},\delta
u_{2},\hat{\delta}\pi_{12},\hat{\delta}\pi_{11},\hat{\delta}\pi_{22})$, are
dimensionless. Since the background is independent of the transverse
coordinates, it is also natural and convenient to introduce the Fourier
transforms just for the perturbations:

$\displaystyle\hat{\phi}(\tau,\bm{x})=\int\frac{d^{2}k}{(2\pi)^{2}}\
e^{i\bm{k}\cdot\bm{x}}{\hat{\phi}}(\tau,\bm{k}),$ (11)

where we retain the argument $\bm{x}$ or $\bm{k}$ to distinguish perturbation
fields $\hat{\phi}$ in different Fourier spaces.

Figure 2: Evolution of Fourier modes with various values of $k$ between
$0-3\,\text{fm}^{-1}$, with the background on the attractor,
$\mathcal{A}(\tau_{i})=6\alpha$ (left) and off it,
$\mathcal{A}(\tau_{i})=30\alpha$ (right). The initial time is $\tau_{i}=0.01$
fm.

Linearization of the full MIS equations (1) around an attractor solution leads
to a system of six linear partial differential equations for the
perturbations. For each value of $\bm{k}$ the set of six modes
$\hat{\phi}(\tau,\bm{k})\equiv(\delta\hat{T},\delta u_{1},\delta
u_{2},\delta\hat{\pi}_{12},\delta\hat{\pi}_{11},\delta\hat{\pi}_{22})$
satisfies a linear system of evolution equations which can be explicitly
written as

$\displaystyle
4\left(\tau\partial_{\tau}+\frac{2\mathcal{A}}{9}\right)\delta\hat{T}+\frac{i(12+\mathcal{A})}{9}\tau
k_{i}\delta u^{i}-\delta\hat{\pi}^{i}_{i}$ $\displaystyle=0,$ (12a)
$\displaystyle\frac{i4}{3}\tau
k_{i}\delta\hat{T}+\frac{1}{9}\left[(12+\mathcal{A})\tau\partial_{\tau}-4+\frac{8C_{\eta}}{C_{\tau}}+\left(\frac{7}{3}-\frac{\tau
T}{C_{\tau}}\right)\mathcal{A}\right]\delta u_{i}+i\tau
k^{j}\delta\hat{\pi}_{ij}$ $\displaystyle=0,$ (12b)
$\displaystyle\frac{1}{9}\left(\frac{\tau
T}{C_{\tau}}\mathcal{A}-\frac{32C_{\eta}}{C_{\tau}}\right)\delta_{ij}\delta\hat{T}+\frac{i4\tau}{27}\left[\left(\mathcal{A}-\frac{6C_{\eta}}{C_{\tau}}\right)\delta_{ij}k_{\ell}+\frac{9C_{\eta}}{C_{\tau}}(\delta_{i\ell}k_{j}+\delta_{j\ell}k_{i})\right]\delta
u^{\ell}+\left(\tau\partial_{\tau}+\frac{\tau
T}{C_{\tau}}+\frac{2\mathcal{A}}{9}\right)\delta\hat{\pi}_{ij}$
$\displaystyle=0.$ (12c)

Since modes with different $\bm{k}$ are decoupled at the linearized level, we
omit the arguments $\tau$ and $\bm{k}$ unless needed.

It is straightforward to solve the system given in Eqs. (12) together with
Eqs. (3) numerically for any given wave vector $k$. The most obvious
consequence of such a calculation is the fact that modes with higher $k$ are
damped more strongly than modes with low $k$, as illustrated in the left panel
of Fig. 2. The origin of this phenomenon can be understood analytically, as
will be discussed in Sec. III.2. Moreover, the perturbations starting away
from the attractor are damped more significantly compared to perturbations
starting on the attractor, as illustrated in the right panel of Fig. 2. This
demonstrates the stability of the linearization around the attractor
background.

###  III.2 Late-time behaviour of perturbations

In this Section we will try to understand analytically the late proper-time
behaviour of the modes $\hat{\phi}$. The set of six ODEs describing the
evolution of these modes can be analysed by standard asymptotic methods,
although the complexity of this problem makes the calculation technically
nontrivial.

For the purpose of calculating the asymptotics it is convenient to introduce
the transverse divergence $\delta\theta(\tau,\bm{x})\equiv\partial_{i}\delta
u^{i}$ and longitudinal vorticity
$\delta\omega(\tau,\bm{x})\equiv\epsilon_{ij}\partial^{i}\delta u^{j}$ where
$\epsilon_{ij}$ is the Levi-Civita symbol. In $\bm{k}$ space we will use the
following dimensionless quantities, normalised by $k\equiv|\bm{k}|$:

$\displaystyle\delta\hat{\theta}(\tau,\bm{k})=\frac{\delta{\theta}(\tau,\bm{k})}{k}=i\hat{k}_{i}\delta
u^{i}(\tau,\bm{k}),\quad\delta\hat{\omega}(\tau,\bm{k})=\frac{\delta{\omega}(\tau,\bm{k})}{k}=i\epsilon_{ij}\hat{k}^{i}\delta
u^{j}(\tau,\bm{k}),$ (13)

where $\hat{k}_{i}\equiv k_{i}/k$. For such a choice of variables, the linear
system of evolution equations reads

$\displaystyle
4\left(\tau\partial_{\tau}+\frac{2\mathcal{A}}{9}\right)\delta\hat{T}+\frac{(12+\mathcal{A})}{9}\tau
k\delta\hat{\theta}-\delta\hat{\pi}^{i}_{i}$ $\displaystyle=0,$ (14a)
$\displaystyle-\frac{4}{3}k^{2}\delta\hat{T}+\frac{1}{9\tau}\left[(12+\mathcal{A})\tau\partial_{\tau}-4+\frac{8C_{\eta}}{C_{\tau}}+\left(\frac{7}{3}-\frac{\tau
T}{C_{\tau}}\right)\mathcal{A}\right]k\delta\hat{\theta}-k_{i}k_{j}\delta\hat{\pi}^{ij}$
$\displaystyle=0,$ (14b)
$\displaystyle\frac{1}{9\tau}\left[(12+\mathcal{A})\tau\partial_{\tau}-4+\frac{8C_{\eta}}{C_{\tau}}+\left(\frac{7}{3}-\frac{\tau
T}{C_{\tau}}\right)\mathcal{A}\right]k\delta\hat{\omega}-\epsilon_{ij}k^{j}k^{\ell}\delta\hat{\pi}^{i}_{\ell}$
$\displaystyle=0,$ (14c) $\displaystyle\frac{1}{9}\left(\frac{\tau
T}{C_{\tau}}\mathcal{A}-\frac{32C_{\eta}}{C_{\tau}}\right)\delta_{ij}\delta\hat{T}+\frac{4\tau}{27k}\left[\frac{18C_{\eta}}{C_{\tau}}k_{i}k_{j}+\left(\mathcal{A}-\frac{6C_{\eta}}{C_{\tau}}\right)k^{2}\delta_{ij}\right]\delta\hat{\theta}$
$\displaystyle-\frac{4C_{\eta}(k_{i}\epsilon_{j\ell}+k_{j}\epsilon_{i\ell})k^{\ell}\tau}{3C_{\tau}k}\delta\hat{\omega}+\left(\tau\partial_{\tau}+\frac{\tau
T}{C_{\tau}}+\frac{2\mathcal{A}}{9}\right)\delta\hat{\pi}_{ij}$
$\displaystyle=0.$ (14d)

At this juncture, one way to proceed is to rewrite this system in terms of
higher-order ODEs; remarkably, it can be written as a set of three second-
order ODEs for $\delta\hat{T},\delta\hat{\theta},\delta\hat{\omega}$ — they
are given in Appendix A. Two of these equations couple $\delta\hat{T}$ and
$\delta\hat{\theta}$, while the third involves $\delta\hat{\omega}$ alone. The
two coupled equations can be combined into a fourth order ODE for
$\delta\hat{T}$, whose asymptotic behaviour can be studied by standard methods
(see e.g. Ref [17]). This is in principle straightforward, but is somewhat
challenging in practice due to the complexity of the coefficients which appear
in this equation. Alternatively, one can analyse the system of six first order
ODEs directly using the methods developed in Ref. [18]. Both approaches gives
rise to the same late-time asymptotic solutions whose leading terms take the
form:

$\displaystyle\delta\hat{T}$
$\displaystyle=\sum_{i=1}^{4}C_{i}(\Lambda\tau)^{\beta_{i}}e^{-i\omega_{i}\tau-
A_{i}(\Lambda\tau)^{2/3}},$ (15a) $\displaystyle\delta\hat{\theta}$
$\displaystyle=\sum_{i=1}^{4}C^{\prime}_{i}(\Lambda\tau)^{\beta^{\prime}_{i}}e^{-i\omega_{i}\tau-
A_{i}(\Lambda\tau)^{2/3}},$ (15b) $\displaystyle\delta\hat{\omega}$
$\displaystyle=\sum_{i=5}^{6}C_{i}(\Lambda\tau)^{\beta_{i}}e^{-i\omega_{i}\tau-
A_{i}(\Lambda\tau)^{2/3}}.$ (15c)

The initial conditions are accounted for by the amplitudes $C_{1},\dots,C_{6}$
and $C_{1}^{\prime},\dots,C_{4}^{\prime}$. The primed integration constants
are related to the unprimed ones by the relations

$\displaystyle C_{1}^{\prime}=3ic_{\alpha}C_{1},\quad
C_{2}^{\prime}=-3ic_{\alpha}C_{2},\quad
C_{3}^{\prime}=-\frac{(1-3\alpha^{2})\Lambda}{2\alpha^{2}k}C_{3},\quad
C_{4}^{\prime}=\frac{\Lambda}{c_{\alpha}^{2}C_{\tau}k}C_{4},$ (16)

so that only the unprimed integration constants $C_{1},\dots,C_{6}$ are
independent. Recall also that each of the perturbations appearing above
depends on the wave vector $\bm{k}$, and so do the coefficients
$C_{1},\dots,C_{6}$.

Figure 3: These figures show the analytic late-time asymptotic solution (solid
line) and the numerical results (discrete points) for various values of wave
vector $k$. They were obtained by fitting the amplitudes $C_{i}$. The inset
plots are provided to give a better view of the small-$\tau$ region.

The remaining quantities appearing in Eqs. (15) above are determined in terms
of the parameters of the theory and do not depend on the initial state. The
parameters appearing in the exponentials are given by

$\displaystyle A_{1}=A_{2}=\frac{\alpha^{2}}{C_{\tau}c_{\alpha}^{2}},\quad
A_{3}=\frac{3}{2C_{\tau}},\quad A_{4}=\frac{1}{2C_{\tau}c_{\alpha}^{2}},\quad
A_{5}=A_{6}=\frac{3}{4C_{\tau}},$
$\displaystyle\omega_{1}=-\omega_{2}=c_{\alpha}k\left[1+\frac{2\alpha^{2}}{3c_{\alpha}^{2}}\left(2C_{\tau}(1-\alpha^{2})-\frac{(1+\alpha^{2})\Lambda^{2}}{C_{\tau}^{2}c_{\alpha}^{4}k^{2}}\right)(\Lambda\tau)^{-2/3}\right],\quad\omega_{3}=\omega_{4}=0,$
$\displaystyle\omega_{5}=-\omega_{6}=\alpha
k\left[1-\left(C_{\tau}\alpha^{2}+\frac{3\Lambda^{2}}{8C_{\tau}^{2}\alpha^{2}k^{2}}\right)(\Lambda\tau)^{-2/3}\right],$
(17)

where $c_{\alpha}=\sqrt{(1+4\alpha^{2})/3}$ is the modified (asymptotic) speed
of sound. The constant coefficients $A_{i}$ are the leading contributions to
eigenvalues of the coupled linear equations describing perturbations around
the attractor. They are real, positive and independent of the wave vectors
$k$, implying stability of the attractor against transverse perturbations
regardless of their scales. Note that the quantities $\omega_{1}=-\omega_{2}$
and $\omega_{5}=-\omega_{6}$ are actually functions of $\tau$, which also have
a nontrivial dependence on the wave vectors $k$. At late times these
quantities approach frequencies of oscillation which become harmonic in that
domain.

The coefficients appearing in the power-law factors in Eqs. (15) are given by

$\displaystyle\beta_{1}=\beta_{2}=\beta^{\prime}_{1}=\beta^{\prime}_{2}=\frac{1}{54c_{\alpha}^{4}}\left(1+8\alpha^{2}+64\alpha^{4}+32\alpha^{6}+\frac{4\alpha^{2}\Lambda^{2}}{C_{\tau}^{3}c_{\alpha}^{4}k^{2}}\right),\quad\beta_{3}=\beta^{\prime}_{3}+1=-\frac{2}{3}(1-\alpha^{2}),$
$\displaystyle\beta_{4}=\beta^{\prime}_{4}+\frac{1}{3}=\frac{2\alpha^{2}}{27c_{\alpha}^{4}}\left(1-16\alpha^{2}-\frac{2\Lambda^{2}}{C_{\tau}^{3}c_{\alpha}^{4}k^{2}}\right),\quad\beta_{5}=\beta_{6}=\frac{1}{6}(1+2\alpha^{2}),$
(18)

A key point which follows from these relations is the suppression of modes
with large $k$ relative to modes with smaller values of $k$. This fact was
already noted in the numerical solutions discussed in Sec. III.1, but the
asymptotic formulae Eqs. (15) in conjunction with Eq. (18) reveal the precise
form of this suppression. Note also that the solution labelled by $i=3$
reproduces the transseries solution describing the homogeneous background
(with the integration constant independent of $k$ in that case).

This highly nontrivial structure reproduces the behaviour seen in numerical
simulations very accurately. At late time, the dominant asymptotic solutions
for $\delta\hat{T}$ and $\delta\hat{\theta}$ are those which are least damped,
i.e. with the smallest value of $A_{i}$. For the $\mathcal{N}=4$
supersymmetric Yang-Mills theory these are terms labeled by $i=1,2$ in the
solutions of Eqs. (15). These formulae then capture only the leading
exponential behaviour, but still give an excellent account of the numerical
solutions, as illustrated in Fig. 3. It is also quite interesting to compare
these results with the asymptotics of Navier-Stokes theory, as well as with
the case of perfect fluids. This is discussed in Appendix B.

To summarise, at late times the characteristics of the initial state are
mapped to the scale $\Lambda$ and dimensionless amplitudes $C_{1}(k)$ —
$C_{6}(k)$ which carry information about the structure of the initial data in
the transverse plane. These numbers can be matched to a given numerical
solution. All physical observables can be expressed in terms of this
asymptotic data, in manner described explicitly in the following Section.

###  III.3 Spacetime evolution of QGP

Figure 4:  The evolution of temperature profile in the transverse plane for
different at different instants of proper time (shown at the bottom of each
image). The initial conditions assume the nuclear radius of $6.62$ fm
(corresponding to lead). The impact parameter was taken to be $b=6$ fm. The
dimensions of the region covered by the images are $60$ fm by $60$ fm.

The initial state resulting from a heavy ion collision can described by
providing values of all the fields at the initial proper time as functions of
the coordinates in the transverse plane. This initial data can be equivalently
represented in terms of the Fourier modes $\hat{\phi}$ at the initial time.
Realistic initial conditions are superpositions of infinitely many Fourier
modes, but due to the strong damping of modes with high $k$, one can
approximate the initial state by neglecting modes corresponding to large wave
vectors. In practice, one needs to evaluate the initial conditions on a finite
grid in the transverse plane, calculate the Fourier transform, evolve the
modes by solving Eqs. (12) and then compute the inverse Fourier transform to
reconstruct the spacetime history of the QGP. The problem can thus be
described by a finite set of modes which encode the transverse structure of
the initial state. The procedure can be made highly efficient by applying the
Fast Fourier Transform (see e.g. Ref. [19]).

In our pilot study we consider a computational domain which is a $60$ fm by
$60$ fm square region in the transverse plane, described by a regular grid of
$100$ by $100$ points. A corresponding (conjugate) grid is then constructed in
$\bm{k}$-space — this also introduces a cutoff on large $k$. We have
implemented the steps outlined above for the case of very simple initial
conditions which assume $\delta T$ in the form of a Gaussian distribution in
transverse coordinates with a specified impact parameter $b$. The remaining
fields are taken to vanish at the initial time, which is a natural choice at
least in the case of the transverse velocity perturbation. This leads to
results shown Fig. 4. The data generated this way was also used to obtain some
more quantitative results described in the following Section. One can of
course apply this approach to initial states generated in other ways, such as
events generated using Glauber Monte-Carlo.

##  IV QGP Physics

In this paper we will not undertake any systematic studies aiming to establish
which physical effects are captured by the set of approximations which we
described here. However, in this Section we would like to emphasise that all
such effects are encoded in the exponentially damped contributions appearing
in Eqs. (15). To this end, we will describe some quantitative analysis of the
numerical calculation described in Sec. III.3.

###  IV.1 Transverse anisotropy

Figure 5: The dependence of transverse anisotropy $\mathcal{A}_{T}$ on $\tau$
for different values of $b$ (left) and on $b$ for different values of $\tau$.
The initial conditions for the perturbations are set at $\tau=0.4\mathrm{fm}$
and are otherwise the same as those used to produce Fig. 4.

As a first indication that we can at least qualitatively account for some
nontrivial physics related to the transverse expansion, we would like to
address the issue of collectivity at the level of elliptic flow. Similarly to
the pressure anisotropy which characterizes the difference between the
longitudinal and two transverse stress tensor components, one can introduce
the notion of transverse anisotropy (also referred to the momentum
anisotropy), which characterizes the difference between the two transverse
diagonal components themselves. This anisotropy is nonvanishing when
rotational invariance is broken and is approximately twice as large as the
flow harmonic coefficient $v_{2}$ which serves as an important indicator of
collective flow in heavy-ion collision experiments[20, 21, 22]. The leading
terms in the transverse anisotropy are

$\displaystyle\mathcal{A}_{T}(\tau)\equiv\frac{\langle{T_{11}-T_{22}}\rangle_{\perp}}{\langle{T_{11}+T_{22}}\rangle_{\perp}}=\frac{12\langle{\delta
u_{1}^{2}-\delta
u_{2}^{2}}\rangle_{\perp}+9\langle{\delta\hat{\pi}_{11}-\delta\hat{\pi}_{22}}\rangle_{\perp}}{2(3+\mathcal{A})},$
(19)

where $\langle{\dots}\rangle_{\perp}=\int dx_{1}dx_{2}(\dots)/\int
dx_{1}dx_{2}$ denotes the spacial average over the region occupied by the QGP
in the transverse plane at a given instant of proper time. To calculate this
quantity it is necessary to determine the transverse extent of the fluid.
Since the background is homogeneous in the transverse plane, the boundary of
the plasma drop is defined by the perturbation. Due to the finite resolution
of the calculation, this means in practice that this boundary is taken to be
the locus where the energy density deviates from the value at the periphery of
the computational domain. The outcomes of this computation are illustrated in
Fig. 5. Note also that the finiteness of the freeze-out time $\tau_{f}$
implies a finite spacial extent of the QGP drop.

###  IV.2 Freeze-out

Figure 6: The left panel shows the transverse momentum dependence of pion
distribution without and with perturbations (represented by dashed and solid
lines), while the right panel shows the ratio of the perturbative contribution
and the background.

The numerical solutions discussed in Sec. III show that at asymptotically late
times the modes representing the perturbations are damped away and eventually
the system follows the Bjorken background solution. The QGP never reaches that
stage however due to hadronization: at some time the system described by the
hydrodynamic variables is converted into a stream of hadrons. To calculate
their distribution one needs to capture the fine details of the flow during
the freeze-out epoch, when the hydrodynamic fields are converted into a set of
outgoing particles. Within the linearization framework discussed here it is
reasonable to model freeze-out by defining a single instant of proper time
$\tau_{f}=\tau_{f}(\tau)$ which is independent of the transverse coordinates
and only set by the homogeneous background, the corrections to which are
beyond the linear approximation we consider here.

In heavy-ion collision experiments the observables refer to hadronized
particles whose momentum distribution is given by the well-known Cooper–Frye
formula [23]. For the Bjorken background
$u^{\mu}=(\cosh\eta,\bm{0},\sinh\eta)$ with $\eta$ being the pseudo-rapidity
of the fluid, and a given particle species with mass $m$ and momentum
$p^{\mu}=(m_{\perp}\cosh y,\bm{p}_{\perp},m_{\perp}\sinh y)$, where
$m_{\perp}=\sqrt{m^{2}+\bm{p}_{\perp}^{2}}$, and $y$ is the kinematic rapidity
of the particle, the momentum distribution of particles created on the freeze-
out surface is given by

$\displaystyle\frac{dN_{\rm B}}{p_{\perp}dp_{\perp}d\phi
dy}=\frac{1}{(2\pi)^{3}}\int
d^{3}\sigma_{\mu}p^{\mu}f(x,p)=\frac{m_{\perp}\tau_{f}R^{2}}{8\pi^{2}}F_{0},\quad
F_{0}=2K_{1}(\hat{m}_{\perp})+\frac{1}{12}\left[\hat{p}_{\perp}^{2}K_{1}(\hat{m}_{\perp})-2\hat{m}_{\perp}K_{2}(\hat{m}_{\perp})\right]\mathcal{A}$
(20)

where $\hat{p}_{\perp}\equiv p_{\perp}/T\equiv|\bm{p}_{\perp}|/T$,
$\hat{m}_{\perp}\equiv m_{\perp}/T$,
$d^{3}\sigma_{\mu}=(\cosh\eta,\bm{0},-\sinh\eta)\tau_{f}d\eta d^{2}x$ is the
area element of the freeze-out surface,
$f(x,p)=e^{u\cdot\hat{p}}(1+\epsilon_{\mu\nu}\hat{p}^{\mu}\hat{p}^{\nu})$,
with $\epsilon_{\mu\nu}=\pi_{\mu\nu}/2(\mathcal{E}+p)$, is the non-equilibrium
phase space distribution for classical particle, $R$ is the transverse radius
cutoff for the flow and $K_{n}(x)$ is the Bessel function of the second kind.
At large $\hat{p}_{\perp}$,
$F_{0}\sim\frac{\sqrt{2\pi}}{24}\mathcal{A}\,\hat{p}_{\perp}^{3/2}e^{-\hat{p}_{\perp}}$
(cf. the ideal limit where
$F_{0}\sim\sqrt{2\pi}\hat{p}_{\perp}^{-1/2}e^{-\hat{p}_{\perp}}$) and thus Eq.
(20) behaves like $\hat{p}_{\perp}^{5/2}e^{-\hat{p}_{\perp}}$. The ratio of
the dissipative and ideal parts in Eq. (20) is asymptotically
$\mathcal{A}\,\hat{p}_{\perp}^{2}/24$, which, at late time (NS limit), sets
the maximum value of $p_{\perp}\sim\sqrt{3/C_{\eta}}\Lambda~{}(\sim 1-2$ GeV)
where the hydrodynamic description breaks down. Indeed, Eq. (20) reduces to
its NS limit by substituting Eq. (4b) [24].

The correction $\Delta N$ due to the contribution of the linearized
perturbations up to quadratic order is given by

$\frac{d\Delta N}{p_{\perp}dp_{\perp}d\phi
dy}=\frac{m_{\perp}\tau_{f}R^{2}}{8\pi^{2}}\Big{[}F_{1}\langle{\delta\hat{T}}\rangle_{\perp}+F_{2}\,\hat{p}_{\perp}^{i}\langle{\delta
u_{i}}\rangle_{\perp}+F_{3}\,\hat{p}_{\perp}^{i}\hat{p}_{\perp}^{j}\langle{\delta\hat{\pi}_{ij}}\rangle_{\perp}+F_{11}\langle{\delta\hat{T}\delta\hat{T}}\rangle_{\perp}+F_{12}\,\hat{p}_{\perp}^{i}\langle{\delta
u_{i}\delta\hat{T}}\rangle_{\perp}\\\
+F_{13}\,\hat{p}_{\perp}^{i}\hat{p}_{\perp}^{j}\langle{\delta\hat{\pi}_{ij}\delta\hat{T}}\rangle_{\perp}+F_{22}\,\hat{p}_{\perp}^{i}\hat{p}_{\perp}^{j}\langle{\delta
u_{i}\delta
u_{j}}\rangle_{\perp}+F_{23}\,\hat{p}_{\perp}^{i}\hat{p}_{\perp}^{j}\hat{p}_{\perp}^{k}\langle{\delta\hat{\pi}_{ij}\delta
u_{k}}\rangle_{\perp}\Big{]},$ (21)

where $\hat{p}_{\perp}^{1}=\hat{p}_{\perp}\cos\phi$,
$\hat{p}_{\perp}^{2}=\hat{p}_{\perp}\sin\phi$, and

$\displaystyle
F_{1}=\hat{m}_{\perp}\left[K_{0}(\hat{m}_{\perp})+K_{2}(\hat{m}_{\perp})\right]-\frac{1}{6}\left[(\hat{m}_{\perp}^{3}+3\hat{p}_{\perp}^{2})K_{1}(\hat{m}_{\perp})+3\hat{m}_{\perp}(\hat{m}_{\perp}-2)K_{2}(\hat{m}_{\perp})\right]\mathcal{A},\quad
F_{2}=2F_{22}=F_{0},$ $\displaystyle
F_{3}=F_{23}=\frac{3}{4}K_{1}(\hat{m}_{\perp}),\quad
F_{13}=\frac{3}{4}\left[\hat{m}_{\perp}K_{0}(\hat{m}_{\perp})-5K_{1}(\hat{m}_{\perp})\right],\quad
F_{11}=-\hat{m}_{\perp}K_{0}(\hat{m}_{\perp})+\hat{m}_{\perp}^{2}K_{1}(\hat{m}_{\perp})$
$\displaystyle+\frac{1}{24}\left[-14\hat{m}_{\perp}\hat{p}_{\perp}^{2}K_{0}(\hat{m}_{\perp})+2(14\hat{m}_{\perp}^{2}+42\hat{p}_{\perp}^{2}+\hat{m}_{\perp}^{2}\hat{p}_{\perp}^{2})K_{1}(\hat{m}_{\perp})-(\hat{m}_{\perp}^{3}+12\hat{p}_{\perp}^{2})K_{2}(\hat{m}_{\perp})-\hat{m}_{\perp}^{3}K_{4}(\hat{m}_{\perp})\right]\mathcal{A},$
$\displaystyle
F_{12}=2\hat{m}_{\perp}K_{0}(\hat{m}_{\perp})+\frac{1}{24}\left[\hat{m}_{\perp}(16+\hat{p}_{\perp}^{2})K_{0}(\hat{m}_{\perp})+2(16-2\hat{m}_{\perp}^{2}-7\hat{p}_{\perp}^{2})K_{1}(\hat{m}_{\perp})+\hat{m}_{\perp}\hat{p}_{\perp}^{2}K_{2}(\hat{m}_{\perp})\right]\mathcal{A}.$
(22)

At large $\hat{p}_{\perp}$, Eq. (21) is asymptotically dominated by terms
$\langle{\delta\hat{T}\delta\hat{T}}\rangle_{\perp}$ and $\langle{\delta
u_{i}\delta u_{j}}\rangle_{\perp}$, with coefficients $F_{11}\sim
F_{22}\hat{p}_{\perp}^{2}\sim\frac{\sqrt{2\pi}}{24}\mathcal{A}\,{\hat{p}_{\perp}}^{7/2}e^{-\hat{p}_{\perp}}$,
so that Eq. (21) behaves as $\hat{p}_{\perp}^{9/2}e^{-\hat{p}_{\perp}}$.
Comparing this to the asymptotic behavior of Eq. (20) we infer that the
transverse perturbation is small relative to the background only as long as
$\hat{p}_{\perp}\lesssim 1$ (see also the right panel of Fig. 6). When
$\hat{p}_{\perp}$ becomes of order one, it would be natural to resum the
series in Eq. (21) using a Padé approximant. Nevertheless, we note that
although the large $p_{\perp}$ behavior is nonperturbative, its the magnitude
is also negligible compared to the lower $p_{\perp}$ behavior. The
distribution of total particle number $N=N_{\rm B}+\Delta N$ obtained by
adding Eqs. (20) and (21) is shown in the left panel of Fig. 6.

The dependence of the multiplicity on both $p_{\perp}$ and $\phi$ is shown in
Fig. 7. In particular, we compare the dependence of the multiplicity
distribution on $\phi$ for two different impact parameters $b$. The
distribution resulting from a central collision ($b=0$) does not depend on
$\phi$, as one would expect, due to the isotropy of the system. On the other
hand, the right panel clearly shows the dependence on the asymuthal angle for
a noncentral collision. This demonstrates the importance of the contribution
from perturbations in both magnitude of particle yields and angular
distribution in transverse plane.

##  V Outlook

Figure 7: This figure shows the 2D contour plot of the multiplicity
distribution $dN/p_{\perp}dp_{\perp}d\phi dy$ at $\tau=10$ fm with different
impact parameters, i.e., $b=0$ fm (left) and $b=6$ fm (right) respectively.

Despite the radical simplifications made, the original Bjorken model provided
a very useful analytic formula describing the dynamics of the energy density
and has lead to many important insights. The incorporation of viscous effects
and higher orders in the large proper time expansion retained some of the
simplicity of the original Bjorken model, allowing for approximate analytic
calculations and significantly simplifying numerical computations. However,
the stringent symmetry requirements — boost invariance and homogeneity in the
transverse plane – placed severe limits on how much of the physics could be
described by such models. In this article we have aimed to balance simplicity
with the ability to capture more of the interesting physics of QGP. To this
end, we have partially relaxed the symmetry requirements, insisting only on
longitudinal boost invariance while incorporating transverse dynamics by
linearizing around a boost-invariant attractor background, characterised by a
choice of the initial temperature. In this way we have formulated a
description which goes beyond toy models, but still allows some analytic
insights which can be extracted by applying known asymptotic techniques.

The initial state is encoded in the initial conditions for the Fourier modes
of the transverse perturbations and the initial temperature of the attractor
background. The evolution equations map this initial-state information into a
set of six exponentially suppressed scale-dependent amplitudes (transseries
coefficients) and a scale $\Lambda$. From the perspective of modern asymptotic
analysis our work shows that the dynamics of QGP created in heavy-ion
collision experiments provides a physical situation where the exponentially
suppressed corrections to asymptotic power series solutions carry almost all
the information which is actually detected in experiments. We have shown how
these exponential corrections are translated into the physics post-freeze-out.

Aside from providing a clear, semi-analytic picture of QGP evolution, our
approach can be implemented numerically in a very straightforward and
efficient fashion, since it relies on the discrete Fourier transforms and
solving systems of coupled linear ODEs. The details of the initial state can
be incorporated scale by scale in a controlled way: to account for finer
detail of the initial state we can increase the number of modes in the
calculation whose computational complexity then scales linearly. Apart from
the conceptual utility, our approach can thus be used as a laboratory for
studying models of the initial state, as well as novel characterisations of
fluid behaviour [25, 26]. Our model can also be extended in a number of ways;
two most prominent are the description of jets [27, 28, 29, 30], and the
incorporation of noise [31, 32, 33, 34, 35, 36]. We hope to return to these
matters in the future.

We have provided some examples which demonstrate that our approach is able to
reproduce basic features which rely on an effective description of the
transverse dynamics. It is however not yet clear which physical effects are
captured by this model, and which require a fully-nonlinear treatment. In our
pilot study we have considered the buildup of elliptic flow and have found
results consistent with earlier work which relied on solving the fully-
nonlinear problem. It would be very interesting to better understand the
limitations of our approach.

##  VI Acknowledgements

It is a pleasure to thank V. Ambrus, L. Du, W. Florkowski and S. Mrówczyński,
M. Stephanov for helpful discussions. This work was supported by the National
Science Centre, Poland, under Grant No. 2021/41/B/ST2/02909. The authors would
also like to thank the Isaac Newton Institute for Mathematical Sciences,
Cambridge, for support and hospitality during the program “Applicable
resurgent asymptotics: towards a universal theory” supported by EPSRC Grant
No. EP/R014604/1.

##  VII Appendix

##  Appendix A Linearized ODEs

The six first-order linear differential equations in (14) can be written as
three second-order linear differential equations by eliminating the variables
$\pi_{ij}$, the resulting equations take the form

$\displaystyle\partial^{2}_{\tau}\delta\hat{T}+P^{(1)}_{TT}\partial_{\tau}\delta\hat{T}(\tau)+P^{(1)}_{T\theta}\partial_{\tau}\delta\hat{\theta}+P^{(0)}_{TT}\delta\hat{T}+P^{(0)}_{T\theta}\delta\hat{\theta}=0,$
(23a) $\displaystyle\partial^{2}_{\tau}\delta\hat{\theta}+P^{(1)}_{\theta
T}\partial_{\tau}\delta\hat{T}+P^{(1)}_{\theta\theta}\partial_{\tau}\delta\hat{\theta}+P^{(0)}_{\theta
T}\delta\hat{T}+P^{(0)}_{\theta\theta}\delta\hat{\theta}=0,$ (23b)
$\displaystyle\partial^{2}_{\tau}\delta\hat{\omega}+P^{(1)}_{\omega\omega}\partial_{\tau}\delta\hat{\omega}+P^{(0)}_{\omega\omega}\delta\hat{\omega}=0,$
(23c)

with the matrix coefficients $P_{AB}^{(n)}$ where $A,B=T,\theta,\omega$ and
$(n)$ denotes the differentiation order of the corresponding term:

$\displaystyle P^{(1)}_{TT}=\frac{9(1+w)+4\mathcal{A}}{9\tau},\quad
P^{(1)}_{T\theta}=\frac{k(12+\mathcal{A})}{36},\quad
P^{(0)}_{TT}=\frac{-288\alpha^{2}+\left(45w+8\mathcal{A}+36\tau\partial_{\tau}\right)\mathcal{A}}{162\tau^{2}},$
$\displaystyle
P^{(0)}_{T\theta}=\frac{k}{324\tau}\left[36(3+2\alpha^{2}+3w)+(57+9w+2\mathcal{A}+9\tau\partial_{\tau})\mathcal{A}\right],\quad
P^{(1)}_{\theta T}=\frac{12k}{12+\mathcal{A}},$ $\displaystyle
P^{(1)}_{\theta\theta}=P^{(1)}_{\omega\omega}=\frac{-36(1-2\alpha^{2}-3w)+(45+2\mathcal{A}+9\tau\partial_{\tau})\mathcal{A}}{9(12+\mathcal{A})\tau},\quad
P^{(0)}_{\theta
T}=-\frac{k[12(8\alpha^{2}+3w)+\left(8-3w\right)\mathcal{A}]}{3(12+\mathcal{A})\tau},$
$\displaystyle
P^{(0)}_{\theta\theta}=\frac{1}{27(12+\mathcal{A})\tau^{2}}\Big{\\{}108\left[(1-2\alpha^{2})(1-w)+4\alpha^{2}k^{2}\tau^{2}\right]$
$\displaystyle+\left[-3(29-16\alpha^{2}-12k^{2}\tau^{2}+9w\tau\partial_{\tau}\ln
T)+(7-3w)(9w+2\mathcal{A}+9\tau\partial_{\tau})\right]\mathcal{A}\Big{\\}},$
$\displaystyle
P^{(0)}_{\omega\omega}=\frac{1}{27(12+\mathcal{A})\tau^{2}}\Big{\\{}108\left[(1-2\alpha^{2})(1-w)+3\alpha^{2}k^{2}\tau^{2}\right]$
$\displaystyle+\left[-3(29-16\alpha^{2}+9w\tau\partial_{\tau}\ln
T)+(7-3w)\left(9w+2\mathcal{A}+9\tau\partial_{\tau}\right)\right]\mathcal{A}\Big{\\}},$
(24)

where $w\equiv\tau T/C_{\tau}$. One immediately notices that the equation for
vorticity $\delta\omega$ in Eq. (23) decouples from the other two equations.
Given a solution to these equations, one can recover the independent shear-
stress tensor components algebraically.

Substituting Eq. (7) into Eqs. (23), and taking $C_{\tau}=0$ and
$C_{\eta}/\tau T\ll 1$, Eqs. (23) reduces to its the NS limit as those
presented also in Ref. [37]:

$\displaystyle\partial_{\tau}\delta\hat{T}+\frac{k}{3}\left(1+\frac{4C_{\eta}}{3\tau
T}\right)\delta\hat{\theta}=0,$ (25a)
$\displaystyle\partial_{\tau}\delta\hat{\theta}-\frac{1}{3\tau}\left(1-\frac{4C_{\eta}k^{2}\tau}{T}\right)\delta\hat{\theta}-k\left(1+\frac{4C_{\eta}}{3\tau
T}\right)\delta\hat{T}=0,$ (25b)
$\displaystyle\partial_{\tau}\delta\hat{\omega}-\frac{1}{3\tau}\left(1-\frac{3C_{\eta}k^{2}\tau}{T}\right)\delta\hat{\omega}=0.$
(25c)

The equation for the ideal fluid is simply obtained by setting $C_{\eta}=0$ in
Eqs. (25).

##  Appendix B The late-time asymptotic solutions in Navier-Stokes and ideal
limits

The asymptotic solutions shall be obtained using Eqs. (14). In the NS limit,
$C_{\tau}=0$, we find

$\displaystyle\delta\hat{T}_{\rm NS}$
$\displaystyle=D_{1}(\Lambda\tau)^{-\beta_{4}^{\rm
NS}}\exp\left[-\frac{C_{\eta}k^{2}}{\Lambda^{2}}(\Lambda\tau)^{4/3}+\frac{3}{8C_{\eta}}(\Lambda\tau)^{2/3}\right]+D_{2}(\Lambda\tau)^{\beta_{4}^{\rm
NS}}\exp\left[-\frac{3}{8C_{\eta}}(\Lambda\tau)^{2/3}\right],$ (26a)
$\displaystyle\delta\hat{\theta}_{\rm NS}$
$\displaystyle=D_{1}^{\prime}(\Lambda\tau)^{-\beta_{4}^{\rm
NS}+\frac{1}{3}}\exp\left[-\frac{C_{\eta}k^{2}}{\Lambda^{2}}(\Lambda\tau)^{4/3}+\frac{3}{8C_{\eta}}(\Lambda\tau)^{2/3}\right]+D_{2}^{\prime}(\Lambda\tau)^{\beta_{4}^{\rm
NS}-\frac{1}{3}}\exp\left[-\frac{3}{8C_{\eta}}(\Lambda\tau)^{2/3}\right],$
(26b) $\displaystyle\delta\hat{\omega}_{\rm NS}$
$\displaystyle=D_{3}(\Lambda\tau)^{1/3}\exp\left[-\frac{3C_{\eta}k^{2}}{4\Lambda^{2}}(\Lambda\tau)^{4/3}+\frac{8C_{\eta}}{3}(\Lambda\tau)^{-2/3}\right],$
(26c)

where $\beta_{4}^{\rm
NS}=-\frac{2}{3}-\frac{3\Lambda^{2}}{64C_{\eta}^{3}k^{2}}$ and $D$’s are
integration constants.

In the ideal limit, $C_{\tau}=C_{\eta}=0$, from Eqs. (25) with $C_{\eta}=0$
one immediately obtains, with $E$’s being integration constants, that

$\displaystyle\delta\hat{T}_{\rm
I}=(\Lambda\tau)^{1/6}\left(E_{1}e^{-\frac{i}{\sqrt{3}}k\tau}+E_{2}e^{\frac{i}{\sqrt{3}}k\tau}\right),\quad\delta\hat{\theta}_{\rm
I}=i\sqrt{3}(\Lambda\tau)^{1/6}\left(E_{1}e^{-\frac{i}{\sqrt{3}}k\tau}-E_{2}e^{\frac{i}{\sqrt{3}}k\tau}\right),\quad\delta\hat{\omega}_{\rm
I}=E_{3}(\Lambda\tau)^{1/3}.$ (27)

## References

  * [1] J.D. Bjorken, _Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region_ , _Phys. Rev. D_ 27 (1983) 140.
  * [2] M.P. Heller and M. Spaliński, _Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation_ , _Phys. Rev. Lett._ 115 (2015) 072501 [1503.07514].
  * [3] P. Romatschke, _Relativistic Fluid Dynamics Far From Local Equilibrium_ , _Phys. Rev. Lett._ 120 (2018) 012301 [1704.08699].
  * [4] A. Kurkela, W. van der Schee, U.A. Wiedemann and B. Wu, _Early- and late-time behavior of attractors in heavy-ion collisions_ , _Physical Review Letters_ 124 (2020) .
  * [5] J.-P. Blaizot and L. Yan, _Fluid dynamics of out of equilibrium boost invariant plasmas_ , _Phys. Lett. B_ 780 (2018) 283 [1712.03856].
  * [6] I. Muller, _Zum Paradoxon der Warmeleitungstheorie_ , _Z. Phys._ 198 (1967) 329.
  * [7] W. Israel, _Nonstationary irreversible thermodynamics: A Causal relativistic theory_ , _Annals Phys._ 100 (1976) 310.
  * [8] R. Baier, P. Romatschke, D.T. Son, A.O. Starinets and M.A. Stephanov, _Relativistic viscous hydrodynamics, conformal invariance, and holography_ , _JHEP_ 04 (2008) 100 [0712.2451].
  * [9] M. Spaliński, _Small systems and regulator dependence in relativistic hydrodynamics_ , _Phys. Rev. D_ 94 (2016) 085002 [1607.06381].
  * [10] W. Florkowski, M.P. Heller and M. Spaliński, _New theories of relativistic hydrodynamics in the LHC era_ , _Rept. Prog. Phys._ 81 (2018) 046001 [1707.02282].
  * [11] J. Jankowski and M. Spaliński, _Hydrodynamic Attractors in Ultrarelativistic Nuclear Collisions_ , 2303.09414.
  * [12] A. Soloviev, _Hydrodynamic attractors in heavy ion collisions: a review_ , _Eur. Phys. J. C_ 82 (2022) 319 [2109.15081].
  * [13] M. Spaliński, _Initial State and Approach to Equilibrium_ , _Acta Phys. Polon. Supp._ 16 (2023) 9 [2209.13849].
  * [14] I. Aniceto, D. Hasenbichler and A.O. Daalhuis, _The late to early time behaviour of an expanding plasma: hydrodynamisation from exponential asymptotics_ , _J. Phys. A_ 56 (2023) 195201 [2207.02868].
  * [15] M.P. Heller, R. Jefferson, M. Spaliński and V. Svensson, _Hydrodynamic Attractors in Phase Space_ , _Phys. Rev. Lett._ 125 (2020) 132301 [2003.07368].
  * [16] I. Aniceto and M. Spaliński, _Resurgence in Extended Hydrodynamics_ , _Phys. Rev. D_ 93 (2016) 085008 [1511.06358].
  * [17] C.M. Bender and S.A. Orszag, _Advanced Mathematical Methods for Scientists and Engineers_ , McGraw-Hill (1978).
  * [18] W. Wasow, _Asymptotic expansions for ordinary differential equations_ , Pure and Applied Mathematics, Vol. XIV, Interscience Publishers John Wiley & Sons, Inc., New York-London-Sydney (1965).
  * [19] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, _Numerical Recipes in C_ , Cambridge University Press, Cambridge, USA, second ed. (1992).
  * [20] P.F. Kolb and U. Heinz, _Hydrodynamic description of ultrarelativistic heavy-ion collisions_ , 2003.
  * [21] U. Heinz, _Early collective expansion: Relativistic hydrodynamics and the transport properties of QCD matter_ , in _Relativistic Heavy Ion Physics_ , pp. 240–292, Springer Berlin Heidelberg (2010), DOI.
  * [22] D.A. Teaney, _Viscous Hydrodynamics and the Quark Gluon Plasma_ , in _Quark-gluon plasma 4_ , R.C. Hwa and X.-N. Wang, eds., pp. 207–266 (2010), DOI [0905.2433].
  * [23] F. Cooper and G. Frye, _Single-particle distribution in the hydrodynamic and statistical thermodynamic models of multiparticle production_ , _Phys. Rev. D_ 10 (1974) 186.
  * [24] D. Teaney, _The Effects of viscosity on spectra, elliptic flow, and HBT radii_ , _Phys. Rev. C_ 68 (2003) 034913 [nucl-th/0301099].
  * [25] J.-Y. Ollitrault, _Measures of azimuthal anisotropy in high-energy collisions_ , _Eur. Phys. J. A_ 59 (2023) 236 [2308.11674].
  * [26] V.E. Ambrus, S. Schlichting and C. Werthmann, _Establishing the Range of Applicability of Hydrodynamics in High-Energy Collisions_ , _Phys. Rev. Lett._ 130 (2023) 152301 [2211.14356].
  * [27] A.K. Chaudhuri and U. Heinz, _Effects of jet quenching on the hydrodynamical evolution of quark-gluon plasma_ , _Physical Review Letters_ 97 (2006) .
  * [28] P.M. Chesler and L.G. Yaffe, _The Wake of a quark moving through a strongly-coupled plasma_ , _Phys. Rev. Lett._ 99 (2007) 152001 [0706.0368].
  * [29] J. Casalderrey-Solana, J.G. Milhano, D. Pablos, K. Rajagopal and X. Yao, _Jet Wake from Linearized Hydrodynamics_ , _JHEP_ 05 (2021) 230 [2010.01140].
  * [30] Z. Yang, T. Luo, W. Chen, L. Pang and X.-N. Wang, _3d structure of jet-induced diffusion wake in an expanding quark-gluon plasma_ , _Physical Review Letters_ 130 (2023) .
  * [31] Y. Akamatsu, A. Mazeliauskas and D. Teaney, _Kinetic regime of hydrodynamic fluctuations and long time tails for a bjorken expansion_ , _Phys. Rev. C_ 95 (2017) 014909.
  * [32] X. An, G. Başar, M. Stephanov and H.-U. Yee, _Relativistic Hydrodynamic Fluctuations_ , _Phys. Rev. C_ 100 (2019) 024910 [1902.09517].
  * [33] X. An, G. Başar, M. Stephanov and H.-U. Yee, _Fluctuation dynamics in a relativistic fluid with a critical point_ , _Phys. Rev. C_ 102 (2020) 034901 [1912.13456].
  * [34] X. An, G. Başar, M. Stephanov and H.-U. Yee, _Evolution of non-gaussian hydrodynamic fluctuations_ , _Phys. Rev. Lett._ 127 (2021) 072301 [2009.10742].
  * [35] X. An, G. Başar, M. Stephanov and H.-U. Yee, _Non-Gaussian fluctuation dynamics in relativistic fluids_ , _Phys. Rev. C_ 108 (2023) 034910 [2212.14029].
  * [36] Z. Chen, D. Teaney and L. Yan, _Hydrodynamic attractor of noisy plasmas_ , 2206.12778.
  * [37] S. Floerchinger and U.A. Wiedemann, _Fluctuations around Bjorken Flow and the onset of turbulent phenomena_ , _JHEP_ 11 (2011) 100 [1108.5535].