1 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^\mu (x)\) or temperature T(x) [3,4,5,6,7]. As the fluid expands and cools down below the cross-over temperature \(T_c\approx 155\,\text {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.

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 \(\sim 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 freeze-out 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 frame-independent 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 hydrodynamics. We end with discussion of further extensions and applications in Sect. 6. Finally Appendix A gives the derivation of the irreducible decay spectrum components.

2 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\leftrightarrow n\) particle scatterings. The probability for such an event is given by a Lorentz invariant integral over the scattering matrix squared \(|\mathcal 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\rightarrow 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.Footnote 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.Footnote 2

$$\begin{aligned} E_\mathbf{p}\frac{d N_b}{d^3 \mathbf{p}} = \int \! \frac{d^3 \mathbf{q}}{(2\pi )^3 2 E_\mathbf{q}} D^a_b(\mathbf{p},\mathbf{q}) \; E_\mathbf{q}\frac{d N_a}{d^3 \mathbf{q}}. \end{aligned}$$
(1)

The linear decay map \(D^a_b(\mathbf{p},\mathbf{q})\) simply gives the Lorentz invariant probability of particle a with momentum \(\mathbf{q}\) to decay to a particle b with corresponding momentum \(\mathbf{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_b^a(\mathbf{p},\mathbf{q})\) fulfils certain sum rules as a consequence of conservation laws for energy, momentum, net baryon number, electric charge, etc.

In general the linear map \(D^a_b(\mathbf{p},\mathbf{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 3-body 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\rightarrow 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^\mu q_\mu \) [19, 21, 24]

$$\begin{aligned} D^a_{b|c}(p^\mu q_\mu ) =B\frac{4\pi ^2 m_a}{p^a_{b|c}} \delta ( q^\mu p_\mu + m_a E^a_{b|c}), \end{aligned}$$
(2)

where B is the branching ratio for this process. In the rest-frame of particle a, Eq. (2) is simply a uniform probability distribution on a sphere with radius \(|\mathbf{p}|=p^{a}_{b|c}\) fixed by energy conservation

$$\begin{aligned} p^a_{b|c} \equiv \frac{1}{2m_a} \sqrt{((m_a+m_b)^2-m_c^2)((m_a-m_b)^2-m_c^2}),\nonumber \\ \end{aligned}$$
(3)

and we also use \(E^a_{b|c}\equiv \sqrt{m_b^2+(p^a_{b|c})^2}\).

Isotropic three body decays \(a\rightarrow 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 particle \(\tilde{c}\) with an effective mass \(m_{\tilde{c}}^2=-(p_c+p_d)^2\), the three body decay map \(D_{b|cd}^a\) can be written as an integral of the 2-body decay map for the allowed values of \(m_{\tilde{c}}\) [19, 21, 24]

$$\begin{aligned} D^a_{b|cd}(p^\mu q_\mu )= \frac{\int _{m_c+m_d}^{m_a-m_b}\!\!dm_{\tilde{c}}\,\, p^{a}_{b|\tilde{c}} p^{\tilde{c}}_{c| d} D^a_{b|\tilde{c}}(p^\mu q_\mu )}{\int _{m_c+m_d}^{m_a-m_b}\!\!dm_{\tilde{c}}\,\, p^{a}_{b|\tilde{c}} p^{\tilde{c}}_{c|d}}. \end{aligned}$$
(4)

where \(p^{a}_{b|\tilde{c}}\) is the momentum of b in the rest-frame of a and \(p^{\tilde{c}}_{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

$$\begin{aligned}&D^a_{b}(p^\mu q_\mu )=\int \frac{d^3\mathbf{k}}{(2\pi )^32 E_\mathbf{k}}D^e_{b}(p^\mu k_\mu )D^a_{e}(k^\nu q_\nu ) \end{aligned}$$
(5)

a map for an arbitrary decay cascade can be constructed. Importantly, this \(a\rightarrow b+X\) decay map will be itself only a function of the invariant product \(p^\mu q_\mu \) 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.Footnote 3

$$\begin{aligned} D^a_{b}(p^\mu q_\mu )=&B \frac{1}{2} \frac{ m_e^2}{m_b^2} \int _{-1}^1\! dw\, D^a_e\left( q^\mu p_\mu \frac{m_e E^{e}_{b|c}}{m_b^2}\right. \nonumber \\&\left. + w\sqrt{(q^\mu p_\mu )^2-m_a^2 m_b^2} \frac{m_e p^e_{b|c}}{m_b^2} \right) \end{aligned}$$
(6)

The extension to three body decays follows immediately by the application of Eq. (4).

The decay chain map \(D^a_{b}(q^\mu p_\mu )\) 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.

3 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 \(\sigma \) according to the Cooper–Frye formula [8]Footnote 4

$$\begin{aligned} E_\mathbf{p}\frac{d N_a}{d^3 \mathbf{p}} = \frac{\nu _a}{(2\pi )^3}\int _\sigma f_a(-u^\nu p_\nu , T, \mu ) p^\mu d\sigma _\mu , \end{aligned}$$
(7)

where \(\nu _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^\mu (x)\), and chemical potential \(\mu (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

$$\begin{aligned} E_\mathbf{p}\frac{d N_b}{d^3 \mathbf{p}} = \frac{\nu _b}{(2\pi )^3}\int _\sigma g^\mu _b(p, u, T, \mu ) d\sigma _\mu , \end{aligned}$$
(8)

where we define vector distribution function \(g^\mu \), which for the primary resonances is \(g^\mu _a = f_a p^\mu \), while for the decay products it is given by

$$\begin{aligned} g^\mu _b(p, u)\equiv \sum _a\frac{\nu _a}{\nu _b}\ \int \! \frac{d^3 q}{(2\pi )^3 2 E_\mathbf{q}} D^a_b(p^\nu q_\nu )f_a(-u^\sigma q_\sigma ) q^\mu . \end{aligned}$$
(9)

Once the function \(g^\mu _b(p, u, T, \mu )\) 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 \(\bar{E}_\mathbf{q}=-q^\mu u_\mu \) in the reference frame moving with velocity \(u^\mu \)Footnote 5 and some Lorentz scalars, e.g. temperature T or chemical potentials \(\mu \), 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

$$\begin{aligned} g^\mu _b(p,u) = f^{\text {eq}}_{1,b}(\bar{E}_\mathbf{p})\left( p^\mu -\bar{E}_\mathbf{p}u^\mu \right) + f^\text {eq}_{2,b}(\bar{E}_\mathbf{p}) \bar{E}_\mathbf{p}u^\mu . \nonumber \\ \end{aligned}$$
(10)

Here \(p^\mu \), and \(\bar{E}_\mathbf{p}u^\mu \) are the only available Lorentz vectors, and \(f^\text {eq}_1\) and \(f^\text {eq}_2\) are functions only depending on the Lorentz scalar \(\bar{E}_\mathbf{p}\), and (implicitly) \(\mu \), T, and decay parameters. In the fluid-restframe, 4-vectors \(\bar{E}_\mathbf{p}u^\mu =(\bar{E}_\mathbf{p}, \mathbf {0})\) and \(p^\mu -\bar{E}_\mathbf{p}u^\mu =(0, \bar{\mathbf{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^\text {eq}_1\) and \(f^\text {eq}_2\) components do not mix during (isotropic) decays. The initial hadrons on the freeze-out surface are initialized by \(g^\mu _a = f_a p^\mu \) and for the equilibrium distribution function both components \(f_1^\text {eq}\) and \(f_2^\text {eq}\) are initialized to be either Bose-Einstein or Fermi-Dirac distributions

$$\begin{aligned} f_\text {eq}(- u^\mu p_\mu , T, \mu ) = \left( e^{-u_\mu p^\mu /T-\mu /T}\mp 1\right) ^{-1}, \end{aligned}$$
(11)

where \(\mu =\sum _Q \mu _Q Q\) represents the sum over the product of all relevant chemical potentials and corresponding charges.

Instead of applying the full decay map \(D_b^a(q^\mu p_\mu )\) in Eq. (9) and calculating immediately the final vector distribution function \(g^\mu _b\) (from which its components \(f^\text {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\rightarrow b+c\) in Eq. (2), the transformation rule between the parent and child components \(f_i^a\) and \(f_i^b\) is simply a one dimensional integral. Leaving the details for Appendix A, the iterative relation between the components is (c.f. Eq. (6))

$$\begin{aligned} f_{1,b}^\text {eq}(\bar{E}_\mathbf{p})&= B \frac{\nu _a}{\nu _b} \frac{m_a^2}{m_b^2}\frac{1}{2}\int _{-1}^1 dw f_{1,a}^\text {eq}\left( E(w)\right) \frac{Q(w)}{|\bar{\mathbf{p}}|}, \end{aligned}$$
(12a)
$$\begin{aligned} f_{2,b}^\text {eq}(\bar{E}_\mathbf{p})&= B \frac{\nu _a}{\nu _b} \frac{m_a^2}{m_b^2}\frac{1}{2}\int _{-1}^1\! dw\, f_{2,a}^\text {eq}\left( E(w)\right) \frac{E(w)}{\bar{E}_\mathbf{p}}, \end{aligned}$$
(12b)

where we used the abbreviations

$$\begin{aligned} E(w)&\equiv \frac{m_a E^a_{b|c}\bar{E}_\mathbf{p}}{m_b^2}-w \frac{m_a p^a_{b|c}|\bar{\mathbf{p}}|}{m_b^2}, \end{aligned}$$
(13a)
$$\begin{aligned} Q(w)&= \frac{m_a E^a_{b|c}|\bar{\mathbf{p}}|}{m_b^2}-w \frac{m_a p^a_{b|c}\bar{E}_\mathbf{p}}{m_b^2}, \end{aligned}$$
(13b)

and \(\bar{E}_\mathbf{p}\) and \(\bar{\mathbf{p}}\) are the energy and momentum of particle b in the fluid-restframe. The isotropic three body decays \(a\rightarrow b+c+d\) can be easily incorporated by integrating the 2-body transformation rules Eq. (12) over the effective decay partner mass \(m_{\tilde{c}}\) 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^\text {eq}_{i,h_1}\) for \(h_1\) meson are initialized by the Bose-Einstein distribution depending on temperature and chemical potential

$$\begin{aligned} f^\text {eq}_{i,h_1}&=f_{\text {eq},h_1}(\bar{E}_\mathbf{p}, T,\mu ). \end{aligned}$$
(14)

Then the two body decay \(h_1\rightarrow \rho \pi \) produces contributions to \(\rho \) and \(\pi \) mesons, \(f^{h_1}_{i,\rho }\) and \(f^{h_1}_{i,\pi }\), according to Eq. (12). These have to be added to the corresponding thermal distributions and the feed-down from other resonances.

$$\begin{aligned} f_{i,\rho }^\text {eq}&= f_{\text {eq},\rho }(\bar{E}_\mathbf{p}, T,\mu ) + f^{h_1}_{i,\rho }+f^\text {other}_{i,\rho }, \end{aligned}$$
(15)
$$\begin{aligned} f_{i,\pi }^\text {eq}&= f_{\text {eq},\pi }(\bar{E}_\mathbf{p}, T,\mu ) + f^{h_1}_{i,\pi }+f^{\rho }_{i,\pi }+f^\text {other}_{i,\pi }. \end{aligned}$$
(16)

Here \(f^{\rho }_{i,\pi }\) represents the total feed-down from \(\rho \rightarrow \pi \pi \) decay irrespective of \(\rho \)’s origin and which is calculated according to Eq. (12) from parent particle distribution \(f_{i,\rho }^\text {eq}\). 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.

Fig. 1
figure 1

Decay cascade \(h_1\rightarrow \rho \pi \rightarrow \pi \pi \pi \). Here \(h_1\) meson has only the initial thermal distribution, while \(\rho \) and \(\pi \) receive feed-down from resonances’ decays

Full size image

The physical meaning of the irreducible components \(f^\text {eq}_{i,b}\) of the vector distribution function \(g^\mu _b\) can be clarified by considering particle spectrum in the fluid-restframe

$$\begin{aligned} E_\mathbf{p}\frac{d N_b}{d^3 \mathbf{p}}&= \frac{\nu _b}{(2\pi )^3}\int _\sigma \left[ f^{\text {eq}}_{1,b}\left( p^\mu -\bar{E}_\mathbf{p}u^\mu \right) + f^\text {eq}_{2,b} \bar{E}_\mathbf{p}u^\mu \right] d\sigma _\mu \nonumber \\&=\left. \frac{\nu _b}{(2\pi )^3}\int _\sigma ( f^\text {eq}_{1,b} p^i d\sigma _i +f^\text {eq}_{2,b} E_\mathbf{p}d\sigma _0)\right| _{u^\mu =(1,\mathbf {0})}, \end{aligned}$$
(17)

where now \(f^\text {eq}_{1,b}\) is the part of particle spectrum proportional to \(p^i d\sigma _i\) element of the freeze-out surface, while \(f^\text {eq}_{2,b}\) contributes to the spectrum with \(E_\mathbf{p}d\sigma _0\) weight. In Fig. 2 we show the irreducible components \(f^\text {eq}_1\) and \(f^\text {eq}_2\) for the final pion \(\pi ^+\) spectrum from a completely decayed thermal \(T_\text {fo}=145\,\text {MeV}\) distributions of hadron resonances as used in Monte-Carlo decay chain generators [17].Footnote 6 We see that the \(f^\text {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^\text {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.

Fig. 2
figure 2

Lorentz invariant pion \(\pi ^+\) weight functions \(f_i^\text {eq}(-u^\mu p_\mu )\), Eq. (10), which determine the final pion spectrum at each space-time point on the freeze-out surface, Eq. (17). Direct decays of \(\sim 300\) hadron resonances were computed from equilibrium distributions with temperature \(T_\text {fo}=145\,\text {MeV}\) and zero chemical potential \(\mu _Q=0^{6}\). The explicit prefactors indicate the required freeze-out surface elements for the Cooper–Frye integration in the fluid-restframe. Thermal pion distribution function shown for comparison

Full size image
Fig. 3
figure 3

The decay pion \(\pi ^+\) 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]

Full size image

In Fig. 3 we plot the final pion spectra \(\pi ^+\) for a simple freeze-out surface with a constant Bjorken time, freeze-out temperature and radial velocity.Footnote 7 In addition to the total pion spectrum (which includes all decay chains producing \(\pi ^+\)), we also show the pion spectrum from the dominant decay channels \(\rho ^{+,-}\rightarrow \pi ^++\pi ^{0,-}\) and \(\omega ^0\rightarrow \pi ^++\pi ^-+\pi ^0\) (where \(\rho ^{+,0}\) and \(\omega ^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 procedure Eq. (8). The vector distribution function components \(f_i^\text {eq}\) shown in Fig. 2 only need to calculated once for a particular freeze-out temperature \(T_\text {fo}\) and then can applied to any shape of the freeze-out surface or the fluid velocity field \(u^\mu (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_i^\text {eq}\) will need to be tabulated not only in the fluid-restframe energy \(\bar{E}_\mathbf{p}=-u^\mu p_\mu \), 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.

Fig. 4
figure 4

Lorentz invariant pion \(\pi ^+\) weight functions \(f_i^\text {bulk}(-u^\mu p_\mu )\) and \(f_i^\text {shear}(-u^\mu p_\mu )\), Eqs. (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/|\mathbf{p}|^2\) was factored out from \(f_i^\text {shear}\) components. Bulk/shear perturbation of (initial) pion distribution function shown for comparison

Full size image

4 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 \(\pi ^{\mu \nu }(x)\) and the bulk-viscous pressure \(\Pi (x)\), so that the initial vector distribution function in the Cooper–Frye formula, Eq. (7), is

$$\begin{aligned} g^\mu (p,u,\pi , \Pi )= f(u^\mu (x), \pi ^{\rho \sigma }(x) , \Pi (x), p^\mu )p^\mu . \end{aligned}$$
(18)

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 \(\delta f\) corrections used in modern hydrodynamic simulations [32]

$$\begin{aligned} \delta f^\text {bulk}(\bar{E}_\mathbf{p}, \Pi ) =&f_\text {eq}(1\pm f_\text {eq})\left[ \frac{\bar{E}_p}{T}\left( \frac{1}{3}-c_s^2\right) \right. \nonumber \\&\left. -\frac{1}{3}\frac{m^2}{T\bar{E}_p} \right] \frac{\tau _\Pi \Pi }{\zeta }, \end{aligned}$$
(19)
$$\begin{aligned} \delta f^\text {shear}(\bar{E}_\mathbf{p}, \pi _{\rho \nu }p^\rho p^\nu ) =&f_\text {eq}(1\pm f_\text {eq})\frac{\pi _{\rho \nu }p^\rho p^\nu }{2(e+p)T^2}. \end{aligned}$$
(20)

Here \(c_s(T)\) is the speed of sound of the medium at the freeze-out temperature, m – mass of the primary resonance, and \(\tau _\Pi /\zeta \) 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

$$\begin{aligned} g^\mu = g^\mu _\text {eq}+ g_\text {shear}^{\mu \rho \sigma }\pi _{\rho \sigma }+g^\mu _\text {bulk} \Pi +\ldots \end{aligned}$$
(21)

where the derivatives \(g_\text {shear}^{\mu \rho \sigma }\equiv \partial g^{\mu }/\partial {\pi _{\rho \sigma }}\) and \(g_\text {bulk}^\mu \equiv \partial g^{\mu }/\partial \Pi \) can only be functions of 4-vectors \(p^\mu \) and \(\bar{E}_\mathbf{p}u^\mu \), and Lorentz scalars like temperature, chemical potential or resonance mass. Initially they are given by Eqs. (19) and (20)

$$\begin{aligned} g^\mu _\text {bulk} \Pi = \delta f^\text {bulk} p ^\mu , \quad g_\text {shear}^{\mu \rho \sigma }\pi _{\rho \sigma }= \delta f^\text {shear}p^\mu . \end{aligned}$$
(22)

After the decays \(g^\mu _\text {bulk}\) and \(g_\text {shear}^{\mu \rho \sigma }\) 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

$$\begin{aligned} g^{\mu }_\text {bulk}\Pi =&\Big [\left( p^\mu -\bar{E}_\mathbf{p}u^\mu \right) f^\text {bulk}_1(\bar{E}_\mathbf{p})\nonumber \\&+ \bar{E}_\mathbf{p}u^\mu f^\text {bulk}_2(\bar{E}_\mathbf{p})\Big ]\times \frac{-\tau _\pi \Pi }{\zeta }, \end{aligned}$$
(23)

and shear perturbations

$$\begin{aligned} g^{\mu \nu \rho }_\text {shear} \pi _{\nu \rho }&= \Big \{ [\eta ^{\rho \sigma }(p^\mu -\bar{E}_\mathbf{p}u^\mu )-\frac{2}{5} \eta ^{\rho \mu } \Delta ^{\sigma \alpha }p_\alpha ] f^\text {shear}_1(\bar{E}_\mathbf{p})\nonumber \\&\quad + \frac{2}{5}\eta ^{\rho \mu } \Delta ^{\sigma \alpha }p_\alpha f^\text {shear}_2(\bar{E}_\mathbf{p})\nonumber \\&\quad + \eta ^{\rho \sigma } \bar{E}_\mathbf{p}u^\mu f^\text {shear}_3(\bar{E}_\mathbf{p})\Big \}\times \frac{p^\nu \pi _{\nu \rho }p_\sigma }{2(e+p)T^2}. \end{aligned}$$
(24)

The bulk pressure perturbation does not introduce new tensor structures and the decomposition is the same as for the equilibrium distribution in Eq. (10), but the initial distribution functions \(f^\text {bulk}_i\) are, of course, different and can be read off from Eq. (19). The linear perturbations in the shear-stress tensor induces a rank-3 tensor distribution function \(g^{\mu \nu \rho }_{\text {shear}}\), which has three non-vanishing irreducible components \(f^\text {shear}_1\), \(f^\text {shear}_2\) and \(f^\text {shear}_3\) 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)

$$\begin{aligned} f_i^b(\bar{E}_\mathbf{p})&= B \frac{\nu _a}{\nu _b}\frac{m_a^2}{m_b^2 }\frac{1}{2}\int _{-1}^1 dw A_i(w) f_i^a(E(w)). \end{aligned}$$
(25)

where integration measures \(A_i(w)\) are listed in Appendix A.

In Fig. 4 we show the final decayed pion \(\pi ^+\) spectrum in the fluid-restframe due to viscous perturbations at \(T_\text {fo}=145\,\text {MeV}\).Footnote 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 fluid-restframe. Note that we factored in \(|\bar{\mathbf{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.

Fig. 5
figure 5

Lorentz invariant pion \(\pi ^+\) weight functions \(f_i^\text {temp.}(-u^\mu p_\mu )\) and \(f_i^\text {velocity}(-u^\mu p_\mu )\), Eqs. (23) and (23), which determine the final pion spectrum from temperature and velocity perturbations of the distribution functions of primary resonances, Eqs. (28) and (29). The explicit prefactors indicate the required freeze-out surface elements for the Cooper–Frye integration in the fluid-restframe. Note that a term \(\bar{T}/|\mathbf{p}|\) was factored out from \(f_i^\text {velocity}\) components. Temperature/velocity perturbation of (initial) pion distribution function shown for comparison

Full size image

Similarly to viscous perturbations, one can also consider linear perturbations of fluid velocity \(\delta u^\mu \), or temperature \(\delta T\) around the background fields \(\bar{u}^\mu \) and \(\bar{T}\),Footnote 9

$$\begin{aligned} \delta f^\text {temp.}(\bar{E}_\mathbf{p}, \delta T)&= f_\text {eq}(1\pm f_\text {eq})\frac{\bar{E}_\mathbf{p}}{\bar{T}}\frac{\delta T}{\bar{T}}, \end{aligned}$$
(26)
$$\begin{aligned} \delta f^\text {velocity}(\bar{E}_\mathbf{p}, \delta u_\mu p^\mu )&= f_\text {eq}(1\pm f_\text {eq})\frac{\delta u_\nu p^\nu }{\bar{T}}. \end{aligned}$$
(27)

The irreducible decomposition for temperature perturbation is the same as for the equilibrium distribution

$$\begin{aligned} g^{\mu }_\text {temp.} \delta T = \Big [&\left( p^\mu -\bar{E}_\mathbf{p}u^\mu \right) f_1^\text {temp.}(\bar{E}_\mathbf{p}) \nonumber \\&+ \bar{E}_\mathbf{p}u^\mu f_2^\text {temp.}(\bar{E}_\mathbf{p})\Big ]\times \frac{\delta T}{\bar{T}} . \end{aligned}$$
(28)

The irreducible decomposition for velocity perturbations consists of symmetric traceless tensor, vector and scalar representations with corresponding weight functions \(f_i^\text {velocity}\)

$$\begin{aligned} g^{\mu \nu }_\text {vel.}\delta u_\nu =&\Big \{[\eta ^{\nu \rho }(p^\mu -\bar{E}_\mathbf{p}u^\mu )-\frac{1}{3} \eta ^{\nu \mu }p_\sigma \Delta ^{\sigma \rho } )]f_1^\text {velocity}(\bar{E}_\mathbf{p})\nonumber \\&+\frac{1}{3} \eta ^{\nu \mu }p_\sigma \Delta ^{\sigma \rho } f_2^\text {velocity}(\bar{E}_\mathbf{p})\nonumber \\&+ \eta ^{\nu \rho }\bar{E}_\mathbf{p}u^\mu f_3^\text {velocity}(\bar{E}_\mathbf{p})\Big \}\times \frac{\delta u_{\nu } p_\rho }{\bar{T}}. \end{aligned}$$
(29)

The total vector distribution function is then given by

$$\begin{aligned} g^\mu _\text {ideal} = \bar{g}^\mu + g_\text {vel.}^{\mu \rho }\delta u_{\rho }+g^\mu _\text {temp.} \delta T+\ldots \end{aligned}$$
(30)

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^\mu \), 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 \(\mu _Q(x)\) and one could consider linear perturbations of the chemical potential \(\delta \mu _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^\mu =n_B u^\mu +j^\mu _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)

$$\begin{aligned} \delta f^\text {diffusion}(\bar{E}_\mathbf{p}, j^\mu _D p_\mu )=f_\text {eq}(1\pm f_\text {eq})\left[ \frac{n_B}{e+p}-\frac{Q_B}{\bar{E}_\mathbf{p}}\right] \frac{j_D^\nu p_\nu }{\hat{\kappa }}. \end{aligned}$$
(31)

where \(n_B\) is the local baryon density, \(Q_B\) is the baryon charge and \(\hat{\kappa }_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).

5 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\sigma _\mu \)

$$\begin{aligned} E_\mathbf{p}\frac{d N}{d^3 \mathbf{p}} = \frac{\nu }{(2\pi )^3} \int _\sigma d\sigma _\mu {\Big \{} F p^\mu + G \bar{E}_\mathbf{p}u^\mu + H \frac{p^\nu \pi _\nu ^{\;\;\mu }|\bar{\mathbf{p}}|^2}{2(e+p)T^2} {\Big \}}, \end{aligned}$$
(32)

where explicitly these terms are

$$\begin{aligned} F=&f^\text {eq}_1(\bar{E}_\mathbf{p}) + f_1^\text {shear}(\bar{E}_\mathbf{p}) \frac{\pi _{\rho \sigma }p^\rho p^\sigma }{2(e+p)T^2} + f_1^\text {bulk}(\bar{E}_\mathbf{p}) \frac{-\tau _\pi \Pi }{\zeta } ,\nonumber \\ G=&f^\text {eq}_2(\bar{E}_\mathbf{p})-f^{\text {eq}}_1(\bar{E}_\mathbf{p}) + \left( f_2^\text {bulk}(\bar{E}_\mathbf{p}) - f_1^\text {bulk}(\bar{E}_\mathbf{p}) \right) \frac{-\tau _\pi \Pi }{\zeta }\nonumber \\&+ \left( f_3^\text {shear}(\bar{E}_\mathbf{p}) - f_1^\text {shear}(\bar{E}_\mathbf{p})\right) \frac{\pi _{\rho \sigma }p^\rho p^\sigma }{2(e+p)T^2},\nonumber \\ H=&\left( f_2^\text {shear}(\bar{E}_\mathbf{p}) -f_1^\text {shear}(\bar{E}_\mathbf{p}) \right) \frac{2}{5} . \end{aligned}$$
(33)

The required values of fluid velocity \(u^\mu \), shear and bulk stresses \(\pi ^{\mu \nu }\), \(\Pi \) on the freeze-out surface must be provided by the hydrodynamic model of the QGP fireball or freeze-out 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^\mu = \bar{u}^\mu +\delta u^\mu \). In high energy collisions, a boost invariance is a good symmetry and the collision is often parametrized in Bjorken coordinates

$$\begin{aligned} ds^2=-d\tau ^2 + dr^2 + r^2 d\phi ^2 + \tau ^2 d\eta ^2. \end{aligned}$$
(34)

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 \(\tau \text {--}r\) plane, which can be conveniently parametrized by \((\tau (\alpha ), r(\alpha ) )\) where \(\alpha \in (0,1)\) and

$$\begin{aligned} d\sigma _\mu = \tau (\alpha ) r(\alpha ) \left( \frac{\partial r}{\partial \alpha }, -\frac{\partial \tau }{\partial \alpha }, 0, 0 \right) d\alpha d\phi d\eta \end{aligned}$$
(35)

Similarly by symmetry the fluid velocity has only two components and can be written in terms of a radial fluid rapidity \(\bar{\chi }\),

$$\begin{aligned} \bar{u}^\mu = \left( \cosh (\bar{\chi }), \sinh (\bar{\chi }), 0, 0 \right) . \end{aligned}$$
(36)

Note that then the particle energy in fluid-restframe is

$$\begin{aligned} \bar{E}_\mathbf{p}= m_T \cosh (\bar{\chi }) \cosh (\eta -y) - p_T \sinh (\bar{\chi }) \cos (\phi -\varphi ). \end{aligned}$$
(37)

where in the coordinate system of Eq. (34) and at space time point \((\tau ,r,\phi ,\eta )\) the particle momentum components are given by

$$\begin{aligned} \begin{aligned} p^{\tau }&=m_T \cosh (\eta -y),\quad p^{r}=p_T \cos (\phi -\varphi ), \\ p^{\phi }&=\frac{p_T}{r} \sin (\phi - \varphi ),\quad p^{\eta }=\frac{m_T}{\tau } \sinh (\eta -y). \end{aligned} \end{aligned}$$
(38)

Here we use the transverse momentum \(p_T\) (and transverse mass \(m_T=\sqrt{m^2+p_T^2}\)), the particle momentum angle \(\varphi \), 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 \(\bar{\pi }^\phi _\phi \) and \(\bar{\pi }^\eta _\eta \)

$$\begin{aligned} \begin{aligned} \bar{\pi }_{\tau \tau }&= -\bar{u}^r \bar{u}^r\left[ \bar{\pi }^\phi _\phi + \bar{\pi }^\eta _\eta \right] , \quad \bar{\pi }_{\tau r } = \bar{\pi }_{r\tau } =\bar{u}^r \bar{u}^\tau \left[ \bar{\pi }^\phi _\phi + \bar{\pi }^\eta _\eta \right] , \\ \bar{\pi }_{rr}&= -\bar{u}^\tau \bar{u}^\tau \left[ \bar{\pi }^\phi _\phi + \bar{\pi }^\eta _\eta \right] , \quad \bar{\pi }_{\phi \phi } = r^2 \bar{\pi }^\phi _\phi , \quad \bar{\pi }_{\eta \eta } = \tau ^2 \bar{\pi }^\eta _\eta , \end{aligned} \end{aligned}$$
(39)

and the bulk viscous pressure is simply \(\bar{\Pi }\).

The decay hadron spectrum for azimuthally and boost invariant freeze-out surface then reduces to a single integral

$$\begin{aligned} \frac{d N}{2\pi p_Tdp_Tdy}= & {} \frac{\nu }{(2\pi )^3} \int _0^1\! d\alpha \; \tau (\alpha ) r(\alpha )\nonumber \\&{\Bigg \{} \frac{\partial r}{\partial \alpha } \Big [ K^\text {eq}_1 + \frac{\bar{\pi }^\eta _\eta }{2(e+p)T^2} \; K_1^\text {shear}\nonumber \\&+ \frac{\bar{\pi }^\phi _\phi }{2(e+p)T^2} \; K_3^\text {shear}+ \frac{-\tau _\pi \bar{\Pi }}{\zeta }\; K_1^\text {bulk} \Big ]\nonumber \\&-\frac{\partial \tau }{\partial \alpha } \Big [ K^\text {eq}_2 + \frac{\bar{\pi }^\eta _\eta }{2(e+p)T^2} \; K_2^\text {shear}\nonumber \\&+ \frac{\bar{\pi }^\phi _\phi }{2(e+p)T^2} \; K_4^\text {shear} + \frac{-\tau _\pi \bar{\Pi }}{\zeta } \; K_2^\text {bulk} \Big ]{\Bigg \}}, \nonumber \\ \end{aligned}$$
(40)

were the freeze-out kernels \(K_i(p_T, \bar{\chi })\) are solely functions of the transverse particle momentum and radial fluid velocity \(\bar{u}^r=\sinh \bar{\chi }\), 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 (\(\tau (\alpha ), r(\alpha )\)) and radial fluid rapidity profile \(\bar{\chi }(\alpha )\). For the equilibrium components \(f_i^\text {eq}\) these kernels are given by the following rapidity and angle integrals

$$\begin{aligned} \begin{aligned} K^\text {eq}_1(p_T,\bar{\chi }) =&\int _0^{2\pi } \!d\phi \int _{-\infty }^\infty \! d\eta \left\{ f^\text {eq}_1(\bar{E}_\mathbf{p}) m_T \cosh (\eta -y) \right. \\&\left. + \left( f^\text {eq}_2(\bar{E}_\mathbf{p}) - f^\text {eq}_1(\bar{E}_\mathbf{p})\right) \bar{E}_\mathbf{p}\cosh (\bar{\chi }) \right\} , \\ K^\text {eq}_2(p_T,\bar{\chi }) =&\int _0^{2\pi }\! d\phi \int _{-\infty }^\infty \! d\eta \left\{ f^\text {eq}_1(\bar{E}_\mathbf{p}) p_T \cos (\phi -\varphi ) \right. \\&\left. + \left( f^\text {eq}_2(\bar{E}_\mathbf{p}) - f^\text {eq}_1(\bar{E}_\mathbf{p})\right) \bar{E}_\mathbf{p}\sinh (\bar{\chi }) \right\} . \\ \end{aligned} \end{aligned}$$
(41)

Recall that \(\bar{E}_\mathbf{p}\) also depends only on the differences of \(\eta -y\) and \(\phi -\varphi \), Eq. (37), so the dependence on momentum rapidity y and angle \(\varphi \) 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_1^\text {eq}(p_T,\bar{\chi }=\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 \(\delta u^\mu \) the corresponding perturbations of the distribution functions are given explicitly in Eq. (27). Because the decay operator Eq. (1) is linear, if \(\delta u^\mu \) is written as linear superposition of Fourier modes in azimuthal angle \(\phi \) and space-time rapidity \(\eta \)

$$\begin{aligned} \delta u^\mu (x) = \left( \tanh (\bar{\chi }) \delta u^r, \delta u^r, \delta u^\phi , \delta u^\eta \right) e^{im\phi + i k \eta }, \end{aligned}$$
(42)

one finds that the distribution of hadrons after kinetic freeze-out and resonance decays depends on momentum space azimuthal angle \(\varphi \) and momentum space rapidity y via the combination \(e^{im\varphi +iky}\) with the same azimuthal wavenumber m and rapidity wavenumber k. This is a direct consequence of \(\text {U}(1)\times \mathbb {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].

6 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_\text {chem}\) for each species i of long lived hadrons by introducing an appropriate chemical potentials \(\mu _i(T)\) for \(T<T_\text {chem}\). Subsequently, the kinetic freeze-out may take place at some lower temperature \(T_\text {kin}\). Because the primary resonance spectra are still described only by temperature \(T_\text {kin}\) and the chemical potentials \( \mu _i(T_\text {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.