Non-equilibrium photons from the bottom-up thermalization scenario
Oscar Garcia-Montero










garcia@thphys.uni-heidelberg.de



              
 
In this work, I calculate the p_⊥ resolved spectra for the three stages of the bottom-up scenario, which are comparable to the thermal contribution, particularly at higher values of the saturation scale Q_S^2. Analytical solutions are obtained by including a parametrization of scaling solutions from far-from-equilibrium classical statistical lattice simulations into a small angle kinetic rate. Furthermore, a theoretically motivated ansatz is used to account for near-collinear enhancement of the low-p_⊥ radiation. The system is phenomenologically constrained using the charge hadron multiplicities from LHC and RHIC as in previous parametric estimates and fair agreement with the data available for photons was found. I find that for this realistic set of parameters, the contribution from the thermalizing glasma dominates the excess photons.


Non-equilibrium photons from the bottom-up thermalization scenario
    Oscar Garcia-Montero  

 Institut für Theoretische Physik,
                     Universität Heidelberg,
 Philosophenweg 16,
                     69120 Heidelberg, Germany

    December 30, 2023
===============================================================================================================================================================================



§ INTRODUCTION



Direct photons are radiated throughout the evolution of a heavy ion collision (HIC) and, due to a lack of final-state interactions, can escape the medium virtually unscathed.  As a consequence, these probes are sensitive to the different stages of the rapidly expanding fireball. In small systems, such as p+p or p+A collisions, direct photons produced are mostly prompt photons, whose invariant yield can be calculated perturbatively <cit.> or using hybrid approaches <cit.> to account for nuclear modification factors. However, in collisions of large nuclear systems, direct photons in the small transverse momentum range exhibit exponential enhancement, commonly explained by  hydrodynamical models. 
 
  In addition to those findings, the transverse plane anisotropy of the photon multiplicities has been studied <cit.>, finding non-vanishing flow coefficients <cit.>. This anisotropy is thought to arise from the space-time evolution of the underlying medium. To compute such quantities, hydrodynamical quark-gluon-plasma (QGP) and transport models have been compared to the available data.  Unfortunately, the simultaneous reproduction of the yields and the photon flow coefficients, v_n,  has been out of reach for theoretical models <cit.>. This challenge has been named the direct photon puzzle <cit.>. 
  
  Nevertheless, in those calculations the pre-equilibrium physics of the medium is not accounted for, which leads to the introduction of several uncertainties including the initial conditions for the hydrodynamical evolution. Pre-equilibrium sourced photons are also omitted, and while 
 it is the traditional idea that, because of the small early-times space-time volumes, a pre-equilibrium contribution to the direct photon spectrum is negligible, new results seem to suggest that such a source may contribute on the same order of magnitude as the thermal stages <cit.>. In a novel estimate, <cit.>, thermal and Glasma total photon yields were computed parametrically using the bottom-up thermalization scenario by Baier, Mueller, Schiff and Son (BMSS) <cit.>. They were found to be comparable to thermal total yields. For this, a  phenomenological matching was performed to account for the energy scale Q_S in the system as well as for the thermalization time, τ_th and temperature, T_th.

While a full phenomenological simulation that links the initial stage of the collision with the onset of hydrodynamics is still out of reach, broad progress has been achieved to understand its dynamics. The initial stage of the collision is classical and highly non-linear in the gluon fields and after a parametrically short time, the evolution leads to instabilities which overpopulate the gluon fields <cit.>. 
Using classical statistical simulations, it was shown that an over-occupied Glasma approaches a non-thermal fixed point <cit.>, and by doing so loses its memory of the details about the initial conditions. In this simulations the system goes through a universal scaling regime, where the gluon distribution function behaves as 

    f_g(τ;  p_⊥, p_z)= 1/α_s τ^α f_S (p_⊥ τ^β,  p_z τ^γ) .
 

Here, f_S  is a time independent function, whose shape is given by non-perturbative physics of the theory. The exponents α=-2/3, β=0 and γ=-1/3 confirm the parametrical descriptions of the BMSS scenario <cit.>, thus identifying its approach to thermalization as the correct description of the expanding Glasma. 

In this paper, I use the 2↔ 2 kinetic photon rate to calculate the p_⊥ resolved spectra following the assumptions for the estimates of ref. <cit.>. This is done in the context of the bottom-up thermalization scenario. The rate is further simplified using a small-angle approximation.
For this calculation, the momentum dependent non-thermal distribution of quarks is needed, which is sampled via the hard g→ qq̅ approximation. This means that f_q ∼α_s f_g, where f_g  is taken to be the non-equilibrium scaling solution from eq. (<ref>). The non-equilibrium rates are enhanced using a bremsbremstrahlung ansatz analogous to the complete leading order (LO) thermal rate <cit.>. The different contributions from the Glasma, as well as the later thermal stage, will be compared, to establish their relative dominance in this model. Finally, I will fix the model's parameters phenomenologically to present a qualitative comparison with data. 

This paper is organized as follows. In section <ref> the reader can find a small account of the kinetic framework used, as well as the low-p_⊥ enhancement ansatz used for the non-equilibrium case. In section  <ref> I summarize the BMSS thermalization scenario and present the 2↔ 2 leading log (LL) results for the Glasma in this context. In sec. <ref>, I will review the parameter fixing of Ref. <cit.> and apply it to the p_⊥ resolved spectra. The main body of the results achieved in this work is given in sec. <ref>. Here, a comparison of the different contributions to direct photons will be presented for both the LL and the LO case. Furthermore, I will also compare photon production in the BMSS scenario with its early thermalization counterpart. Finally, a qualitative comparison to ALICE and PHENIX data will be presented.  This will be followed by conclusions and outlook. In the Appendices <ref> and <ref> I present the derivation of the pre-equilibrium and thermal yields, respectively. 



§ APPROXIMATE KINETIC DESCRIPTION


Following Ref. <cit.>, the photon rate for a thermalizing colored medium will be calculated using a kinetic description. The emission rate of an on-shell photon with three-momentum p=(p_x,p_y,p_z) at a space-time point X=(t, x,y,z) from two-to-two scatterings is given generally by <cit.>


    EdN /d^4X d^3p   =    1/2 (2 π)^12∫d^3 p_3/2E_3d^3 p_2/2E_2d^3 p_1/2E_1 |ℳ|^2
         ×  (2 π)^4 δ^4(P_1+P_2-P_3-P) 
         ×  f_1(p_1) f_2(p_2) [1± f_3(p_3)]  ,

with P_i=(E_i,p_i), i=1,2,3.
The total squared amplitude |ℳ|^2 is understood as summed over spins, colors and flavors of all in- and outgoing particles. For massless quarks, the annihilation of a quark-antiquark pair into a photon and a gluon yields the squared amplitude

    |ℳ_anni|^2 =160/9 16π^2 α α_s u^2+t̅^2/u t̅,

with the strong interaction coupling α_s and the electromagnetic coupling α. 
The squared amplitude for mixed Compton scattering, where a gluon kicks a (anti)quark producing a photon is  

    |ℳ_Comp|^2 =320/9 16π^2 α α_s u^2+s^2/-u s.

These are given in terms of the Mandelstam variables , t̅=(p_1-p)^2, and u=(p_3-p_1)^2.

To estimate the scattering rates, I will consider these two processes. Furthermore, the production of photons may be simplified using the small-angle approximation <cit.>. In this case, one finds the rate 

    EdN /d^4X d^3p = 40/9π^2 α α_S  ℒ   f_q (𝐩)  I_g   ,

where f_q is the quark distribution. The I_g integral is given by 

    I_g,q= ∫d^3p/(2π)^31/p f_g,q(p) .

The Coulomb logarithm ℒ in (<ref>) serves as a regulator and quantifies the ratio between the infrared and hard scales of the system, 

    ℒ= log(Λ_UV/Λ_IR)   .

As an example, one can take an isotropic thermal medium, where the hard scale Λ_UV is given by the temperature T, while the infrared scale for the kinetic description is defined by the screening mass, m_D ∼ g T, with α_s ≡ g^2/(4π). More generally, the Debye screening mass squared of the medium is estimated by m_D^2 = 4 g^2 (N_c I_g + N_f I_q).

The rate (<ref>) still neglects the resummation of multiple interactions with the medium that would contribute at leading order in α and α_s <cit.>. In thermal equilibrium the naive rate (<ref>) differs from the full leading-order result by about a factor of two in the relevant photon momentum range. For simplicity, and to also effectively take this into account, I adopt the prescription of Refs. <cit.> to replace the Coulomb logarithm in (<ref>) by 

    ℒ⟶ 2 log(1+ 2.912/g^2)

to match the leading log (LL) thermal result. This will be employed for all the small-angle estimates shown in this work. 




 §.§ Bremsstrahlung Ansatz




The small angle approximation gives the correct limit for thermal radiation at the LL level, once the coulomb logarithm has been identified as in eq. (<ref>). Nonetheless, it has been shown in ref. <cit.> that in a thermal medium, near-collinear bremsstrahlung dominates the rates for photon energies of p≲ 2 T, while at intermediate photon momenta, 2 T≲ p ≲ 10  T, the 2↔ 2 contribution is comparable to the near-collinear ones. 
One would expect this simple result to be the case in the non-equilibrium setting of the Glasma, with one general caveat. In the Glasma, the characteristic momentum scale is given by Q_s, making the near-collinear contributions during the early stages dominant at p≲ 2 Q_s which for RHIC and LHC energies covers most of the kinematic window at which excess has been observed (0.5-3 GeV).

Following this argument, the bremmstrahlung contribution has to be included, and we will do so by changing the total constant under the log in eq. (<ref>) where the temperature has to be substituted for the characteristic scale of the Glasma. Following the result from <cit.> the leading order (LO) thermal rate can be expressed as 



    Ed N /d^4 X d^3 p =A(p) ν(E/T) ,

where  

    A(p)= 2α d_F [ ∑_c q_c^2] m^2_D   f_q,eq(E/T) ,

with the thermal screening mass m^2_D=C_F  g^2_s  T^2 /4 and the Casimir operators   d_F=3 and C_F=4/3 for SU(3). The function ν is here the total constant under the log, and can be represented by different functions for both the LL and LO cases, 


    LL: ν(E/T) →ν_LL(E/T) →ℒ
       LO: ν(E/T) →ν_LL(E/T) + C_bremss(E/T)+C_anni(E/T) .


The explicit forms of C_bremss and C_anni can be found in ref. <cit.>.
In the non-equilibrium case, one can expand the former results by using the same LO function, while changing the temperature dependence for the appropriate characteristic scale of the system,  Q_s. Using this change, one can write down the function as follows 

    LL: ν(E/Q_s)→ν_LL(E/Q_s)  →ℒ 
       LO: ν(E/Q_s) →ν_LL(E/Q_s) + C_bremss(E/Q_s)+C_anni(E/Q_s) .


In the next sections, both the 2↔ 2 and the LO ansatz results will be presented to allow an appropriate comparison with the estimates presented by Berges et al <cit.>, as well as to show a better case scenario for photons coming from the BMSS scenario. For this, the non-equilibrium rates of Stages (i) and (ii) will be computed via eq. (<ref>) but would receive near-collinear enhancements thanks to the substitution from above.
In the case of the third stage, the rate is already given by the LL thermal rate, which means it can be just upgraded to the LO thermal rate, as it is parametrized in ref. <cit.>, for the appropriate space-time dependence of the temperature.


§ NON-EQUILIBRIUM SOURCES


The main assumption of the bottom-up scenario <cit.>, and also this work, is that gluonic saturation physics takes place for energies at LHC and RHIC. In this scheme, the non-thermal colored medium undergoes  three stages in its path for thermalization, which can be parametrically separated  as

    (i)   1≪ Q_s τ≪α_s^-3/2
       (ii)   α_s^-3/2≪ Q_s τ≪α_s^-5/2
       (iii)   α_s^-5/2≪ Q_s τ≪α_s^-13/5


In the beginning of the collision, at τ<Q_s^-1 the physics of the glue is highly  non-perturbative, and the states can be characterized by very non-linear macroscopic fields. During this stage, instabilities highly populate modes with p_⊥≲ Q_s <cit.>. After these modes have been occupied, the system is completely dominated by hard modes, for which  p_⊥∼ Q_s. 
These modes are approximately conserved, but due to Bjorken expansion their number density is diluted as n_h∼ Q_s^3/(Q_s τ). 
During this stage, gluons interact via 2↔ 2 hard scatterings, with a very small momentum exchange. This produces a broadening, or melting, of the distribution of longitudinal momentum p_z∼ Q_s (Q_s τ)^-1/3. The produced effect is a decrease of the typical occupation number as f_g ∼α_s^-1 (Q_s τ)^-2/3. 

The second stage starts when the typical occupation f_g falls below unity. This happens parametrically at Q_s τ∼α_s^-3/2.  During this stage, the number of soft gluons rise rapidly via collinear splitting.  Nonetheless, hard gluons still dominate the total number, with their number densities given by 

    n_h∼Q_s^3/Q_s τ  and   n_s∼α_s^1/4 Q_S^3/(Q_S τ)^1/2 .
  

In this stage, soft gluons possess a typical momentum of p_soft∼α^1/2 Q_s. It can be seen that n_soft/p_soft≫ n_hard/p_hard∼ n_hard/Q_s. This makes the Debye integral peaked strongly around the soft sector, from which m_D can be found to be m_D^2∼α_s  n_soft/p_soft∼α_s^3/4 Q_s^2  (Q_s τ)^-1/2. The typical longitudinal momentum of hard gluon stops decreasing in (ii) and converges to p_z∼α_s  Q_s, which means that the anisotropy of the system saturates at finite value.

The thermalization stage (iii) starts around the time Q_s τ∼α^-5/2, where n_h and n_s become comparable, while p_soft≪ p_hard. This signals that soft modes are now dominant both in number densities and in the screening mass. Soft gluons thermalize very fast via  2↔ 2 soft scatterings, and act as a thermal bath to which the hard sector looses energy via mini-jet quenching <cit.>. Since hard gluons act as a source of energy to the bath, the temperature rises with T=c_T α^3_s   Q_s^2 τ, to finally achieve full thermalization of the medium at 


    τ_th∼ c_eq α^-13/5 Q_s^-1  and   T_th∼ c_T c_eq α^2/5 Q_s .
  

Using this model to estimate the evolution of the Glasma, I can use these results to calculate the photon spectra produced on the road to thermalization. 

In the following, I use the variables τ=√(t^2 -z^2), η=arctanh(z/t) and y= arctanh (p_z/E). For the transverse plane, I employ a polar parametrization, p_x=p_⊥cosϕ and p_y=p_⊥sinϕ in terms of the transverse momentum p_⊥, longitudinal momentum p_z and azimuthal angle ϕ. The photon multiplicities will be obtained by integrating the eq. (<ref>) over the four-volume of the evolution for each stage, using d^4X= τdτdη d^2 x_⊥ and d^3p/E = dyd^2p_⊥. 



 §.§ Glasma, Stage I  


In the first stage of the evolution of the Glasma, quarks are taken to inherit their properties from the gauge sector via hard gluon splitting. This means that f_q∼α_s  f_g, where f_g is the gluon distribution found in classical statistical simulations. A small caveat has to be noted here regarding the quantum statistics of this function. This approximation is only valid while α_s  f_g≪ 1 which  will be the case for realistic parameters.
The distribution exhibits self-similarity, and during the scaling regime,  dynamics is given by a time-independent function f_S, from which one gets the gluon distribution via the relation 

    f_g(τ;  p_⊥, p_z)= 1/α_S(Qτ)^-2/3 f_S (p_⊥,  p_z (Qτ)^1/3) .


This scaling solution, f_S, was found in numerical studies for Bjorken expanding lattices <cit.>, and it is given by the form 

    f_S(p_⊥, p_z)=f_0 Q/p_⊥ exp[-1/2p_z^2/σ^2_0]  W_r[p_⊥-Q_s] ,

which was fitted from the results of ref. <cit.>,

where I define  W_r[p_⊥-Q_s] as the function that guarantees the suppression of the distribution function around p_⊥ = Q_s, as observed by simulations in <cit.>. The parametrization taken from the fits for W_r[p_⊥-Q_s] is as follows,

    = θ(Q_s - p_⊥) +  θ(p_⊥-Q_s )   
    e^-1/2(p_⊥-Q_s /r Q_s)^2 ,

where r is a free parameter  of 𝒪(1) that allows the correct suppression at higher momenta. During this stage, the gluon occupancy is dominated by hard gluons, which are approximately conserved, up to the expansion dilution factor τ^-1. This behaviour determines the time dependence of I_g, which can also be found by simply using the scaling properties in eq. (<ref>). An overall normalization constant, κ_g, is used as a  proportionality constant. This constant was found in ref. <cit.> to be given κ_g = c/(2 N_c) where c=1.1 is the gluon liberation coefficient <cit.>. Having this together, the gluon integral as 


    I_g(τ)=Q^2_s/4 π^2α_sκ_g/(Q_s τ) .


Using the aforementioned ingredients, the rate for the first stage is


    EdN /d^4 x d^3  p = 10/9π^4 α ℒ κ_g Q^2_s/(Q_s τ)  f_q (𝐩) .


For the computation of the photon multiplicity of the first stage, I assume that the photon momentum p is on shell, which means p_0= p_⊥cosh(y-η), as well as  p_z = p_⊥sinh(y-η). Integrating over the spatial variables, I get the expression 

    1/S_⊥dN_γ/d^2 p_⊥dy =  √(π/2)20/9π^4α_e κ_g     f_0 σ_0/Q_s log(τ_1/τ_0) 
       × Q^2_s/p^2_⊥  W_r[p_⊥-Q_s]  ,

using the limit of σ_0/Q_s → 0 [The details of the derivation of the Glasma yields can be found in Appendix A.  The derivation of the thermal QGP yield can be found in Appendix B. ].
The initial and final times are then substituted for τ_0 = c_0  Q_s^-1 and τ_1= c_1 α_s^-3/2 Q_s^-1, where c_0 and c_1 are unknown proportionality constants of order 𝒪(1). As it was stated before, a precise description of the Glasma will include such coefficients, but it is out of the scope of this work. In this model, these coefficients will be set to unity. This is in fact, supported by slow, logarithmic dependence of the coefficient ratio c_1/c_0  in eq. (<ref>).
  
The expansion in eq. (<ref>) is safe as long as one is interested in the p_⊥>σ_0 portion of the spectrum. Experimentally, this is the case, as we will be interested in the kinematic window p_⊥≥ 1 GeV,  and given that in RHIC and LHC, the characteristic saturation scales are thought to be  Q_s≳ 1 GeV for heavy ions, we can see that the interesting momenta for this setup satisfy this condition, for σ_0≲ 0.1 Q_s. Finally one can further integrate the multiplicities over the radial momenta to find the multiplicity of photons per unit rapidity. In the totally anisotropic limit, i.e. vanishing σ_0/Q_s, the results can be simplified to 


    .1/S_⊥dN_γ/dy|_p_⊥<Q_s =   √(π/2)20/9 π^3α_e κ_g f_0 σ_0  Q_s  
       ×log(α_s^-3/2) log(Q_s^2 /σ_0^2)  .


Which is identical to the result in ref. <cit.>  after performing matching the normalisation of the distribution function , f_0 in terms of the anisotropy parameter,  σ_0/Q_s,


    f_0 σ_0/Q_s = √(2/π)κ_q [log(Q_s^2/σ_0^2) ]  ^-1  ,


 I will use this matching  from now on this work. The rest of the yield can be found by integrating eq. (<ref>) above the saturation scale Q_s, which gives the result 

    .1/S_⊥dN_γ/dy|_p_⊥<Q_s =   √(π/2)40/9 π^3α_e κ_g κ_q  Q_s^2  χ_r  
       ×log(α_s^-3/2) [log(Q_s^2/σ_0^2) ]  ^-1 ,

where 

    χ_r = 1/2 e^-1/2 r^2(πerfi(1/√(2)
       r)-Ei(1/2 r^2)) .

In the strict limit of full anisotropy, this contribution vanishes. For the set of parameters used in this work it will contribute to around 5% of the total yield from the Glasma. From these results it can be seen that the total yield of photons of the first stage of the Glasma is basically insensitive to the fit parameters used in eq. (<ref>). 



 §.§ Glasma, Stage 2 


After a time τ∼ Q_s^-1 α_s^-3/2, the typical gluon occupation drops below unity, and the rate should be revised. At this point, hard gluons still dominate the total number density, but soft modes take over the behaviour of the Debye mass. This change affects the time dependence of the I_g integral, which behaves as I_g∼α_s^-1 m_D^2 <cit.>. This leads to the expression


    I_g(τ)=κ_g α_s^-1/4/4 π^2 √(c_1)Q^2_s/(Q_s τ)^1/2 ,

where the overall normalization of I_g has been modified to match the expression of stage (i) at τ∼ Q_s^-1 α_s^-3/2. The fermionic sector is always dominated by hard quarks, which can still be described by f_q=α_s  f_g, with f_g as in eqs. (<ref>) and (<ref>). With those changes, the full rate for stage (ii) now is given by 

    EdN /d^4 x d^3  p = 10/9π^4 αα_s^3/4 ℒ κ_g/√(c_1) Q^2_s/(Q_s τ)^1/2  f_q (𝐩) .


The photon multiplicity for this stage can be found by again integrating over the full space-time volume and expanding to leading order in the anisotropy parameter, to get  

    1/S_⊥dN^(ii)_γ/d^2 p_⊥dy =    √(2 π)20/9π^4α_e κ_g f_0 σ_0/Q_s Q^2_s/p^2_⊥ 
       × W_r[p_⊥-Q_s]  (√(c_2/c_1)α_s^-1/2-1) .



In the weak coupling limit, the system stays in (ii) a parametrically long time, which naturally leads to a higher yield than the stage (i). Here, as in the first stage, the dependence on the c_2/c_1 ratio is also slow and, since we expect these coefficients to be of order 𝒪 (1), it will be taken to be unity. 

The known result for the total yield can be found by again expanding in σ_0/Q_s up to lowest order. Applying the matching from eq. (<ref>) I get


    .1/S_⊥ Q^2_SdN_γ/dy|_p_⊥<Q_s = √(π/2)40/9 π^3   α_e ℒ κ_g κ_q
       × (α_s^-1/2-1) .


This result is identical to the one found the previous estimate <cit.>. The yield for photons with  p_⊥> Q_s  can be found by integrating eq. (<ref>), 


    .1/S_⊥ Q^2_SdN_γ/dy|_p_⊥>Q_s =  √(π/2)80/9π^3   α_e ℒ κ_g   f_0 σ_0/Q_s
       ×χ_r   (α_s^-1/2-1) ,

where  χ_r was defined in eq. (<ref>). Once again, it can be observed that the total yield exhibits independence to the fit parameters, just as in stage (i).



 §.§ Glasma, Stage 3 

 
During this stage, soft gluons take over the number density. The beginning of this stage is at time parametrically larger than the relaxation time, τ_2>τ_rel, which means that soft gluons modes  have thermalized already. Therefore, the distribution can be calculated using the rapidity integrated thermal rate (see Appendix <ref>), which can be further integrated in time to find the invariant yield, 


    1/S_⊥dN_γ/d^2 p_⊥dy  =   5/9  /π^2 ℒ∫_τ_2^τ_thτ τ  T^2(τ)K_0(/T(τ))   .



In this stage, the temperature of the soft gluon bath increases linearly in time, and can be parametrized as

    T=c_T α_s^3 Q_s^2 τ

Substituting the time parameter τ in the thermal rate integral, I get

    1/S_⊥dN_γ/d^2 p_⊥dy  =   5/9  ^-5/π^2 c_T^2   ℒ/Q^4∫_T_2^T_th T  T^3  K_0(/T)   .



Once  again it is important to notice that the parametrisation used in this equation works only for values of the transverse momentum such that p_⊥≳ T. However, for the kinematic window of interest, p_⊥≥ 1 GeV, it is also parametrically satisfied that p_⊥ >  T_th, which  means that the rate can be approximated using the asymptotic form 

    K_0(x) ∼√(π/2 x) e^-x  for   x ≫ 1  .
 

Throughout this stage the temperature rises, which makes the assumption safer, since T<T_th< p_⊥. In the next section I will show that for the photon kinematic window and for realistic parameters this is always the case. Using this approximation, the integral for the photon invariant yield in the third stage can be found to be 


    1/S_⊥dN_γ/d^2 p_⊥dy  =√(π/2) 5/9  ^-5/π^2 c_T^2  p_⊥^4/Q^432/945 ℒ( √(π)erf(√(p_⊥/T)).   
      . +√(T/p_⊥) e^-p/T(1
       - 1/2T/p_⊥+3/4T^2/p^2_⊥
       -.15/8T^3/p_⊥^3+105/16T^4/p_⊥^4))|_T_2^T_th .





§ PARAMETRICAL FIXING AND COMPARISON


The calculations presented above depend on a collection of parameters which describe the non-equilibrium dynamics of the BMSS scenario. To reduce the parameter phase-space, one can use a phenomenological matching to pin down some of these parameters. The first coefficient to fix is the semi-hard scale, Q_s, from which all the BMSS dynamics is dependent. It can be found in terms of the energy and N_part using IP-Glasma model <cit.>. The IP-Glasma model combines the IP-Sat model  <cit.>  with the MC-Glauber model dependence on the geometry of the system  <cit.>.  Since only the mixed quantity ⟨ S_⊥  Q^2_s ⟩≡∫ d^2x_⊥ Q_s(x_⊥, √(s)) can be calculated in terms of N_part, I will approximate the transverse average Q_s(x,N_part) by using 

    Q^2_s(x,N_part)=⟨ S_⊥  Q^2_s ⟩/⟨ S_⊥⟩ ,

where ⟨ S_⊥⟩ is calculated using the Glauber model. 

I have adjusted the overall normalization to a reference Q_s^2 = 2 GeV for the highest centrality class, 0-5%, with N_part=353. This can be done because we lack a first-principles determination of the hard scale in the Glasma, which makes it possible to vary the definition while keeping in mind it should be a semihard scale. 



  The next step is to relate the coefficients c_eq and c_T, used in sec. <ref>  and appendix  <ref> to measured quantities. For this the entropy per unit rapidity of QGP was matched to the entropy of produced charged hadrons. For an ideally expanding fluid, entropy is conserved, and one can perform the matching at thermalization time and temperature. The calculation is given in ref. <cit.>, and results in the expression 
 
    c_eq c_T^3/4=[ 45/148π^2k_S/Nα^7/5N_part/Q^2_s S_⊥2/N_partdN_ch/dy]^1/4 .
 

Here, the experimental input is the multiplicity of charged hadrons per unit rapidity, dN_ch/dy,   in terms of N_part, given in ref. <cit.>. The parameter k_S/N, is the proportionality constant that links total entropy of the hadronic phase with the measured multiplicity of charged hadrons. I will adopt the value k_S/N= 7.2  which has been extracted from particle spectra and particle interferometry <cit.>. I have adopted a running coupling parametrisation given by 

    α_s (Q_s) = 12π/(33-2N_f)log(Q_s^2/Λ_QCD^2) ,
 
where I take Λ_QCD=0.2 GeV, and the number of in-medium quarks is taken to be N_f=3. Using this data I get the phenomenologically interesting parameters listed in table <ref>, where Q_s, S_⊥, and c_eq c_T^3/4 are given for different centrality classes, that is N_part, for RHIC and LHC energies. The quantities c_eq and c_T can only be fixed together, as specified in eq. (<ref>). However, the thermalization temperature coefficient was found in ref. <cit.> up to logarithmic accuracy to be c_T = 0.18, and will vary it over an overall factor of 2. 

An important caveat is that, using the expressions for τ_th and T_th, as well for the temperature evolution in a Bjorken expansion, one can find also a phenomenological matching for τ_c, which is given by 

    τ_c= 45/74π^2k_S/N1/S_⊥dN_ch/dη1/T_c .


This renders τ_c insensitive to Q_s and α_s, while dependent uniquely on the number of participants in the collisions, N_part. Because, for N_part≲150 I find that τ_c≲τ_eq I will only take on account events which lie inside the 0-20% centrality range. 



§ RESULTS


In the following, the comparisons will be made by including all the sources for direct photons. That is, apart from the BMSS and QGP radiation, one should include also thermally  produced photons from a hadronic gas (HG) after the hadronization of the QGP phase. Photons produced by hard scattering and annihilation of the participating partons, referred normally as prompt photons have to be also included. The total direct invariant yield is then given by 


    dN /d^2 p_⊥dy= T_AA dσ^pp/d^2 p_⊥dy + K^γ[dN^gl /d^2 p_⊥dy +dN^th/d^2 p_⊥dy] .

where the hadronic  the σ^pp label stands for prompt photon (pQCD) cross section <cit.>, which will be scaled by a centrality dependent factor, T_AA, with A= Au , Pb. This can be calculated directly from the Glauber model. For the extension of the scaled pQCD  to lower p_⊥ values, which is needed to sum this contribution to the in-medium spectra, I have used the following functional form, 

    dσ^pp/d^2 p_⊥dy =  A_pp (1+p^2_⊥/P_0)^-n ,

which was used by PHENIX to fit the p+p results in <cit.>.
In eq. (<ref>), hadronic bremsstrahlung of photons is included as a thermal contribution, and summed over the QGP rate. In this work I included meson and baryonic rates, the π π bremsstrahlung rate, as well as reactions of the π ρ ω  meson system. <cit.>. This is calculated naively for a Bjorken expanding system by switching from the QGP to the HG rate at a freezout temperature of 154 MeV. A mild variation of this parameter did not affect the result. 

To account for the Glasma radiation, I have used σ_0/Q_s = 0.15 where a variation of 50% has been included as error bands. It can be seen parametrically from eqs. (<ref>) and (<ref>) that the dependence on this parameter is parametrically slow. The suppression parameter is chosen to be r=0.35 for a better fit to data. It is always important to note that even when both r and σ_0 are parameters physically motivated by ab-initio calculations, the particular value is picked as free parameter with soft constraints. This is indeed a source of systematic error in the model, but solving this issue with better simulations is outside the scope of this paper. Besides the gluon fitting parameters, other sources of uncertainty come from the thermalization time and temperature. These quantities vary with the constants c_T  and c_eq. As mentioned before, the product c_eq c_T^3/4  has been fixed together using eq. (<ref>), which still leaves freedom to vary c_T around the value calculated in ref. <cit.>. This will have an effect on the thermal invariant yield per transverse area since it is given parametrically by  τ_th^2 T_th^2 parametrically,  which means that after fixing it still depends polynomially on c_T^-10/3. For changes of a factor of 2 in c_T, that still amounts for an extra 𝒪(1) coefficient when one takes on account the changes in the temperature. The Glasma stages (i) and (ii) are also mildly sensitive to the coefficients c_0, c_1 and c_2, which can result in the rates picking up also extra order  𝒪(1) factors. Nevertheless, such fine tuning of the transition times is not the objective of this discussion.
 

I have also included the fitting prefactors K^γ to account for the total uncertainties of the parameters which may amount for overall normalization discrepancies, such as the initial volume and its evolution, as well as uncertainties in the spatial dependence of initial energy density and subsequent translation to the thermalization temperature and time. It can be seen in eqs. (<ref>), (<ref>), as well as for the thermal rate, eq. (<ref>), that these rates could in principle be sensitive to the inclusion of fluctuations of the spatial profile of Q_s. This factor is found to be K^γ =2.8. The result of this fitting can be seen in fig. <ref>, where I get fair agreement, particularly for ALICE data. 




In the previous estimates <cit.>, the difference between the thermal and Glasma total yields was found to be of order 𝒪(1)
In this work I confirm analytically these results (see secs <ref> and appendix <ref>) from the fit, regardless of the choice of parameters, once the matching of the constants has been enforced. Nonetheless,  in the p_⊥ differential result, a dominant structure was found for Glasma photons.  I want to make an emphasis on the fact that this structure is a signal from the thermalization process. Although its particular shape is model dependent as it strongly depends on the quark distribution, a more refined calculation from kinetic theory would give a reshuffling, or stretching of this yields, giving a more exponential look, while keeping the same order of magnitude around Q_s, since the small angle approximation is a good approximation for the hard scatterings of the 2↔ 2 processes in the Glasma. 

For the LL results, this structure is dominant at higher energies, and it is peaked at p_⊥∼ Q_s, while the thermal case strongly dominates at small energies, as can be seen in fig. <ref>. However, once enhanced to the LO rate this is not the case, as Glasma photons completely dominate over the thermal rates. In a realistic simulation, this may become less apparent, as the radial flow from the hydrodynamical expansion will blue-shift the thermal spectrum <cit.>, changing its slope, as well as enhancing the number of photons in this kinematic window. 



 
Apart from a clear-cut comparison between the stages of the evolution of the fireball, the LO version of fig <ref> allows us also to see the BMSS scenario in action, specifically from the curves for the first and second stages of the Glasma. The direct photon fraction for the first stage exhibits, in the very low-p_⊥ limit, a dip compared to its second stage counterpart. This comes directly from the screening mass time-dependence given in eq. (<ref>). We can see that the rise of soft gluons impacts the system such that it enhances the production of low-p_⊥ photons. In a violently evolving Glasma, the overoccupied gluons enhance the number of scatterings, which gives the conditions for photon production. 

As it can be seen in fig. <ref> and was stated above, in this model the direct photon contribution is dominated by the Glasma, while the thermal contribution is relevant only in the deeper infrared part of the kinematic window p_⊥≲ Q_s. This happens because of the BMSS estimate of τ_th and T_th, which gives a late and quite colder thermalization, both in RHIC and LHC. One has to remember that in more refined calculations, this may change, since the BMSS scenario poses as a upper bound for the thermalization time. Nonetheless, the calculation here serves as a proof of concept, showing that photons are extremely relevant and may serve a as extra source which may help solve the photon puzzle. This has to be contrasted, of course, with current phenomenological ideas, regarding the extra photons to come from later times <cit.>. A comparison as such stresses the need for higher correlation functions, especially those which are more sensitive to the space time evolution. Photon correlations, i.e. Hanbury-Brown-Twiss correlations <cit.>, may be able to help us discriminate between the two different scenarios[A manuscript on photon HBT as a way to discern between scenarios is already in preparation <cit.>.].


Aside from the comparisons with the invariant yield from fig.<ref>, I want to compare also to the total number of photons, in a way which total normalization drops off. For this, I will use the
total yield per unit rapidity, taken as a sum of all momenta, with an infrared cutoff of p_0. This quantity in the language of this work is given by  

    N_p_0/ y≡ 2 π ∫_p_0^∞ p_⊥  p_⊥ dN_γ/d^2 p_⊥dy


and will be normalized by its minimum bias counterpart.  For this observable at PHENIX (see fig. <ref>) I observe excellent agreement with the data except for deviations from this trend at p_0=1.2GeV. Furthermore, I checked numerically that for the LO case this observable is completely independent of the choosing r, within reasonable ranges given softly by constrain of the data from <cit.>.





 §.§ Early Hydro



For the phenomenologically matched system the bottom-up scenario gives a thermalization time and temperature in the order of τ_th∼ 2 fm amd T_th∼ 0.2 GeV. This contrasts with how normally hydrodynamical simulations are initialized, with initial times down to τ_i ∼ 0.4-0.6 fm and average initial temperatures as high as ⟨ T⟩_i ∼ 0.4-0.6 GeV. In this section we will compare the Glasma spectrum with a system in which thermalization occurs at Q_s τ_0∼ 1. For this comparison, we will call Early thermalization the integrated QGP rate from τ_0 until τ_th, for a Bjorken expanding fireball. For the results of this section, the hadronic bremsstrahlung is not taken on account, since for this discussion one is mostly interested in what happens before the system arrives to the freeze-out temperature. 

In fig. <ref> the reader can see a comparison of the early thermalization scenario with its Glasma counterpart. As a reference, I also give a curve for the late contribution from the thermal QGP. It is immediately apparent that in all results, the structure at p_⊥∼ Q_s dominates, and it is higher than the early hydro contribution. I do not expect the slope of the latter contribution to change significantly with a full-fledged simulation, since the initial conditions normally assigned require the pressure gradients to build up to start the expansion, which will result in the eventual blueshift of the spectrum. 





In the case of the LL results, the three contributions deliver similar results for the total photon yield, eq. (<ref>), as they have similar dependences at the infrared level. This is, however, not the case for the collinearly enhanced result (LO) in which, even though the early hydro and Glasma contributions are still of the same order, the Glasma dominates. This case is particularly strong for higher energies, in which the results of this paper seem to suggest that thermalizing matter is more relevant for the production of photons. In a higher energy system like Pb-Pb at √(s)=2,76 TeV in ALICE, gluons are generally speaking, harder. Apart from increasing the overall order of the radiation, this means that all the system will be shifted towards higher energies, which allows for the near-collinear radiation to be greatly enhanced for the Glasma. 

 On a final note, it is commonly argued <cit.> that the hydro expansion has to be initialized early on for the pressure gradients to build up, which will cause the blueshift expected from radial flow, as well to account for the transverse anisotropy coefficients, v_n. However, this assumption lays on top of a specific pick of initial conditions for the evolution, starting from an initial energy density, without any pre-equilibrium evolution. For a better understanding, a study of the fluctuations of the non-thermal evolution transverse pressure has to be done to start shedding light into this matter. Furthermore, a comprehensive matching of the resulting effective kinetic theory would give the correct initial conditions for hydro evolution, which, conceivably, would start from non-zero transverse velocities. Thus, the assumption that a late thermalization scenario cannot account for the anisotropy coefficients has to be checked via real-time lattice simulations, as it may not be correct. 
 




§ CONCLUSIONS



In this work, I derived the photon invariant yield from the Glasma using the  bottom-up thermalization framework. My results expand on the parametric estimates for the total photon yields from ref. <cit.> by obtaining analytically a transverse momentum resolution for the leading log spectra. In this scenario, a sizeable contribution coming from early times was found, which dominates the total direct photon spectra for p_⊥∼ Q_s and realistic parameters at LHC and RHIC energies. The LL results fail, however, to properly describe the low p_⊥ window of the measured spectra. The aforementioned results where improved by using the complete LO thermal rates, as well as introducing a bremsstrahlung modification ansatz for the non-equilibrium rates. 
With this ansatz, the non-thermal contribution becomes mostly dominant and a fair comparison to data can be found up to the overall normalization. To account for this uncertainty an order 𝒪(1) parameter is introduced. 

Apart from the overall normalization, this model is very sensitive to the specific shape of the quark distribution. This uncertainty is just an extension of weak coupling methods to assess collective QCD at RHIC and LHC energies. Nonetheless, this simple model serves as a proof of concept to challenge the idea that pre-equilibrium photon production is suppressed. Nonetheless, as was stated in the text, it was found that the total number of photons is affected only by the assumption of scaling in the quark distribution and not by the fit parameters σ_0 and r, and thus both remain free -yet physically motivated- parameters. This issue may be solved by more thorough calculations, in the context of classical statistical simulations or kinetic theory <cit.>. Even if weak coupling assumption has to be kept, this model can be improved by directly fitting  the quark distributions from such more involved calculations. 

As it was stated before, the inclusion of anisotropies is beyond the scope of this work, as I have taken all the distributions to be isotropic and homogeneous in the transverse plane. This approximation, while useful for building this model, relies only on the assumption that flow is suppressed at early times, and its accuracy is not clear in the extreme non-equilibrium setting of the early stages of a HIC. The inclusion of such anisotropies can in principle be added to a simplified model like this one, and lies within the scope of future investigations. 

Even though the BMSS scenario poses as an estimate for the upper limit of the thermalization time, I  would like to argue that this paper presents theoretical evidence that pre-equilibrium photons may be  essential to fully understand the photon spectra recovered from heavy ion collisions. These results, as well as other previous works <cit.>, seem to suggest that the  direct photon puzzle may be solvable by the thorough and 
spatially resolved calculation of non-thermally sourced photons.



§ ACKNOWLEDGEMENTS

I thank Jürgen Berges, Eduardo Grossi, Aleksas Mazeliauskas, Jorge Lopez, David Lafferty, Nicole Martin, Klaus Reygers, Naoto Tanji for fruitful discussions. I would like to thank Prithwish Tribedy for giving me the IP-Glasma saturation scale data, as well as Werner Vogelsang for sharing his pQCD calculation. 
This work is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)". 


§ ANALYTICAL DERIVATION OF 2↔ 2 YIELDS


This Appendix works as a brief guide to the calculation of the LL Glasma yields. I will only include here the derivation of stages (i) and (ii) of the BMSS scenario, since stage (iii) was briefly explained in sec. <ref>.

 §.§ Stage I

For the calculation of the invariant yield of the first stage one starts with the space-time integration of the rate in eq. (<ref>),

    1/S_⊥dN_γ/d^2 p_⊥dy =    Q_s^2 ∫^τ_1_τ_0dτ τ∫^∞_-∞ dη
         ×1/Q_s τ f_q(τ,p_⊥, p_z).

Here, the transverse dependence is taken to be homogenous, and parametrised by S_⊥, which can be taken then from the Glauber model.  By taking p_z/p_⊥=  sinh(y-η) ≡ v, and defining the dimensionless time τ̃≡ Q_s τ one can find 

    1/S_⊥dN_γ/d^2 p_⊥dy =     f_0 Q_s/p_⊥ ℐ_2/3(p_⊥)  ,

where the unitless function ℐ_a(p_⊥) is given by 

    ℐ_a(p_⊥) =∫^τ̃_1_τ̃_0dτ̃/τ̃^-a∫^∞_-∞ d v/√(1+v^2 ) e^-1/2(p_⊥ v/σ(τ̃))^2],.

The v integral can be performed to find the expression

    ∫^∞_-∞ d v/√(1+v^2 ) e^-1/2(p_⊥ v/σ(τ̃))^2 =e^1/4(p_⊥/σ(τ̃))^2 K_0(p^2/4 σ^2(τ̃)) ,

where the ℐ function can be fully integrated to find the invariant yield,

    1/S_⊥dN_γ/d^2 p_⊥dy =   10/3π^4α_e κ_g/√(2π) f_0 σ_0/Q_s  Q_s^2/p^2_⊥ 
         ×.  G ^2,2_2,3( p^2_⊥/2 σ^2(τ)|[      1 ,  1; 1/2, 1/2, 0 ]. ) |_τ_0^τ_1 .


Here, G stands for the Meijer-G function (see Ref. <cit.>). By keeping p_⊥ fixed and expanding up to leading order in σ_0/p_⊥ and substituting Q_s τ_0=1 and  Q_s τ_1=α_s^-3/2, one can obtain the simplified and more meaningful expression for the photon p_⊥ resolved multiplicity. This gives equation (<ref>), namely


    1/S_⊥dN_γ/d^2 p_⊥dy =  √(π/2)20/9π^4α_e κ_g     f_0 σ_0/Q_s log(α_s^-3/2) 
       × Q^2_s/p^2_⊥  W_r[p_⊥-Q_s]  ,
     .


As it  was stated in the main text, this expansion can be trusted for values corresponding to p_⊥ > σ_0. Below that, the approximation fails and, in fact, gives a diverging total yield. To find the  yield per unit rapidity, dN_γ /dy, the full result, eq. (<ref>), has to be integrated. This gives the function 


    .1/S_⊥dN_γ/dy|_p_⊥<Q_s =   10/3π^3α_e κ_g/√(2 π)  f_0 σ_0  Q_s 
       ×.G_3,4^2,3(. p^2_⊥/2 σ^2(τ_1)|
    [       1,1,1; 1/2,1/2,0,0;             ])|_τ_0^τ_1 .



Once again, we can expand this expression to find (<ref>) at LO in terms of the normalized anisotropy parameter, σ_0/Q_s, which once again, we state, confirms the results of the former parametric estimate <cit.>.


 §.§ Glasma II

For the second stage of the Glasma, the integration proceeds in the same fashion. We start now with the rate at stage II,  


    1/S_⊥dN_γ/d^2 p_⊥dy =    Q_s^2 ∫^τ_2_τ_1dτ τ∫^∞_-∞ dη
         ×1/(Q_s τ)^1/2 f_q(τ,p_⊥, p_z) .

By making the same substitutions as before, one can find the general expression 

    1/S_⊥dN_γ/d^2 p_⊥dy =     f_0 Q_s/p_⊥ ℐ_1/2(p_⊥)  ,

where integrating ℐ_1/2(p_⊥) as above, one can find the total result for the invariant yield of the second stage of BMSS scenario,

    1/S_⊥dN_γ/d^2 p_⊥dy =   10/3π^4α_e κ_g/2^5/4√(π)  f_0  σ_0/Q_s  W_r[p_⊥-Q_s]  Q_s^2/p^2_⊥ 
         ×.  G_2,3^2,2(. p_⊥^2/2 σ^2(τ)|
    [     1,7/4; 5/4,5/4,0;           ]) |_τ_1^τ_2 .


To find a simplified version of the rate for the kinematic window of interest, we expand in terms of σ_0/p_⊥ and substitute τ̃_1 = α_s^3/2 and τ̃_2= α_s^5/2  find eq. (<ref>), namely

    1/S_⊥dN_γ/d^2 p_⊥dy =    √(2 π)20/9π^4α_e κ_g f_0 σ_0/Q_s Q^2_s/p^2_⊥ 
       × W_r[p_⊥-Q_s]  (α_s^-1/2-1) .


Just as with stage (i), this yield is divergent at low-p_⊥ and to be able to find the total yield, one has to integrate eq. (<ref>). I find the expression


    .1/S_⊥ Q^2_SdN_γ/dy|_p_⊥<Q_s =   20/3π^3α_e κ_g/2^3/4√(π) f_0 (σ_0/ Q_s)^5/2
       ×.G_3,4^2,3(. 1 /2 σ^2(τ)|
    [     1,7/4,7/4; 5/2,5/2,0,3/4;               ])|_τ_1^τ_2.


From which, once again, I can confirm the parametrical results from <cit.> by expanding in terms of  σ_0/Q_s, as in eq. (<ref>)


§ THERMAL PHOTONS



For photon production from a thermal medium, I may further simplify the computations 
by using thermal Boltzmann distributions for quarks and gluons at high energies in (<ref>). This leads to 

    E dN^th/d^4X d^3p = C 5/9/2π^2T^2 e^-E/T  .


From this expression the photon multiplicity is obtained by integrating (<ref>) over the four-volume of the evolution,

    dN^th/dyd^2 p_⊥  = C 5/9/2π^2∫τ η ^2 x_⊥ τ  T^2  e^-p_⊥cosh(η-y)/T .

Here, I made the equation explicitly Lorentz invariant by 
E → p^μ u_μ = p_⊥cosh(η-y) 
with the comoving four-velocity u_μ = (coshη, 0,0 sinhη). 




The integration itself depends on the spacetime temperature profile. This can be derived from simplified hydrodynamical models or from direct simulation. The simplest scenario would be the evolution for a system invariant under η,  and ϕ transformations. From basic symmetry arguments, one can get the evolution of the energy density, 

    ϵ = ϵ_th[τ_th/τ]^4/3 .

Using that ϵ∼ T^4, one gets the τ-dependent temperature profile of the expansion, 

    T(τ) = T_th[τ_th/τ]^1/3  ,

where T_th is the thermalization temperature, and τ_th is the proper time at which it thermalizes. In the bottom-up thermalization scenario <cit.>, their parameteric dependence is given by

    T_th∼ c_T c_eq α_s^2/5  Q_s  ,
    τ_th∼ c_eq α_s^-13/5  Q_s^-1 .

Here, c_eq is a coefficient of order unity, which arises from the uncertainty from the parametric dependence of the thermalization time, τ_th. The other coefficient, c_T,  is a constant needed to finish to constraint the thermalisation temperature in the BMSS scheme <cit.>. These coefficients are constrained in sec. <ref> using the method in ref. <cit.>.
For the thermal epoch, the system will evolve from the thermalization time, until the critical time, τ_c, which signals the arrival to the critical temperature. At this point, the deconfined quark-gluon-plasma phase transitions via a crossover to the hadronic phase. For this work I will take T_c	= 0.154  GeV <cit.>, while τ_c can be from the temperature profile as follows


    τ_c = τ_th(T_th/T_c) .


To start, I want to derive the photon multiplicity for the thermal case. For this, as it was noted before, I integrate emission function, eq.<ref> ,  over the 4-volume. This gives


    dN_γ/d^2 p_⊥dy = S_⊥ C̃ 5/9  /2π^2∫_τ_th^τ_c   τ τ∫_-∞^∞η τ  T^2(τ) 
        × e^-cosh[η-y]/T(τ) .


Here, the temperature profile is assumed to be homogeneous in the transverse plane, which yields the the x_⊥ integration trivial.  The integration can be continued by the rapidity dependent part of the integrand. Using 

    2  K_0(z) = ∫_-∞^∞η  e^ -z coshη ,

where K_n(x) stands for the modified Bessel function of order n. After this integration, I find 


    1/S_⊥dN_γ/d^2 p_⊥dy  =   5/9  /π^2∫_τ_th^τ_cτ τ  T^2(τ)K_0(/T(τ))   .


By transforming the integration from proper time to inverse temperature, β=1/T, one can get the integral 

    1/S_⊥dN_γ/d^2 p_⊥dy    =   5/3  /π^2 τ^2_th  T^6_th∫_β_th^β_cβ β^3  K_0(β )
        ≡ 5/3  /π^2 τ^2_th  T^6_th B_th() .


The function B_th()  can be found by analytical integration, to  get 

    B_th() =   {β_th^4 [K_1(β_th)/β_th+2 K_2(β_th)/(β_th)^2] 
       -β_c^4 [K_1(β_c)/β_c+2 K_2(β_c)/(β_c)^2] },

It is important to note that this function is reliable when >T, since the emission function I used had the assumption of E>>T. Finally, one can find the total yield by integrating over the transverse momentum 

    1/S_⊥dN_γ/dy    =   10/3  /π τ^2_th  T^6_th∫_0^∞d p_⊥  p_⊥ B_th())
        ≡ 5/3  /π τ^2_th  T^4_th ( T^2_th/T_c^2 -1 ) .

where I find the same yield than in ref. <cit.> 
Here, as in the previous estimates, I have integrated from p_⊥ =0 instead of p_⊥ = T, where the formula is valid. Here I quantify the  error to be of order K_2(1), which will give a relative error of 𝒪(1). 

apsrev4-1


