*Communication* **Low-Momentum Pion Enhancement from Schematic Hadronization of a Gluon-Saturated Initial State**

#### **Elizaveta Nazarova 1,2 , Łukasz Juchnowski 2, David Blaschke 2,3,4,\* and Tobias Fischer <sup>2</sup>**


Received: 18 December 2018; Accepted: 6 March 2019; Published: 11 March 2019

**Abstract:** We study the particle production in the early stage of the ultrarelativistic heavy-ion collisions. To this end the Boltzmann kinetic equations for gluons and pions with elastic rescattering are considered together with a simple model for the parton-hadron conversion process (hadronisation). It is shown that the overpopulation of the gluon phase space in the initial state leads to an intermediate stage of Bose enhancement in the low-momentum gluon sector which due to the gluon-pion conversion process is then reflected in the final distribution function of pions. This pattern is very similar to the experimental finding of a low-momentum pion enhancement in the ALICE experiment at the CERN Large Hadron Collider (LHC). Relations to the thermal statistical model of hadron production and the phenomenon of thermal and chemical freeze-out are discussed in this context.

**Keywords:** Boltzmann equation; gluon saturation; pion enhancement; ALICE; LHC; thermalization; hadronization

#### **1. Introduction**

One of the issues which can be addressed by the kinetic approach is the question of a low-momentum pion enhancement in heavy ion collisions [1]. There are several solutions proposed to explain this effect as, e.g., the hadronization and freeze-out in a chemical non-equilibrium [2–4], the separate freeze-out for strange particles [5], Bose-Einstein condensate (BEC) of pions [6–10], established by elastic rescattering in the final stage [10,11]. However, none of them are commonly accepted yet [8]. We believe, an explanation linked to the presence of non-equilibrium physics and a precursor of pion condensation in heavy ion collisions should be the favorable one, especially after the recent analysis of particle correlations performed by the ALICE collaboration is showing a coherent fraction of charged *π*-meson emission that is reaching 23% [1,9]. Such formation of a Bose condensate is usually described by the introduction of additional non-equilibrium parameters to the statistical approach [10,12], see also [2,8,13].

An alternative scheme may rely on the Boltzmann kinetic equation for gluons and pions with elastic rescattering and a simple model for the parton-hadron conversion process (hadronisation). There are deep physical reasons for the non-equilibrium and pion condensation at the Large Hadron Collider (LHC). It can be due to fast expansion and overcooling of the quark-gluon plasma (QGP), or due to gluon condensation in the color glass condensate (CGC) initial state preceeding subsequent hadronization of the low-momentum gluons into low-momentum pions [8]. A scenario with an

initial state dominated by gluons which subsequently hadronize, eventually via a quarkless evolution through a first order phase transition, has recently been considered in Ref. [14].

In this short communication we investigate the idea that a certain oversaturation of the purely gluonic initial state could lead by elastic rescattering to a precursor of Bose condensation in the gluon sector in the form of a low-momentum gluon enhancement which, however, should be depopulated by the gluon-pion conversion process and thus appear as low-*p* pion enhancement in the pion sector. The gluon-pion conversion process is assumed with a constant matrix element which may be pictured as the local limit of a quark one-loop diagram for the case of large quark mass (quark confinement). We demonstrate the evolution of the coupled gluon and pion distribution functions in this case within a schematic model of coupled kinetic equations.

#### **2. Kinetic Equation Approach to Thermalization and Hadronization**

We start with the kinetic equation in the form of a Boltzmann-Nordheim equation, which for a single particle distribution function *<sup>f</sup>* <sup>=</sup> *<sup>f</sup>*(*x*,*p*, *<sup>t</sup>*) can be written as

$$\frac{\mathrm{d}f}{\mathrm{d}t} = \mathbb{C}[f] \,. \tag{1}$$

where

$$\frac{\mathbf{d}f}{\mathbf{d}t} = \frac{\partial f}{\partial \vec{\mathbf{x}}} \frac{\mathbf{d}\vec{\mathbf{x}}}{\mathbf{d}t} + \frac{\partial f}{\partial \vec{p}} \frac{\mathbf{d}\vec{p}}{\mathbf{d}t} + \frac{\partial f}{\partial t} \tag{2}$$

and *C*[ *f* ] represents the collision integral. In this study we restrict ourselves to the case of a uniform (*<sup>∂</sup> <sup>f</sup>* /*∂<sup>x</sup>* <sup>=</sup> 0) system in a non-expanding box ( *<sup>F</sup>* <sup>=</sup> <sup>d</sup>*p*/d*<sup>t</sup>* <sup>=</sup> 0), therefore only the explicit time-dependence remains: d*f* /d*t* = *∂ f* /*∂t* ≡ *∂<sup>t</sup> f* .

On the other hand, the collision integral for the 1 + 2 → 3 + 4 process is defined as:

$$\mathcal{L}\left[f(t,\vec{p}\_1)\right] = \frac{(2\pi)^4}{2\mathcal{E}\_1} \int \delta^4(\sum\_i P\_i) |\mathcal{M}|^2 F[f] \prod\_{k=2}^4 \frac{\mathbf{d}^3 \vec{p}\_k}{(2\pi)^3 2E\_k} \,\mathrm{}\tag{3}$$

so that the Equation (1) will take the following form:

$$\delta\partial\_{l}f(t,\vec{p}\_{1}) = \frac{(2\pi)^{4}}{2\mathcal{E}\_{1}} \int \mathcal{S}^{4} (\sum\_{i} P\_{i}) |\mathsf{M}|^{2} F[f] \prod\_{k=2}^{4} \frac{\mathsf{d}^{3} \vec{p}\_{k}}{(2\pi)^{3} 2E\_{k}} \,. \tag{4}$$

describing elastic scattering of the system of particles of one type, e.g., gluons. Here for the process <sup>1</sup> <sup>+</sup> <sup>2</sup> <sup>→</sup> <sup>3</sup> <sup>+</sup> 4 we define as *fi* the distribution function of particle *<sup>i</sup>* with 4-momentum *Pi* = (*Ei*,*pi*), |*M*| as the transition amplitude of the process, and *F*[ *f* ]=(1 + *f*1)(1 + *f*2)*f*<sup>3</sup> *f*<sup>4</sup> − *f*<sup>1</sup> *f*2(1 + *f*3)(1 + *f*4) represents the gain and loss terms in the collision integral. In the current study we consider the distribution function to be isotropic through the whole evolution. Moreover, the matrix elements of all the processes involved are taken to be constant:

$$|M|\_{12 \to 34} = \text{const},\tag{5}$$

following [11], where the case of a system of pions was considered. Albeit this work describes the academic study with only constant matrix elements, the ongoing project involving momentum- and angle-dependent transition amplitudes is discussed in the Section 4.

As we consider an isotropic, uniform, non-expanding system and constant matrix elements for the processes, and taking into account the 4-momentum conservation (*P*<sup>1</sup> + *P*<sup>2</sup> = *P*<sup>3</sup> + *P*4), the Equation (4) takes the form

$$\left|\partial\_{l}f(\varepsilon\_{1})\right| = \frac{\left|M\right|^{2}}{64\pi^{3}\varepsilon\_{1}} \int \int \mathrm{d}\varepsilon\_{3} \mathrm{d}\varepsilon\_{4} \mathrm{DF}[f] \,\,,\tag{6}$$

where D = min{*p*1, *p*2, *p*3, *p*4} and *pi* are now the radial components of the three-momenta. Details of the derivation are shown in the Appendix A. For future investigations it is helpful to rewrite the Equation (6) in terms of a momentum integration, as we would like to extend the approach to the angle-dependent collision integral, as well as to non-uniform systems. Therefore, in the current work we use the following formula:

$$\left|\partial\_t f(p\_1)\right| = \frac{\left|M\right|^2}{64\pi^3 \varepsilon\_1} \int \int \frac{p\_3 p\_4}{\varepsilon\_3 \varepsilon\_4} \mathrm{d}p\_3 \mathrm{d}p\_4 \mathrm{DF}[f] \,. \tag{7}$$

Obviously, elastic scattering is necessary but not sufficient to achieve low-*p* pion enhancement. The second required process which needs to be accounted for is hadronization. In this exploratory work we connect the gluon sector directly with the pion one. For such system there are three contributing channels: *ππ* → *ππ*, *gg* → *gg*, and *gg* ↔ *ππ*. Therefore at the end we have a coupled system of equations:

$$\begin{split} \frac{\partial f\_{\pi}}{\partial t}(t, \vec{p}\_{1}) &= \int \int \frac{|M\_{\pi \pi \to \pi \pi \pi}|^{2}}{64\pi^{3}\varepsilon\_{1}} \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}F[f\_{\pi}] \\ &+ (1 + f\_{\pi}(t, p\_{1})) \int \int \frac{|M\_{\mathrm{S}\Sigma} \rightarrow \pi\pi|}{64\pi^{3}\varepsilon\_{1}} \Big| \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}(1 + f\_{\pi}(t, p\_{2})) f\_{\mathcal{S}}(t, p\_{3}) f\_{\mathcal{S}}(t, p\_{4}) \\ &- f\_{\pi}(t, p\_{1}) \int \int \frac{|M\_{\mathrm{T}\pi \rightarrow \pi \frac{\pi}{2}}|^{2}}{64\pi^{3}\varepsilon\_{1}} \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}f\_{\pi}(t, p\_{2}) \left(1 + f\_{\mathcal{S}}(t, p\_{3})\right) \left(1 + f\_{\mathcal{S}}(t, p\_{4})\right) \end{split} \tag{8a}$$

$$\begin{split} \frac{\partial f\_{\mathcal{S}}}{\partial t}(t, \vec{p}\_{1}) &= \int \int \frac{|M\_{\mathcal{S}\mathcal{S}} - \boldsymbol{p}\_{\mathcal{S}}|^{2}}{64\pi^{3}\varepsilon\_{1}} \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}\mathcal{F}[f\_{\mathcal{S}}] \\ &+ \left(1 + f\_{\mathcal{S}}(t, p\_{1})\right) \int \int \frac{|M\_{\pi\pi\rightarrow\gamma\mathcal{S}\mathcal{S}}|^{2}}{64\pi^{3}\varepsilon\_{1}} \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}\left(1 + f\_{\mathcal{S}}(t, p\_{2})\right) f\_{\pi}(t, p\_{3}) f\_{\pi}(t, p\_{4}) \\ &- f\_{\mathcal{S}}(t, p\_{1}) \int \int \frac{|M\_{\mathcal{S}\mathcal{S}\rightarrow\pi\pi\pi\big|}|^{2}}{64\pi^{3}\varepsilon\_{1}} \frac{p\_{3}p\_{4}}{\varepsilon\_{3}\varepsilon\_{4}} \mathrm{d}p\_{3} \mathrm{d}p\_{4} \mathrm{D}f\_{\mathcal{S}}(t, p\_{2}) \left(1 + f\_{\pi}(t, p\_{3})\right) \left(1 + f\_{\pi}(t, p\_{4})\right) \end{split} \tag{8b}$$

where M*gg*→*ππ* and M*ππ*→*gg* are matrices for hadronization channels. Note, that due to the momentum conservation *p*<sup>2</sup> = *p*<sup>3</sup> + *p*<sup>4</sup> − *p*<sup>1</sup> in Equation (8). In this study we set M*ππ*→*gg* = 0, which is motivated by the threshold for this process due to the large value of the gluon mass: *mg* = 0.7 GeV. The value of M*gg*→*ππ* is set to be constant and should be seen as an academic example.

As the initial condition of the system we take an oversaturated gluon distribution given by a step-like function [15,16] inspired by the CGC picture of the initial state which is assumed to have no pions

$$\left.f\_{\pi}(t,p)\right|\_{t=0} = 0, \qquad \qquad \left.f\_{\S}(t,p)\right|\_{t=0} = f\_0\theta(1-p/Q\_s) \,. \tag{9}$$

By *Qs* we denote the saturation scale. However, in order to avoid numerical problems that would occur with the step-function distribution, we use instead the following smooth function [15]

$$f\_{\mathcal{S}}(t, p)\Big|\_{t=0} = f\_0 \left[\theta(1 - p/Q\_s) + \theta(p/Q\_s - 1)e^{-a\left(p/Q\_s - 1\right)^2}\right], \quad a = 10\tag{10}$$

to define the initial conditions.

We keep our model simple and therefore do not introduce an extra timescale for the start of hadronization. However, we keep in mind that the underlying microphysical process is, e.g., a quark-box diagram, which consists of the Breit-Wheeler type process of 2*g* → *qq*¯ and subsequent hadronization cross section *qq*¯ → *ππ*. In the future we plan to investigate the problem of the gluon-to-pion conversion in detail, for instance within a Nambu–Jona-Lasinio model [17–20] and/or by exploiting dynamical schemes of hadronization that would address the confinement aspect as well [21–24].

#### **3. Results**

In Figure 1 we show the evolution of the gluon distribution function from a CGC motivated initial (over-)saturated gluon state to a thermal distribution due to elastic scattering according to the *gg* → *gg* process. The timescale to reach a thermalized final state is of the order of *t*final ∼ 250 fm/c and thus exceeds the typical duration evolution towards freeze-out of the fireball created in a heavy-ion collision. This is mainly due to the fact that the value of the matrix element taken in this example calculation as |*M*| = 4.5 is unrealistically small.

**Figure 1.** (Color online) The evolution of gluon distribution function *f*(*p*) with time in a system of massive gluons (*m* = 0.7 GeV). The final distribution is shown as a bold red line, while the initial function is drawn as a bold black line. Thin black lines represent the intermediate stages of the gluon distribution function. The final time of the evolution represents the point when the pion distribution reaches equilibrium.

In Figure 2 we show the same evolution of the gluon distribution function for three different values of the matrix element. The value *M* = 140 leads to a thermalization time scale which nicely corresponds to the result of a calculation by Shuryak [25].

**Figure 2.** (Color online) Same as Figure 1 for different matrix elements |*M*| = 1, 4.5, 140.

When the coupling to the pion sector is switched on, the gluon conversion proceeds and the initially empty pion phase space gets populated at the expense of the gluon one. Due to the relation of the gluon and pion masses the reverse process (the pion annihilation to two gluons) does not practically take place. In Figure 3 the evolution from the initially pure gluon saturated state to the thermal pion state without gluons is shown. The pion distribution shows clearly the low-momentum enhancement typical for a precursor of Bose condensation. This is the fact observed in the ALICE experiment at CERN for which we wanted to give a qualitative explanation with the simple kinetic model presented here. It should be noted that here we used as a test the equal values for the three transition amplitudes: |*Mgg*→*gg*| = |*Mππ*→*ππ*| = |*Mgg*→*ππ*| = 4.5.

**Figure 3.** (Color online) The evolution of pion (*m* = 0.14 GeV) and gluon (*m* = 0.7 GeV) distribution functions *f*(*p*) with time in a coupled pion-gluon system. Blue line represents the initial gluon distribution, while the final distributions are shown as bold black and red lines for pion and gluon distribution functions, respectively. Thin black lines represent the intermediate stages of the distribution functions. The final time of the evolution represents the point when the pion distribution reaches equilibrium.

Our simplified model shows, under the assumption of gluon dominance in the initial state, the quarkless evolution of the system towards a pion gas with low-momentum pion enhancement as a precursor of Bose condensation. According to (8) both particle species (gluons and pions) undergo two main processes: conversion and elastic scattering. Both of them are responsible for low-momentum (low-*p*) pion enhancement.

The fist process turns *π*-mesons to gluons and vice versa and its rate is defined by two matrix elements *Mππ*→*gg*, *Mgg*→*ππ*, which in the simplest case considered here are constant numbers. Particle conversion can take place only when energy of incoming particles is at least equal to mass of outgoing ones. Consequently, in case of massless gluons kinematics restricts *gg* → *ππ* reaction to higher energy region of a spectrum, making the whole process slow.

The impact of the second process (elastic scattering ) is more subtle. It lowers momentum of particles through subsequent collisions, leading them to "pile-up" near zero momentum mode. The effect is especially strong for bosons due to the statistical factor (1 + *fi*) in (8) and allows pre-condensate formation even before thermalization. In normal circumstances for long enough times the distribution should become an equilibrium Bose function. However, in our model massive gluons undergo a complete conversion to pions before thermalization because *mg* <sup>&</sup>gt; *<sup>m</sup><sup>π</sup>* (*mg* <sup>=</sup> 0.7 GeV, while *m<sup>π</sup>* = 0.14 GeV).

#### **4. Discussion**

The present model, albeit quite simple, shows the formation of the pion condensation precursor emerging from an oversaturated purely gluonic state. The process takes place before the system reaches equilibrium. The model can be improved by the use of non-constant matrix elements and thus taking into account scattering angle in collision kinematics. Such an improved model would allow us to discuss the different scales and their evolution, e.g., the Debye scale, the UV and IR scale, see Refs. [7,11,15,26].

These improved matrix elements should also bear the confining aspects of gluon-gluon interactions which ultimately should be responsible for the absence of gluons from the final state. The assumption of a constant gluon mass, exceeding the value of the pion mass is a rather schematic realization of this concept which provides ample room for improvement. Here it would be beneficial to make a comparison with the study in Ref. [15], where a system of massless gluons undergoes the evolution due to elastic scattering with similar restrictions as used in the current paper. However, the Equation (7) will no longer be valid in the case of non-constant matrix elements and angle-dependence, and thus will need to be rederived.

Another room for advancement lies in direct handling of the kinetics of Bose condensation (see, e.g., Ref. [10]). One way to do that is the separation of the distribution function into two parts:

$$\bar{f}\_{\pi}(p) = f\_{\pi}(p) + (2\pi)^3 n\_c^{\pi} \delta(p) \tag{11}$$

$$\bar{f}\_{\mathcal{S}}(p) = f\_{\mathcal{S}}(p) + (2\pi)^3 n\_c^{\mathcal{S}} \delta(p) \tag{12}$$

where the first term represents the "gas" and the second describes BEC. This ansatz has been discussed for the oversaturated pion gas in Ref. [11] and recently also for the gluon plasma in Ref. [27]. We hope to achieve manifest energy and particle number conservation with such an improved formulation of the particle kinetics in the presence or precursory development of a Bose condensate in the system.

The model can be extended towards a more realistic description of a hadronizing gluon-dominated initial state in high-energy heavy-ion collisions by including more hadronic species as they are observed in those experiments in good agreement with the thermal statistical model [28]. This calls then for an extension of the collision integrals in our kinetic model to other classes of processes than just 2 → 2 processes as, e.g., the three-meson conversion to a baryon-antibaryon pair and its reverse [29].

Last, but not least we want to mention that the assumed absence of dynamical quarks is only a simplifying assumption. In an improved model, their kinetics shall be coupled to that of the gluons and all considered hadron species. Their absence in the final state shall be realised due to a confining mechanism. The one already tested in the framework of a kinetic theory is the Gribov-Zwanziger confinement realized by an infrared-divergent selfenergy [22–24]. We shall come back to these issues in a subsequent, more elaborate work on the subject.

**Author Contributions:** The original conceptualization was done by D.B. and T.F., while E.N. and L.J. are responsible for the development of methodology, software, the formal analysis and preparation of the original draft. All of the authors (E.N., L.J., D.B., T.F.) contributed to reviewing and editing of the manuscript.

**Funding:** This research was funded by the Polish Narodowe Centrum Nauki under grant number UMO-2014/15/ B/ST2/03752 (E.N., L.J., D.B.) and grant number UMO-2016/23/B/ST2/00720 (E.N.,T.F.).

**Acknowledgments:** We acknowledge discussions with Viktor Begun, Marcus Bleicher, Wojciech Florkowski, Carsten Greiner, Gerd Röpke, Ludwik Turko and Dmitry Voskresensky on the topic of the nonequilibrium nature of the pion distribution emerging from an oversaturated initial state. Further we thank Niels-Uwe Friedrich Bastian for his support in numerical and algebraic questions. Brent Harrison gave a nice presentation of his ongoing work with Andre Peshier on "Bose-Einstein Condensation from a Gluon Transport Equation" at the Symposium on "Nonequilibrium Phenomena in Strongly Correlated Systems" in Dubna, 16 April 2018, which inspired us to prepare this communication.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Collision Integral Derivation**

Here we will simplify the collision integral (3):

$$\mathbb{C}\_{1}[f] = \frac{(2\pi)^{4}}{2\mathcal{E}\_{1}} \int \delta^{4}(\sum\_{i} P\_{i}) |\mathsf{M}|^{2} F[f] \prod\_{k=2}^{4} \frac{\mathsf{d}^{3} \vec{p}\_{k}}{(2\pi)^{3} 2E\_{k}} \tag{A1}$$

Using the identity:

$$\delta^3(\sum \vec{p}\_i) = \int \exp\left(i(\vec{\lambda}, \vec{p}\_1 + \vec{p}\_2 - \vec{p}\_3 - \vec{p}\_4)\right) \cdot \frac{\mathbf{d}^3 \vec{\lambda}}{(2\pi)^3},\tag{A2}$$

and separating the angle integrations:

$$\mathbf{d}\,\vec{p}\_i = \mathbf{d}\,p\_i \mathbf{d}\,\cos\theta\_i p\_i^2 \mathbf{d}\,p\_i = \varepsilon\_i p\_i \mathbf{d}\Omega\_i \mathbf{d}\varepsilon\_i \tag{A3}$$

the integral takes the following form:

$$\mathbb{C}\_{1}[f] = \frac{\left|M\right|^{2}}{64\pi^{3}\varepsilon\_{1}} \int \delta(\varepsilon\_{1} + \varepsilon\_{2} - \varepsilon\_{3} - \varepsilon\_{4}) \mathrm{DF}[f] \mathrm{d}\varepsilon\_{3} \mathrm{d}\varepsilon\_{4} \mathrm{d}\varepsilon\_{2} \tag{A4}$$

where D is defined as follows:

$$\mathbf{D} = \frac{p\_2 p\_3 p\_4}{64\pi^5} \int \lambda^2 \mathrm{d}\lambda \int \epsilon^{i(\vec{p}\_1, \vec{\lambda})} \mathrm{d}\Omega\_\lambda \int \epsilon^{i(\vec{p}\_2, \vec{\lambda})} \mathrm{d}\Omega\_2 \int \epsilon^{i(\vec{p}\_3, \vec{\lambda})} \mathrm{d}\Omega\_3 \int \epsilon^{i(\vec{p}\_4, \vec{\lambda})} \mathrm{d}\Omega\_4. \tag{A5}$$

Taking into account that

$$\begin{split} \int \epsilon^{i(\vec{p}\_1, \vec{\lambda})} \mathrm{d}\Omega\_{\lambda} &= \int \epsilon^{i(p\_1 \lambda \cos \theta\_{\lambda})} \mathrm{d}\Omega\_{\lambda} = \int\_0^{2\pi} \mathrm{d}\varrho \int\_{-1}^1 \mathrm{d}\cos \theta \epsilon^{i(p\_1 \lambda \cos \theta\_{\lambda})} = \\ &= \frac{2\pi}{ip\_1 \lambda} \epsilon^{ip\_1 \lambda x} \Big|\_{x=-1}^{x=1} = \frac{2\pi}{p\_1 \lambda} \frac{\epsilon^{ip\_1 \lambda} - \epsilon^{-ip\_1 \lambda}}{2i} \cdot 2 = \frac{4\pi}{p\_1 \lambda} \sin(p\_1 \lambda) \end{split} \tag{A6}$$

we can rewrite D as:

$$\mathbf{D} = \frac{4}{\pi p\_1} \int \frac{\mathbf{d}\lambda}{\lambda^2} \sin(p\_1 \lambda) \sin(p\_2 \lambda) \sin(p\_3 \lambda) \sin(p\_4 \lambda) \tag{A7}$$

Using the Fourier transformation:

$$\begin{split} & -\sqrt{\frac{\pi}{2}} w \text{Sign}(w) = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} \frac{1}{\mathbf{x}^2} e^{iw\mathbf{x}} \mathbf{dx} = \\ & \frac{1}{\sqrt{2\pi}} \left( \int\_0^{\infty} \frac{1}{\mathbf{x}^2} e^{iw\mathbf{x}} \mathbf{dx} + \int\_0^{\infty} \frac{1}{(-\mathbf{x})^2} e^{-iw\mathbf{x}} \mathbf{dx} \right) = \frac{1}{\sqrt{2\pi}} \int\_0^{\infty} \frac{1}{\mathbf{x}^2} \left( e^{iw\mathbf{x}} + e^{-iw\mathbf{x}} \right) \mathbf{dx} \end{split} \tag{A8}$$

we can simplify the integral in the formula for D:

<sup>D</sup> <sup>=</sup> <sup>4</sup> *πp*<sup>1</sup> ∞ 0 d*λ λ*2 *<sup>e</sup>ip*1*<sup>λ</sup>* <sup>−</sup> *<sup>e</sup>*−*ip*1*<sup>λ</sup>* 2*i <sup>e</sup>ip*2*<sup>λ</sup>* <sup>−</sup> *<sup>e</sup>*−*ip*2*<sup>λ</sup>* 2*i <sup>e</sup>ip*3*<sup>λ</sup>* <sup>−</sup> *<sup>e</sup>*−*ip*3*<sup>λ</sup>* 2*i <sup>e</sup>ip*4*<sup>λ</sup>* <sup>−</sup> *<sup>e</sup>*−*ip*4*<sup>λ</sup>* <sup>2</sup>*<sup>i</sup>* <sup>=</sup> <sup>=</sup> <sup>1</sup> 4*πp*<sup>1</sup> ∞ 0 d*λ λ*2 (*e <sup>i</sup>λ*(*p*1+*p*2) <sup>−</sup> *<sup>e</sup> <sup>i</sup>λ*(*p*2−*p*1) − −*<sup>e</sup> <sup>i</sup>λ*(*p*1−*p*2) + *e iλ*(−*p*1−*p*2) )(*e iλ*(*p*3+*p*4) − − *e <sup>i</sup>λ*(*p*4−*p*3) <sup>−</sup> *<sup>e</sup> <sup>i</sup>λ*(*p*3−*p*4) + *e iλ*(−*p*3−*p*4) ) = <sup>=</sup> <sup>1</sup> 4*πp*<sup>1</sup> ∞ 0 d*λ λ*2 *e <sup>i</sup>λ*(*p*1+*p*2+*p*3+*p*4) + *e* <sup>−</sup>*iλ*(*p*1+*p*2+*p*3+*p*4) <sup>−</sup> *<sup>e</sup> iλ*(*p*3+*p*4+*p*2−*p*1) − *e* <sup>−</sup>*iλ*(*p*3+*p*4+*p*2−*p*1) <sup>−</sup> *<sup>e</sup> <sup>i</sup>λ*(*p*3+*p*4+*p*1−*p*2) <sup>−</sup> *<sup>e</sup>* <sup>−</sup>*iλ*(*p*3+*p*4+*p*1−*p*2) + *e iλ*(*p*3+*p*4−*p*1−*p*2) + + *e* <sup>−</sup>*iλ*(*p*3+*p*4−*p*1−*p*2) <sup>−</sup> *<sup>e</sup> <sup>i</sup>λ*(*p*1+*p*2+*p*4−*p*3) <sup>−</sup> *<sup>e</sup>* <sup>−</sup>*iλ*(*p*1+*p*2+*p*4−*p*3) + *e iλ*(*p*4−*p*3+*p*2−*p*1) + + *e* <sup>−</sup>*iλ*(*p*4−*p*3+*p*2−*p*1) + *e <sup>i</sup>λ*(*p*4−*p*3+*p*1−*p*2) + *e* <sup>−</sup>*iλ*(*p*4−*p*3+*p*1−*p*2) <sup>−</sup> *<sup>e</sup> iλ*(*p*4−*p*3−*p*1−*p*2) − *e* −*iλ*(*p*4−*p*3−*p*1−*p*2) = <sup>=</sup> <sup>1</sup> 4*πp*<sup>1</sup> (−*π*)(|*p*<sup>1</sup> + *p*<sup>2</sup> + *p*<sup>3</sup> + *p*4|−|*p*<sup>3</sup> + *p*<sup>4</sup> + *p*<sup>2</sup> − *p*1|−|*p*<sup>3</sup> + *p*<sup>4</sup> + *p*<sup>1</sup> − *p*2| + |*p*<sup>3</sup> + *p*<sup>4</sup> − *p*<sup>1</sup> − *p*2|−|*p*<sup>1</sup> + *p*<sup>2</sup> + *p*<sup>4</sup> − *p*3| + |*p*<sup>4</sup> − *p*<sup>3</sup> + *p*<sup>2</sup> − *p*1| + |*p*<sup>4</sup> − *p*<sup>3</sup> + *p*<sup>1</sup> − *p*2|−|*p*<sup>4</sup> − *p*<sup>3</sup> − *p*<sup>1</sup> − *p*2|) = <sup>=</sup> <sup>1</sup> 4*πp*<sup>1</sup> (−*π*)(−4 min{*p*1, *<sup>p</sup>*2, *<sup>p</sup>*3, *<sup>p</sup>*4}) = min{*p*1, *<sup>p</sup>*2, *<sup>p</sup>*3, *<sup>p</sup>*4} *p*1 (A9)

The last step in this Equation (A9) (changing to the minimum function between the four momenta) can be easily done by checking one of the possibilities—for example, the case when *<sup>p</sup>*<sup>1</sup> <sup>=</sup> min{*p*1, *<sup>p</sup>*2, *<sup>p</sup>*3, *<sup>p</sup>*4} (or *<sup>p</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>p</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>p</sup>*<sup>3</sup> <sup>&</sup>lt; *<sup>p</sup>*4). Taking into account the 4-momentum conservation: *P*<sup>1</sup> + *P*<sup>2</sup> = *P*<sup>3</sup> + *P*4, we get the final result:

$$\mathcal{K}[f(\varepsilon\_1)] = \frac{|M|^2}{64\pi^3 \varepsilon\_1} \int \int \mathrm{d}\varepsilon\_3 \mathrm{d}\varepsilon\_4 \mathrm{DF}[f] \tag{A10}$$

where *D* = <sup>1</sup> *<sup>p</sup>*<sup>1</sup> min{*p*1, *p*2, *p*3, *p*4}. In order to change the Formula (A10) to the integration over momentum, we can use the connection between energy and momentum in the relativistic case: *<sup>ε</sup>*<sup>2</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>+</sup> *<sup>m</sup>*<sup>2</sup> <sup>→</sup> *<sup>ε</sup>*d*<sup>ε</sup>* <sup>=</sup> *<sup>p</sup>*d*p*, so that the Equation (A10) takes form:

$$\mathbb{C}[f(p\_1)] = \frac{\left|M\right|^2}{64\pi^3\varepsilon\_1} \int \int \frac{p\_3p\_4}{\varepsilon\_3\varepsilon\_4} \mathrm{d}p\_3 \mathrm{d}p\_4 \mathrm{DF}[f] \,. \tag{A11}$$

#### **References**



c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
