Fast resonance decays in nuclear collisions

In the context of ultra-relativistic nuclear collisions, we present a fast method for calculating the final particle spectra after the direct decay of resonances from a Cooper–Frye integral over the freeze-out surface. The method is based on identifying components of the final particle spectrum that transform in an irreducible way under rotations in the fluid-restframe. Corresponding distribution functions can be pre-computed including all resonance decays. Just a few of easily tabulated scalar functions then determine the Lorentz invariant decay spectrum from each space-time point, and simple integrals of these scalar functions over the freeze-out surface determine the final decay products. This by-passes numerically costly event-by-event calculations of the intermediate resonances. The method is of considerable practical use for making realistic data to model comparisons of the identified particle yields and flow harmonics, and for studying the viscous corrections to the freeze-out distribution function.


Introduction
Ultra-relativistic heavy ion collisions create the deconfined state of matter called the Quark-Gluon Plasma (QGP), which has been under intensive experimental and theoretical research in the last two decades [1,2]. Remarkably, the expanding QGP has been very successfully described as a relativistic fluid, where the system dynamics is completely determined by a few macroscopic fields like fluid velocity u μ (x) or temperature T (x) [3][4][5][6][7]. As the fluid expands and cools down below the cross-over temperature T c ≈ 155 MeV, quarks and gluons are re-confined in hadronic degrees of freedom. Therefore a systematic comparison between the hydrodynamic models of the QGP and experimental data necessitates the conversion of hydrodynamic fields into hadronic degrees of freedom. a e-mail: a.mazeliauskas@thphys.uni-heidelberg.de Various techniques of treating the hadronic phase have been developed over the years. Resonances are sampled at the freeze-out surface using the Cooper-Frye formula [8] and then passed to hadronic transport models, which describe both the decays and possible rescatterings of resonances [9,10]. However direct resonance decays (without rescatterings) are often used in phenomenological studies [11][12][13][14]. The decay processes of resonances are simulated by Monte-Carlo generators [15][16][17], or by semi-analytic treatments of decay integrals [18,19]. From ∼ 300 species of hadronic resonances produced in high energy nuclear collisions, only a handful long-lived hadrons (e.g. pion, kaons and protons) reach the particle detectors and are directly observed [20,21]. In this paper we show how to by-pass the numerically costly procedure of calculating the intermediate resonance decays and to relate directly the hydrodynamic fields on the freezeout surface to the final decay particle spectra.
Let us remark here that semi-analytic description of resonance decays was studied previously [18,19,22,23]. While these works constitute the basis of our approach, our framework is applicable to arbitrary freeze-out surfaces and more general particle distribution functions.
In Sect. 2 we describe a particle decay process as a linear Lorentz invariant map, which transforms the spectrum of initial (or primary) particles to the spectrum of final decay products. In Sect. 3 we argue that the decay map can be applied directly to the particle distribution function before performing the Cooper-Frye integral and we define a distribution function for the decay products. Using group theoretical arguments we find the Lorentz invariant decomposition of the decay particle distribution functions as a sum of frameindependent weight functions, which we calculate. In Sect. 4 we show that the same procedure also applies to viscous and linear perturbations of particle distribution function on the freeze-out hypersurface. Then in Sect. 5 we discuss the implementation of fast resonance decay procedure for general freeze-out surfaces and phenomenologically convenient setups of blast-wave freeze-out and mode-by-mode hydro-dynamics. We end with discussion of further extensions and applications in Sect. 6. Finally Appendix A gives the derivation of the irreducible decay spectrum components.

Lorentz invariant decay map
Relativistic particle decays is a well established subject [18,19,21,[24][25][26] and here we briefly discuss some of the key formulas. In kinetic theory, decays can be treated as 1 ↔ n particle scatterings. The probability for such an event is given by a Lorentz invariant integral over the scattering matrix squared |M| 2 , the available momentum phase space (constrained by 4-momentum conservation) and the phase space densities of initial and final particles, i.e. the gain and loss terms. In chemical and thermal equilibrium both the decay and the reverse process are equally likely, however, if the system becomes dilute and falls out of the detailed balance, multi-particle scatterings become very improbable and the 1 → n decays, which are linear in the initial particle densities, dominates the subsequent phase space evolution of particles. This is exactly what happens for hadron resonances in heavy ion collisions below the freeze-out temperature. At sufficiently late times, all allowed decays will have taken place and the Lorentz invariant spectrum of the final particle species b will be proportional to the primary populations of resonance species a, which decay (directly or through intermediate resonances) to particles of type b. 1 A particle decay is intrinsically a probabilistic process, and the resultant particle spectrum from a decay cascade will fluctuate event by event, but for very large number of initial resonances (or an average over many decay cascades), we can write the 1-body particle spectra of final particles as a Lorentz invariant integral over the primary resonances. 2 The linear decay map D a b (p, q) simply gives the Lorentz invariant probability of particle a with momentum q to decay to a particle b with corresponding momentum p. Summing over all species of primary resonances a then gives the total decayed particle spectrum of particle species b. We note in passing that the decay map D a b (p, q) fulfils certain sum rules as a consequence of conservation laws for energy, momentum, net baryon number, electric charge, etc. 1 This applies to all strong decays of hadron resonances, which have typical lifetimes τ ∼ 10 −23 s and decay before reaching the detector. Some weak decays of strange particles take place within the detector and can be reconstructed experimentally [27,28]. 2 We also note here that the two-particle distribution E p E q d N bc /(d 3 pd 3 q) inherits correlations from a decay a → b + c as a consequence of energy and momentum conservation. However, in the present paper we only concentrate on one-particle distributions.
In general the linear map D a b (p, q) is a composition of phase space integrals, 4-momentum conservation and decay matrix elements for each successive decay in the decay cascade [19,29]. Most of the listed decays or hadron resonances in the Particle Data Group book [21] are 2-body and 3body decays, which in heavy ion simulations are customary approximated as isotropic decays with a branching ratio B. For the simple case of isotropic two-body decay a → b + c the phase space integral of the decay partner c can be done analytically and the map D a b|c is reduced to a Lorentz invariant delta function of the product of initial and final particle momenta p μ q μ [19,21,24] where B is the branching ratio for this process. In the restframe of particle a, Eq. (2) is simply a uniform probability distribution on a sphere with radius |p| = p a b|c fixed by energy conservation and we also use E a b|c ≡ m 2 b + ( p a b|c ) 2 . Isotropic three body decays a → b + c + d have larger phase-space and 4-momentum conservation is not enough to fix the particles' momenta even in the rest-frame of the primary resonance a. However treating the two partner particles c and d as a fictitious particlec with an effective mass m 2 c = −( p c + p d ) 2 , the three body decay map D a b|cd can be written as an integral of the 2-body decay map for the allowed values of mc [19,21,24] where p a b|c is the momentum of b in the rest-frame of a and pc c|d is the momentum of particle c or d in their common rest-frame [21].
Composing 2-body and 3-body decay maps, Eq. (2) and Eq. (4), according to the chain rule a map for an arbitrary decay cascade can be constructed. Importantly, this a → b + X decay map will be itself only a function of the invariant product p μ q μ of initial and final particle momenta. Such a map can be calculated iteratively by applying the 2-body decay map Eq. (2) for which the chain rule Eq. (5) can be reduced to just a single dimensional integral over dimensionless variable w.
The extension to three body decays follows immediately by the application of Eq. (4). The decay chain map D a b (q μ p μ ) is independent of initial particle spectrum and only depends on particle properties and branching ratios listed in the particle data book [21]. This means that the main computational cost is in the evaluation of the decay map, which only needs to be done once, and then the final decay spectrum can be computed from an arbitrary initial particle spectrum according to Eq. (1). In particular, more finer details of resonance decay processes could be thus efficiently treated. The primary example is a finite width of resonances, which can be included in the decay map as an additional integral over resonance mass with, for example, Breit-Wigner distribution [15]. The formalism may also be generalized to anisotropic as well as spin-dependent decays. However, even ignoring these additional complications, the sheer number of primary resonances in heavy ion collisions makes the numerical evaluation of Eq. (1) a burden. Therefore we will now specialize to the decays of initial resonance spectrum specified by a common freeze-out procedure and will leave the inclusions of resonance widths and other improvements of the decay map for future work.

Cooper-Frye for the final decay spectrum
Even after integration of the intermediate resonances in the decay map, evaluation of the total decayed particle spectrum requires the sum over a (possibly large) number of primary resonances and corresponding freeze-out integrals. This sum can be performed explicitly, if we make some assumptions about the initial resonance spectra. In the hydrodynamic description of heavy ion collisions, the initial hadron spectrum is given as an integral over a freeze-out hypersurface σ according to the Cooper-Frye formula [8] 4 3 In the rest-frame of the particle b the variable w has the physical interpretation as the fraction of particle's a momentum in the direction of the fluid velocity. The limit of the massless final state particle m b = 0 can be treated by a simple change of variables u = (1 − w)  4 The hypersurface element is dσ μ = d 3 x √ hn μ , where h is the determinant of the induced metric on the freeze-out surface, d 3 x √ h is the invariant volume element and n μ is a normal vector on the surface, which we take to be pointing inwards. In this work we use mostly positive metric convention.
where ν a is the degeneracy factor of spin/polarization states and f a is a particle distribution function which depends on local fluid temperature T (x), flow velocity u μ (x), and chemical potential μ(x). We will discuss more general initial particle distributions arising in dissipative hydrodynamics in Sect. 4. In the calculation of the decay particle spectrum the order of surface integral and the linear map given by Eq. (1) can be reversed, resulting in the formula for the final decay particle spectrum where we define vector distribution function g μ , which for the primary resonances is g μ a = f a p μ , while for the decay products it is given by Once the function g μ b ( p, u, T, μ) is calculated and stored, the final decay spectra can be straightforwardly obtained by Cooper-Frye integral, Eq. (8), without the need of ever calculating distribution of intermediate hadrons.
If the initial distribution f a is only a function of particle energyĒ q = −q μ u μ in the reference frame moving with velocity u μ5 and some Lorentz scalars, e.g. temperature T or chemical potentials μ, then by Lorentz invariance of the decay process, the vector distribution function before and after the decay integral in Eq. (9) can be uniquely written as a sum of two scalar functions Here p μ , andĒ p u μ are the only available Lorentz vectors, and f eq 1 and f eq 2 are functions only depending on the Lorentz scalarĒ p , and (implicitly) μ, T , and decay parameters. In the fluid-restframe, 4-vectorsĒ p u μ = (Ē p , 0) and p μ −Ē p u μ = (0,p) are two irreducible SO(3) representations transforming under rotations as a scalar and a vector, respectively. The decay operator in Eq. (9) is a linear map and therefore guarantees that f eq 1 and f eq 2 components do not mix during (isotropic) decays. The initial hadrons on the freeze-out surface are initialized by g μ a = f a p μ and for the equilibrium distribution function both components f eq 1 and f eq 2 are initialized to be either Bose-Einstein or Fermi-Dirac distributions where μ = Q μ Q Q represents the sum over the product of all relevant chemical potentials and corresponding charges. Instead of applying the full decay map D a b (q μ p μ ) in Eq. (9) and calculating immediately the final vector distribution function g μ b (from which its components f eq i,b can be determined), one can also apply repeatedly the elementary 2-body and 3-body decay maps. This procedure is simple, because for the isotropic 2-body decay a → b + c in Eq. (2), the transformation rule between the parent and child components f a i and f b i is simply a one dimensional integral. Leaving the details for Appendix A, the iterative relation between the components is (c.f. Eq. (6)) where we used the abbreviations andĒ p andp are the energy and momentum of particle b in the fluid-restframe. The isotropic three body decays a → b+c+d can be easily incorporated by integrating the 2-body transformation rules Eq. (12) over the effective decay partner mass mc as in Eq. (4). Such one-dimensional integrals can be easily done by standard numerical integration routines [30]. For concreteness consider the decay of h 1 mesons illustrated in Fig. 1. The initial irreducible components f eq i,h 1 for h 1 meson are initialized by the Bose-Einstein distribution depending on temperature and chemical potential Then the two body decay h 1 → ρπ produces contributions to ρ and π mesons, f h 1 i,ρ and f h 1 i,π , according to Eq. (12). These have to be added to the corresponding thermal distributions and the feed-down from other resonances. Here f ρ i,π represents the total feed-down from ρ → ππ decay irrespective of ρ's origin and which is calculated according to Eq. (12) from parent particle distribution f eq i,ρ . By starting from the heaviest resonance and summing the thermal and decay contributions of lower mass resonances, the irreducible components of the final stable particles can be calculated with the minimal number of decay integrals.
The physical meaning of the irreducible components f of the vector distribution function g μ b can be clarified by considering particle spectrum in the fluid-restframe where now f eq 1,b is the part of particle spectrum proportional to p i dσ i element of the freeze-out surface, while f eq 2,b contributes to the spectrum with E p dσ 0 weight. In Fig. 2 we show the irreducible components f eq 1 and f eq 2 for the final pion π + spectrum from a completely decayed thermal T fo = 145 MeV distributions of hadron resonances as used in Monte-Carlo decay chain generators [17]. 6 We see that the f eq 2 component, which gives the sole contribution to the particle spectra for time-like (fixed time) freeze-out surface, is larger than thermal pion distribution due to feed-down from resonance decays, while the space-like component f eq 2 remains of the same size. In an arbitrary reference frame the decay pion spectrum can be straightforwardly calculated using frame independent formulas Eq. (10) and Eq. (8). We will discuss this further in Sect. 5.  The decay pion π + spectra for a simple freeze-out surface 7 for different decay channels calculated using irreducible weight functions (see Fig. 2). Results from a Monte-Carlo generator are shown for comparison [17] In Fig. 3 we plot the final pion spectra π + for a simple freeze-out surface with a constant Bjorken time, freeze-out temperature and radial velocity. 7 In addition to the total pion spectrum (which includes all decay chains producing π + ), we also show the pion spectrum from the dominant decay channels ρ +,− → π + + π 0,− and ω 0 → π + + π − + π 0 (where ρ +,0 and ω 0 spectra themselves include decay contributions from yet heavier resonances). We compare our results with the decay pion spectrum generated by a Monte-Carlo resonance decay generator THERMINATOR 2 [17]. All spectra are in excellent agreement, however we would like to stress that the decay pion spectrum in Fig. 3 is obtained immediately from a simple Cooper-Frye freeze-out proce- 7 We used the following THERMINATOR 2 options for the freeze-out surface: τ fo = 8.17 fm, T fo = 145 MeV, radius of the surface R = 8.21 fm and a constant radial velocity v T = 0.341. dure Eq. (8). The vector distribution function components f eq i shown in Fig. 2 only need to calculated once for a particular freeze-out temperature T fo and then can applied to any shape of the freeze-out surface or the fluid velocity field u μ (x), without the need of costly calculations of intermediate resonances.
Although we gave an example of a constant temperature and chemical potential freeze-out surface, our method can be equally well applied for varying freeze-out temperature or chemical potential. In such case irreducible components of the vector distribution function f eq i will need to be tabulated not only in the fluid-restframe energyĒ p = −u μ p μ , but also the additional freeze-out variables. However, since this tabulation needs only to be done once, the freeze-out integral Eq. (8) can be performed essentially without additional computational cost.

Viscous and linear corrections to particle spectrum
In viscous hydrodynamics, the freeze-out distribution function differs from the equilibrium Bose-Einstein or Fermi-Dirac distribution with additional dependences on the dissipative terms like the shear-stress tensor π μν (x) and the bulkviscous pressure (x), so that the initial vector distribution function in the Cooper-Frye formula, Eq. (7), is The functional dependence of the distribution function on viscous corrections for hadron resonance gas is largely unresolved problem even at linear order in the dissipative terms and various parametrizations are used [31,32]. Therefore below we also consider only linear dissipative corrections to the decay particle spectrum, but note that higher order terms, if known, could be also straightforwardly included. For concreteness we consider the form of viscous δ f corrections used in modern hydrodynamic simulations [32] δ Here c s (T ) is the speed of sound of the medium at the freezeout temperature, m -mass of the primary resonance, and τ /ζ is the ratio of bulk relaxation time and bulk viscosity. We want to compute the final decay particle spectrum arising from such viscous components. The procedure is analogous to the one described in the previous section. First we expand Eq. (18) to linear order in viscous corrections  (23) and (23), which determine the final pion spectrum from bulk and shear viscous perturbations of the distribution functions of primary resonances, Eqs. (19) and (20). The explicit prefactors indicate the required freeze-out surface elements for the Cooper-Frye integration in the fluid-restframe. Note that a term T 2 /|p| 2 was factored out from f shear i components. Bulk/shear perturbation of (initial) pion distribution function shown for comparison where the derivatives g μρσ shear ≡ ∂g μ /∂π ρσ and g μ bulk ≡ ∂g μ /∂ can only be functions of 4-vectors p μ andĒ p u μ , and Lorentz scalars like temperature, chemical potential or resonance mass. Initially they are given by Eqs. (19) and (20) After the decays g μ bulk and g μρσ shear can be written uniquely as a certain sum of Lorentz vectors/tensors. In Appendix A we discuss the general irreducible decomposition of Lorentz tensors in terms of representations transforming differently under SO(3) rotations in the fluid-restframe. Here we only reproduce the final result for the bulk and shear perturbations The bulk pressure perturbation does not introduce new tensor structures and the decomposition is the same as for the equilibrium distribution in Eq. corresponding to a symmetric traceless tensor, vector and scalar representations (see Appendix A). The irreducible weight functions f i of final decay particle distribution can be calculated iteratively using similar integrals as for the equilibrium distribution Eq. (12) where integration measures A i (w) are listed in Appendix A.
In Fig. 4 we show the final decayed pion π + spectrum in the fluid-restframe due to viscous perturbations at T fo = 145 MeV. 8 Different lines in Fig. 4 correspond to different contributions stemming from components f i in Eq. (23) and Eq. (24). The labels next to the lines indicate the required factors for the Cooper-Frye freeze-out integral in the fluidrestframe. Note that we factored in |p| 2 for the shear perturbations. We also factored out the terms proportional to the transport coefficients, so any (small) viscous perturbation will produce the same correction to the particle spectrum (up to the magnitude) in the local fluid-restframe. However the presence of such viscous corrections in the hydrodynamic evolution modifies the fluid velocity and temperature fields, and therefore the freeze-out surface will itself be different. Then evaluating the generalized Cooper-Frye freeze-out integral Eq. (8) will yield different particle spectrum. Similarly to viscous perturbations, one can also consider linear perturbations of fluid velocity δu μ , or temperature δT around the background fieldsū μ andT , 9 The irreducible decomposition for temperature perturbation is the same as for the equilibrium distribution The irreducible decomposition for velocity perturbations consists of symmetric traceless tensor, vector and scalar representations with corresponding weight functions f The total vector distribution function is then given by 9 Note that a constant temperature freeze-out surface depends on δT . However, one can also consider a freeze-out at a constant background temperatureT , which is then independent of the perturbation.
In Fig. 5 we show the particle spectrum decomposition due to temperature and velocity perturbation for the same freeze-out temperature. As before we factor out the explicit dependence on the magnitude of the perturbation in the fluid-restframe. Such linearised perturbations can be used to study, for example, the angular velocity modulations around a known freeze-out flow u μ , e.g. provided by a blast-wave model or in the mode by mode description of heavy ion collisions [34,35]. We would like to note in passing that other types of perturbations to the equilibrium spectrum can be considered. At lower collision energies, as in the Beam Energy Scan studies, the freeze-out distribution function will also depend on the chemical potential μ Q (x) and one could consider linear perturbations of the chemical potential δμ Q , which are treated identically to the temperature variations Eq. (28). Also, the hydrodynamics at non-zero baryon density must also evolve the baryon current j μ = n B u μ + j μ D (e.g. see [36]). The transverse (diffusion) part of this current will induce perturbations of the freeze-out surface distribution functions analogous to velocity perturbations in Eq. (29) where n B is the local baryon density, Q B is the baryon charge andκ B is a related transport coefficient [37]. Such additional terms are easy to accommodate in our framework, because the same transformation rules can be used to find the final decay spectrum (see Appendix A).

Application to blast-wave and mode-by-mode freeze-out
In this section we will discuss practical implementation of the fast resonance decay procedure in heavy ion collisions. Collecting together all equilibrium and viscous terms contributing to the particle spectrum, we can group them by how they are contracted with the freeze-out surface element dσ μ where explicitly these terms are The required values of fluid velocity u μ , shear and bulk stresses π μν , on the freeze-out surface must be provided by the hydrodynamic model of the QGP fireball or freezeout parametrization, e.g. blast-wave model [38]. Then the complete final decay particle spectrum (for a particle species or even for the total sum of charged particles) can be computed according to Eq. (32) by using the tabulated values of irreducible weight functions f i . Up to now the freeze-out surface was left completely general. An interesting application of our formalism arises from the mode-by-mode solution to the fluid dynamic expansion of a fireball [34,35]. In that formalism, one decomposes the fluid fields (temperature, chemical potentials, fluid velocity, shear stress and bulk viscous pressure) into a background part and a fluctuating part, e.g. u μ =ū μ + δu μ . In high energy collisions, a boost invariance is a good symmetry and the collision is often parametrized in Bjorken coordinates It is particularly convenient to take the background part as symmetric with respect to azimuthal rotations and longitudinal Bjorken boosts (see e.g. [39] for the fluid equations of motion in this situation). Then the hadron spectrum resulting from the fluid fields after freeze-out can be also split into a background, which is invariant under these symmetries (now in momentum space), and a fluctuating part. The freeze-out surface in this case is given by a 1D curve in τ -r plane, which can be conveniently parametrized by (τ (α), r (α)) where α ∈ (0, 1) and Similarly by symmetry the fluid velocity has only two components and can be written in terms of a radial fluid rapiditȳ χ, Note that then the particle energy in fluid-restframe is where in the coordinate system of Eq. (34) and at space time point (τ, r, φ, η) the particle momentum components are given by Here we use the transverse momentum p T (and transverse mass m T = m 2 + p 2 T ), the particle momentum angle ϕ, and the particle momentum rapidity y.
The background contribution to the shear stress tensor can be parametrized in terms of two independent components, here taken to beπ φ φ andπ η η π τ τ = −ū rūr π φ φ +π η η ,π τr =π r τ =ū rūτ π φ φ +π η η , and the bulk viscous pressure is simply¯ . The decay hadron spectrum for azimuthally and boost invariant freeze-out surface then reduces to a single integral were the freeze-out kernels K i ( p T ,χ) are solely functions of the transverse particle momentum and radial fluid velocitȳ u r = sinhχ , and the terms proportional to viscous tensors are factored out. Analogously to the original irreducible weights f i , the freeze-out kernels can be precomputed and applied to an arbitrary freeze-out surface (τ (α), r (α)) and radial fluid rapidity profileχ(α). For the equilibrium components f eq i these kernels are given by the following rapidity and angle integrals Recall thatĒ p also depends only on the differences of η − y and φ−ϕ, Eq. (37), so the dependence on momentum rapidity y and angle ϕ disappears after the integration. The integral expressions for viscous kernels are given in Appendix B. For the simple constant time freeze-out surface used in Fig. 3 only the temporal part of the freeze-out surface contributes and the decay pion spectrum is proportional to K eq 1 ( p T ,χ = arctan v T ).
All deviations from the azimuthally and boost invariant background in mode-by-mode hydrodynamics is carried by the fluctuating part of the fluid fields and to first approximation only the linear part in these perturbations contribute to the final particle spectra. For perturbations in e.g. fluid velocity δu μ the corresponding perturbations of the distribution functions are given explicitly in Eq. (27). Because the decay operator Eq. (1) is linear, if δu μ is written as linear superposition of Fourier modes in azimuthal angle φ and space-time rapidity η one finds that the distribution of hadrons after kinetic freezeout and resonance decays depends on momentum space azimuthal angle ϕ and momentum space rapidity y via the combination e imϕ+iky with the same azimuthal wavenumber m and rapidity wavenumber k. This is a direct consequence of U(1) × R 1 symmetry, which prevents different representations labelled by m and k from mixing under linear operations. This way arbitrary linear perturbations in fluid fields can be mapped to the modes of the final particle spectra, which can be straightforwardly incorporated in the formalism of mode-by-mode hydrodynamics [34,35].

Conclusions
We presented a method to calculate the final decay spectrum of direct resonance decays directly from hydrodynamic fields on a freeze-out surface. By applying the decay map, Eq. (1), to the distribution function of primary particles before the Cooper-Frye integration, we found the (vector) distribution function of decay products, Eq. (9). By decomposing this distribution into components transforming differently under SO(3) rotations in the fluid-restframe, we expressed the final decay particle spectrum as a sum of a few Lorentz invariant weight functions and known Lorentz vectors. The explicit procedure to determine the irreducible weight functions for an arbitrary decay chain of isotropic 2-body and 3-body decays was derived and a numerical implementation was made public [30]. We considered primary hadron resonances generated by the equilibrium distribution function, viscous shear and bulk perturbations, and linearised temperature and velocity perturbations. Modifications to the particle spectrum due to variations in the chemical potential and the diffusive part of particle current can be also straightforwardly included in this framework. The final 1-body particle spectrum of decay products is then calculated from a general Cooper-Frye-type freeze-out integral, Eq. (32). The most important aspect of our method is that intermediated particle decays do not need to be calculated event-by-event. The irreducible components of the decay particle distribution function Eq. (9) are computed only once, and the spectrum of a few relevant hadron species, which includes feed-down of all direct decays, can be computed for an arbitrary freeze-out surface. This significantly reduces the computational costs of direct resonance decays.
Although our method of calculating direct resonance decays is already competitive with other treatments available on the market, the computational efficiency of our approach makes it practical to include finer details of resonances decays. For example, new hadron resonance states can be easily added to improve the agreement between the lattice QCD and hadronic equation of state [11]. Finite widths of the resonances can be incorporated in the decay map [15,40,41]. This has recently been shown to reduce the discrepancy between measured and predicted proton yield in the statistical hadronization models [42].
In this work we neglected hadronic rescatterings after the chemical freeze-out, which may change the final particle spectra, but the effect is subleading in comparison with the decay feed-down [43]. Elastic scatterings in the hadronic phase can be modeled by a hydrodynamic evolution of hadron fluid in partial chemical equilibrium [44,45]. In this scenario, the particle ratios are fixed at the chemical freeze-out temperature T chem for each species i of long lived hadrons by introducing an appropriate chemical potentials μ i (T ) for T < T chem . Subsequently, the kinetic freeze-out may take place at some lower temperature T kin . Because the primary resonance spectra are still described only by temperature T kin and the chemical potentials μ i (T kin ), the direct decays can be calculated using the techniques proposed in this work.
Another interesting generalization of the framework is to keep track of the particle spin in the decays. This could be particularly useful for the studies of vorticity polarization in heavy ion collisions [46,47]. However, in this case one might need to go beyond isotropic s-wave decays and consider more general momentum dependent decay patterns, which to our knowledge were not included in phenomenological works so far. Finally, we note that the 1-particle distribution function does not have the information of the connected two-particle function, namely the non-flow correlations of particles produced by the same resonance decay. However, the decay map Eq. (1) can be generalized to two-particle spectrum.
In summary, we believe that the computationally efficient way of computing direct resonance decays, which was presented in this paper, will be of great practical utility for phenomenological studies of heavy ion collisions and make realistic particle yield calculations much more affordable.
Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data used in the figures can be generated by the publicly available code [30].] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Tensor decomposition and decay rules
Instead of calculating the decay operator to a particle spectrum, Eq. (1), one may also apply elementary decay operators directly to the distribution function and obtain thus a modified distribution function, Eq. (9), that already includes resonance decays. A technical complication is that one needs to do this based on a Lorentz-vector form of the distribution function, which at the freeze-out surface is given by p μ f ( p, u, T, μ). (43) In ideal hydrodynamics f is simply Bose-Einstein or Fermi-Dirac distribution Note that although the decay operator is Lorentz invariant, the Lorentz boost symmetry is explicitly broken by the fluid velocity, which singles out a reference frame. The residual rotational subgroup allows to decompose g μ into two 4-vectors transforming under different representations of SO(3) in the fluid-restframe where transform as a scalar (1) and a vector (3) respectively. Both representations have identical weight functions f eq 1 = f eq 2 = f eq on the freeze-out surface, but after the resonance decays are taken into account the two functions will differ. Isotropic decays do not mix different SO(3) representations and the new weight functions are found by successively applying the decay maps Eq. (2) and Eq. (4). Before deriving the specific transformation rules for the weight functions f i , lets consider a more general case for the distribution function. In viscous hydrodynamics, the freeze-out distribution function differs from the equilibrium expression Eq. (44) and also depends on shear-stress tensor and bulk viscous pressure Close to equilibrium one may Taylor expand the vector distribution function around the equilibrium distribution function where the derivatives g μρσ shear ≡ ∂g μ /∂π ρσ and g μ bulk ≡ ∂g μ /∂ can only be functions of 4-vectors p μ andĒ p u μ , and Lorentz scalars like temperature, chemical potential or resonance mass. Similarly, we can also consider perturbations of the hydrodynamic fields around an arbitrary background, e.g. u μ =ū μ + δu μ and T =T + δT . Then the vector distribution function decomposes as Thanks to the linearity of the decay map, each term in the Taylor expansion (g μ , g μρ , g μρσ ) can be also decomposed into irreducible representations of SO(3) rotational group. For example, a two-tensor distribution function g μν (u, p) contains irreducible representations of SO(3) according to the tensor decomposition (1+3)⊗(1+3) = 2×1+3×3+5. In terms of the 4-vectorsq μ andp μ they are given as (a) two scalars, However, in practice not all of these terms are needed. Thanks to the orthogonality relationū ρ δu ρ = 0, the first terms in Eq. (50) and Eq. (51) drop out, while the antisymmetric term in Eq. (51) is not present initially and by symmetry reasons is not generated by the decay operator. Only one scalar, one vector and one tensor representation contribute to the contraction g μν velocity δu ν Analogously we list the most general decomposition of the three-tensor distribution function g μνρ into the SO (3) representations (1 + 3) ⊗ (1 + 3) ⊗ (1 + 3) = 5 × 1 + 9 × 3 + 5 × 5 + 7 a) five scalars 1)q μqνqρ , 2) |p| 2qμ νρ , and 2 perm.
Note that a third possible permutations of Eq. (60) is not linearly independent of the other two. d) one symmetric and traceless three-tensor p μpνpρ − 1 5 |p| 2 p μ νρ +p ν μρ +p ρ μν . (61) Of these only three independent SO(3) representations contribute to g μνρ shear π νρ , namely, those which are symmetric, traceless and orthogonal to fluid velocity in the indecies ν and ρ: g μνρ shear π νρ = [ p ν π νρ p ρ ( p μ −Ē p u μ ) − The extension of this procedure to quadratic and higher terms in the Taylor expansions Eq. (48) and Eq. (49) is straightforward, if laborious. Now we will discuss the transformation rules for the weight functions f i due to an elementary 2-body decay, Eq. (2). The new vector distribution function is given by the convolution with the decay map Because the Lorentz vectors corresponding to different SO(3) representations in Eq. (45) are orthogonal both before and after the decay, we can use them to project out the desired component and reduce Eq. (63) to a scalar integral, which can be simplified (the same as Eq. (6)) to a single dimensional integral The calculation is straightforward and the appropriate A i (w) functions for Eq. (45) are where we remind that E(w) and Q(w) were defined as Note that some representations in Eq. (53) and Eq. (62) are just a products of lower dimensional representations and the corresponding functions A i (w) are equal to the product of their A i 's. In summary the Eq. (64) and functions A i (w) defines a simple iterative scheme for calculating weight functions f i for different irreducible components of the vector/tensor distribution function undergoing a 2-body decay Eq. (2), which can be easily extended to a 3-body decay rule by Eq. (4). By repeated application of these transformations and summing over all possible decay chains the final decay particle components f i can be determined. A concrete realization of such scheme is made publicly available [30].

B Decay kernels for azimuthally symmetric and boost-invariant freeze-out surface
For azimuthally and boost invariant freeze-out surface the general decay particle spectrum formula Eq. (32) reduces to one dimensional integral Eq. (40), where the azimuthal and rapidity integrals are factored out in the freeze-out kernels K i . The integral formulas for them are obtained by a straightforward algebra of inserting Eqs. (35), (36), (39), and (38), in Eqs. (32) and (33), and collecting terms proportional to temporal and radial components of the freeze-out surface element. Explicitly these are