Probing the evolution of heavy-ion collisions using direct photon interferometry
Oscar Garcia-Montero Nicole Löher Aleksas Mazeliauskas Jürgen Berges Klaus Reygers


APS/123-QED









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

^bPhysikalisches Institut, Universität Heidelberg, Im Neuenheimer Feld 226, 69120 Heidelberg, Germany




We investigate the
measurement of Hanbury Brown-Twiss (HBT) photon correlations as an experimental tool to
discriminate different sources of photon enhancement, which are proposed to
simultaneously reproduce the direct photon yield and the azimuthal
anisotropy measured in nuclear collisions at RHIC and the LHC. To showcase
this, we consider two different scenarios in which we enhance the yields from standard hydrodynamical simulations.
In the first, additional photons are produced from the early pre-equilibrium
stage computed from the “bottom-up" thermalization scenario. In the second, the thermal rates are enhanced close to the
pseudo-critical temperature T_c≈ 155 MeV using a phenomenological ansatz.
 We compute the correlators for relative momenta
, and  for different  transverse pair momenta, , and find
that the longitudinal correlation is the most sensitive to different photon
sources. Our results also demonstrate that including  anisotropic pre-equilibrium
rates enhances non-Gaussianities in the correlators, which can be quantified
using the kurtosis of the correlators. Finally, we study the feasibility of
measuring a direct photon HBT signal in the upcoming high-luminosity LHC runs. 
Considering only statistical uncertainties, we find that with the projected ∼ 10^10 heavy ion events a
measurement of the
HBT correlations for <1 GeV is statistically significant.







Valid PACS appear here

Probing the evolution of heavy-ion collisions using direct photon interferometry
    Klaus Reygers^b
    December 30, 2023
================================================================================








§ INTRODUCTION


The relativistic nuclear collision experiments explore the physics of 
dense and hot QCD matter, also known as the quark-gluon plasma (QGP) <cit.>. The bulk properties
of this new state of matter are inferred indirectly from the yields and correlations
of the produced hadrons. However the QCD degrees of freedom participate in the 
strong interaction and are subject to the effects of multiple-rescatterings and
non-perturbative physics of hadronization, which tend to erase the information about
the earlier stages of the collision.
Electromagnetic probes, e.g. photons and dilepton production, are therefore often championed as
penetrating probes of the QGP dynamics <cit.>.
Although it is true that photons escape virtually unscathed from the medium,
the continuous electromagnetic emission makes it very hard to discriminate
between different photon sources. 
Furthermore, in the standard hydrodynamical picture, it is  challenging to
simultaneously describe the measured photon yields and their azimuthal
anisotropy, which is commonly referred to as the  direct-photon puzzle
<cit.>.

In this paper, we explore two-photon interferometry, called in this
context femtoscopy, as a tool to untangle the space-time evolution 
of the QGP and in order to shed light on the direct-photon puzzle.
This addresses the question whether direct photons in heavy-ion collision (HIC)
originate predominantly from the early or the late stage of the collision.
Specifically, we use Hanbury Brown-Twiss (HBT) correlations,
which are the only known way how to directly extract space-time information
from the particles measured in heavy-ion collision experiments <cit.>. 
HBT correlations, originally introduced to measure the radii of stars
from the incoming photons <cit.>, have been used
extensively across physics, from atomic gas correlations in cold atom experiments <cit.>, to
pion interferometry in 
heavy ion collisions
experiments <cit.>.
Interferometry of direct photons as a tool to study the space-time
evolution of a heavy–ion collision was theoretically explored by several
authors, see
<cit.>
and references therein. So far only one measurement in Pb–Pb collisions at
√(s_NN) = 17.3 GeV at the CERN SPS was reported
<cit.>. In view of the upcoming high-luminosity runs at the
LHC <cit.>, we expect further photon measurements at the TeV energy scale and
therefore present theoretical and experimental analysis of  the HBT signal.

In this work we study HBT correlators in different scenarios. First,
we compute the yield and HBT correlators for a hydrodynamically expanding
quark-gluon plasma  and the subsequent hadronic stage using realistic 2+1D event-by-event
simulations of a heavy ion collision.
We then consider two additional sources of
photons, coming from early and late stages of the expansion respectively. At early times 
 we supplement the thermal yield by including a pre-equilibrium
contribution, which was found in previous work to be on par with the
thermal one <cit.>. Motivated by the idea that thermal rates might be
enhanced  around the pseudo-critical temperature by confining modes during hadronization <cit.>, we 
add another source for photons at late times.
We present a detailed analysis of photon HBT signal
sensitivity to different photon sources and make a realistic estimate of
experimental
statistics needed to measure these signals by the ALICE detector.



§  HBT CORRELATIONS


Quantum statistical effects can be used to understand the spacetime distribution of particle sources <cit.>. In the context of HIC experiments, we are interested in finding the spatial extension of the photon source in the fireball. For this, we use the normalized HBT correlator, 

    C(1,2) = 12/12  ,

where the numerator is given by the two-photon distribution, which can be expressed in terms of asymptotic states, i.e. creation and annihilation operators of a gauge field

    N/^3 p_1 ^3 p_2  =  ∑_λ_1,λ_2⟨ a^†_1,λ_1 a^†_2,λ_2 
    a^ _2,λ_2 a^ _1,λ_1⟩ .

Here, n and λ_n are, respectively, the spatial momenta of the detected photons and polarization mode of the nth photon. In a field theoretical language, this can be computed generally from a four-point correlator of gauge fields in momentum space, for equal in- and outgoing momenta.
The denominator is the product of the invariant yields, and can be expressed with asymptotic states as follows,

    N/^3 p =  ∑_λ⟨ a^†_,λ a_,λ⟩ .


We can further simplify Eq. (<ref>) by splitting the four-point function into connected and disconnected parts. The photon fields during a HIC are not expected to   be highly occupied in-medium. This renders the electromagnetic sector to be a dilute gas of particles, for which the photon-photon interaction vertex is very small. In this case we can loose the connected part, and Wick's theorem states that

    ⟨ a^†_1,λ_1 a^†_2,λ_2 
    a^ _2,λ_2   a^ _1,λ_1⟩
    ≃    ⟨ a^†_1,λ_1 a_1,λ_1⟩⟨ a^†_2,λ_2 a_2,λ_2⟩
       + ⟨ a^†_1,λ_1 a_2,λ_2⟩⟨ a^†_2,λ_2 a_1,λ_1⟩ .
	

From this it can be seen that the two-photon correlator splits into a trivial (diagonal)  and non-trivial (off diagonal) part. It was shown in Ref. <cit.> that these correlators can be directly related to scalar Wigner density functions S(x,K) (also called emission function in the literature), where the information over polarization of the sources can simply averaged out using the Ward-Takahashi identity. The correlator is found to be

    C(,) = 1 + 1/2| S(,) |^2/S(0, 1) S(0, 2)  ,

where S(,) is the Fourier transform of the emission function,

    S(,) = ∫^4 x  e^ x· q  S(x,) .

The result is a version the scalar HBT correlator, modified only by a relative degeneracy factor of 1/2. The new variables, q =p_1 - p_2 and K=(p_1 + p_2)/2 are the relative and average momenta for two photons, respectively. In what follows, both S(,) and the photon invariant yield, S(0,),  will be calculated by associating the Wigner function with photon emission rates, that is 


    S(x,)↔ E_K  N /^4 x  ^3 K .



 §.§ Variables and Approximations

The detected photons are on-shell, and we express the photons four-momenta

    p^μ = (p_⊥cosh y, p_⊥cosφ, p_⊥sinφ, p_⊥sinh y)

with rapidity y, transverse momentum p_⊥  and azimuthal angle φ.
For the average and relative momentum variables, q and K, defined above, we choose a coordinate system such that 

    K^μ   = (K^0, , 0, K^z )
    
    q^μ   = (q^0, , , ) ,

i.e.  lies in the x-z plane, with z being the beam direction. The q components are called the longitudinal, outwards, and side momenta. We can express them using 

    = (·)/
       = | - (·)/|.


Since both photons are on-shell, both the  pair and relative momenta will be off-shell, and for two identical particles, they satisfy 

    q_μ K^μ = 0    ⇒   q^0 = ·/K^0.


In the case of pion-pion interferometry, two approximations are taken to
further simplify the computation of the HBT correlator. In the literature they
are commonly referred as the on-shell and smoothness
approximation <cit.>. For the former, the pair momenta itself is taken to be
on-shell, K^0 ≈|K| +𝒪(^2). In hadron
interferometry, this can be used to good approximation because of the large
masses of hadrons. Even for pions, the subleading terms are suppressed by E_K
for all average momenta. In general, in such calculations, if the corrections
are not suppressed by powers of ^2/^2, they are controlled by the group
velocity β =/K^0 <cit.>.
However, photons are massless, and this expansion will break at
||/(2||)∼ 1. Unfortunately, in experimentally realizable settings at
the LHC, the direct photon signal is  contaminated by photons from decays,
which form the vast majority of the signal. This leaves direct photons with a
deficiency in statistics. As a consequence, photon pairs cannot be correlated
for infinitesimal ||, with reasonable confidence (see Sec. <ref>). This means that in general || and || will be
on the same order of magnitude. Nonetheless, for a single Gaussian source, the
correlator half-widths can be computed using this approximation without any
problem. For such a source, the correlator will be perfectly Gaussian and the
radii can be directly extracted by fitting the curves <cit.>, or
by computing the curvature of the correlator at =0. In the case of direct
photons, we will be having photons from different sources (stages of the
fireball) which will present different scales. Thus, the condition
||/(2||)≪ 1 cannot be met for all the kinematic regime. Furthermore,
the Wigner function in Eq. (<ref>) is generally given for any
combination of momenta. The function can be evaluated off-shell
<cit.>, and therefore to avoid unexpected
deviations coming from uncontrolled terms we choose to not use it.

The other approximation normally used in the literature is called  the smoothness approximation, which consists of  neglecting the  dependence in the denominator of Eq. (<ref>), via S(0, 1,2) → S(0, ). The correlator is given in this limit 
 
    C(,) = 1 + 1/2|S̃(,)/S̃(0, ) |^2 ≡ 1+1/2⟨ e^i q· x⟩ ,

for which we will introduce the commonly used averaging notation <cit.>: 

    ⟨ f(x) ⟩ = ∫^4 x f(x) S(x,K) /∫^4 x S(x,K)  .


The smoothness approximation is accurate if the curvature logarithm of the one particle distribution is small <cit.>, which is not true for photons at small momentum (see Fig. <ref>). In reference <cit.> it was found that the convergence of the smoothness approximated to the full correlator is restricted for values of ||< 2 ||. This is the same scale that signals the breakdown of the on-shell approximation. We use this approximation for the case of 1D slices for  and , where the other directions of q are set to zero. In this case, the  and  direction look Gaussian, and the explored  values will be larger than the inverse half width of the correlator, which makes this approximation safe. 



 §.§ Homogeneity Radii


We can get a general form of C(,) for an arbitrary Gaussian source around the origin in  space

    C(,) =1 + 1/2exp[-q_μ R̃^μν  q_ν]  ,
 
For sources with relatively small non-Gaussianities this approximation is still valid, since the perturbations around C increase only at high  values <cit.>. 
The half-width tensor, R̃^μν = R̃^μν(), is a function of the pair momentum. To clean the notation, we avoid writing its  dependence. Using the orthogonality relationship, Eq. (<ref>), we can reduce this expression to 

    C(,) =1 + 1/2exp[-q_i  R^ij q_j ] ,

 by redefining R^ij≡β^i β^j R^00 +2 β^i R^0j +R^ij.  Because of symmetry, R^ij =R^ji, we only get 6 independent components. Using the relative momentum parametrization introduced above, we can express it as 

    R_ij()=
    [ ^2 ^2 ^2; ^2 ^2 ^2; ^2 ^2 ^2 ] .


In this work we only focus on the diagonal of this matrix. While it has been shown that for longitudinally expanding sources the  term is relevant <cit.>, it is also true that numerically calculating such cross-terms is more computationally complex.  

To compute the radii in Eq. (<ref>), we use the the method of moments, which is stable for correlators with strong non-Gaussianities <cit.>. We use the moments of the true correlator C(, )-1 in relative momentum space, 

    ⟨⟨ q_i q_j ⟩⟩  = ∫^3q  q_i  q_j   g(;)≡1/2(R^-1)_ij ,

where R^-1 is the inverse matrix of Eq. (<ref>).  We have defined the  distribution function 

    g(;) ≡C(,)-1/∫^3 q [C(,) -1]

to ensure correct normalization. Because of the symmetry properties of the correlator, we can safely assume the one-point functions vanish, ⟨⟨ q_i ⟩⟩ = 0. For simplicity, and because we do not explore the off-diagonals, we will keep the notation one-dimensional.  That means that the homogeneity radii are going to be given by 

    R_i^2 = 1/2⟨⟨ q_i^2 ⟩⟩^-1 ,   with   i∈{l,o,s}.


It is important to clarify that this method requires the correlator to be highly localized around =0, to give sensible results for the characteristic scale. In other words, the correlator needs to decay faster than a power-law. We can use also this method to quantify the deviations from Gaussianity by computing the normalized  excess kurtosis, 

    Δ_i = ⟨⟨ q^4_i ⟩⟩/3⟨⟨ q^2_i⟩⟩^2 -1,

which, as expected, vanishes in the Gaussian limit. In theoretical calculations of HBT correlations, going to higher values of q_i requires only better numerical precision. However, it may be problematic for experiment, where high relative momentum values will suffer from statistic limitations.




§  MODELLING THE PHOTON SOURCES

As it was stated in the introduction, we calculate the thermal photon observables, which are enhanced by the inclusion of early- and late-time photon sources. The thermal base is calculated from hydrodynamic simulation using the VISHNU package <cit.>, from which realistic space-time evolution of temperature and velocity fields was obtained. Using the default model parameters tuned to the experimental data, we simulated 200 Pb-Pb collision events  at the centre of mass energy
 √(s_NN)=2.76 in 0-20% centrality class. The initial conditions at τ_hydro=0.6 were provided by the two-component Monte Carlo Glauber model <cit.>. The relativistic hydrodynamic simulation was then performed using fixed shear viscosity over entropy ratio η/s=0.08 and the decoupling energy density e=0.1/fm^3.
 The space time evolution of transverse velocities v_x and v_y and temperature T was recorded on
 a coarsened grid with spacing dx=dy=0.4 and dτ=0.2 (x_max=y_max=25.2). The final time τ_max varied depending on the initial conditions, but at least 100 recorded events had τ_max≥ 15.8. We calculate photon emission for each event separately and then do the ensemble average.
 
Direct photons can be emitted from the QGP and hadron resonance gas (HRG) epochs of the evolution of the fireball. The transition from the QGP production to the HRG is signaled by a switch at 160 MeV. It is assumed that the emission threshold for thermal photons is at a temperature of 120 MeV. In addition, two possible sources for enhancing the invariant photon yield are discussed. The first is the inclusion of a pre-equilibrium source based on the first stage of the bottom-up thermalization scenario <cit.>. The second source is a phenomenological enhancement of the thermal rates near a pseudo-critical temperature T_pc, presented first in Ref. <cit.>. We discuss these  and other photon contributions below. 
 




 §.§ Prompt photons from the initial stage


During the initial stage of the collision, prompt photons are produced via hard scattering of the partons from the individual nucleons. The photon cross-section for the NN →γ X process can be calculated using perturbative QCD (pQCD) <cit.>, which is then scaled by the number of binary collisions,  N_coll, via the relation

 
    d N_prompt/d^2 p_⊥dy= N_coll/σ^NN_inel dσ^NN →γ X/d^2 p_⊥dy .

Here σ^NN_inel is the total inelastic collision for a collision of two nucleons. We compute N_coll using the optical Glauber model. For the computation of the full photon invariant yield we need to extend the pQCD computation to smaller p_⊥ values. We do so by taking the same parametrization used by PHENIX Ref. <cit.>. The fit function is given by the functional form 

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

Because this contribution takes on account incoherent production of single photons, we do not include prompt photons in the calculation of the HBT correlator, but add them to the total photon yield.




 §.§ Photon emission from the quark-gluon plasma




To compute the photon contribution due to the thermal QGP we use the full leading order (LO) computation, parametrized in Ref. <cit.>. This rate contains not only the two-to-two contributions which dominate at higher momenta, but also near-collinear bremsstrahlung and the inelastic pair annihilation, thereby fully including the Landau-Pomeranchiuk-Migdal effect (LPM), which can be understood as suppression of emission owing to interference of multiple scatterings <cit.>. The parametrization used in this work is given explicitely in Appendix <ref>.




 §.§ Photon emission from the hadron resonance gas


For the thermal photon emission rate from the hadron resonance gas phase the parametrizations of Ref. <cit.> is used. The given parametrizations agree within 20% with the microscopic calculated values. Microscopic calculations have already been performed <cit.>, but, as pointed out in <cit.>, the results cannot be easily used in models like the one described here. Two different parametrizations for the photon emission rate are given: one for the contribution from the in-medium ρ mesons and one for the contribution from bremsstrahlung originating from ππ scattering. They can be applied to photons with energies q_0 between 0.2 and 5 GeV, which are produced from chemically equilibrated matter with a temperature between 100 and 180 MeV and baryon chemical potentials of 0 to 400 MeV. In the case of ALICE, vanishing chemical potential is assumed.





 §.§  Photon production from pre-equilibrium


Using the  “bottom-up" thermalization scenario <cit.>, recent estimates <cit.> show the pre-equilibrium contribution photons to be non-negligible. The central idea in this estimate is that gluon saturation takes place at RHIC and LHC energies, which means that during the initial stage of a collision, the nuclei behave as macroscopic fields, and undergo very strong, non-linear interactions. After a parametrically short time given by the saturation scale, τ_0∼ Q_s^-1, the gluon fields get highly occupied and undergo three stages of relaxation. During the first stage, hard modes with p_⊥∼ Q_s completely dominate the system. These modes are approximately conserved, yet diluted thanks to Bjorken expansion. During this stage, hard modes scatter via two-to-two scatterings, which produces a broadening of the distribution in the p_z direction. 
The second stage starts once the occupation of the gluon modes falls below unity, where the typical longitudinal momentum of hard gluons saturates at a finite value. In this stage, hard gluons still dominate the total gluon number, while the typical interactions are taken over by the soft sector. Finally, we arrive to the third stage of the BMSS scenario, where the number of soft and hard gluons becomes comparable. Soft gluons thermalize rapidly via two-to-two scatterings, which creates a bath to which hard gluons quickly loose energy to, via mini-jet quenching. The system has then fully thermalized. 

We fix the initial characteristic scale IP-Glasma model <cit.>, which combines the
geometry of the MC-Glauber model <cit.> with the IP-Sat model
<cit.>, while the BMSS scenario gives the time
dependence of the rates. We use as well experimental data to constraint the
needed parameters, the thermalization time was found in Ref.
<cit.> to be τ_th∼ 2.4 for LHC and RHIC
energies. Since the bottom-up scenario does not account for the
transverse expansion , such late thermalization poses a phenomenological
problem, as the photons will not be able to build up enough anisotropy,
creating tension with data. To avert this, we will only evolve the
pre-equilibrium stage up to the end of the first stage of the
bottom-up scenario, τ_hydro=0.6
[The photon spectra from all three stages of the bottom-up thermalization is studied in Ref. <cit.>].
From the field
theoretical point of view, in this stage, the gluon medium approaches a
non-thermal fixed point <cit.>,
where the gluon occupation is given by

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

Here, α_S is the strong coupling, and f_S is a scaling function, which can be parametrized from the results of classical statistical simulations <cit.> as follows

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

Here, W_r[p_⊥,Q_s] stands for a suppression function, inspired by the classical statistical simulations. It depends on a free suppression parameter, r, and it is given by 

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



At the end of this stage, the system is assumed to instantaneously thermalize, and we match the energy densities in the pre-equilbrium  and hydro stages at τ_hydro, 

    ϵ_early(τ_hydro,x_⊥) = ϵ_hydro(τ_hydro,x_⊥) ,

which gives also the spatial profile of the saturation scale Q_s(x_⊥). 
In the pre-equilibrium stage, most of the energy density resides in the gluonic sector. Using Eq. (<ref>) and the QGP energy density one can obtain

    Q_s(x_⊥)/T(x_⊥)= [√(2/π)37π^2(2 π)^2α_S/30(1 + √(2π) r  + 2 r^2))τ_hydro/τ_0f_0 s]^1/4

where τ_0 is taken to be the spatially averaged saturation scale, ⟨ Q_s⟩, and can be determined parametrically using the method described in Ref. <cit.>.

For the rate, we will use a kinetic rate, generally given by 

    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)] ,

where the processes included are the two-to-two annihilation, qq→ gγ, and Compton scattering, qg→ qγ. Because the computation at each space-time point of such rate requires a 5-dimensional integral, we simplify the rate using the small angle approximation.  For massless mediators, hard scatterings present collinear enhancement, which will dominate the integrals in Eq. (<ref>). Expanding in the exchange momentum of the mediator and keeping only the leading term one finds the simplified rate <cit.>, 

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

where α is the electromagnetic coupling, κ_g=c (2N_c)^-1, where N_c is the number of colors, and c is the gluon liberation factor described in Ref. <cit.>. The quark distribution, f_q is taken from hard splitting of gluons in-medium, namely  f_q∼α_S   f_g. That is, using this parametrization, we assume the quark distribution inherits the scaling properties of the parent gluons. To avoid breaking fermion statistics, we suppress the quark distribution for low p_⊥ values, so that f_q=1/2 at its highest value. The ℒ term is called the Coulomb logarithm, and it is a regulator, which relates the UV and IR scales, two cutoffs which are needed for this approximation. In the thermal case, the UV scale can be related to the temperature, T, while the IR scale can be related to the Debye mass, m_D ∼ gT. Using this identification, the leading-log (LL) thermal rate from Ref. <cit.> can be found from the small-angle approximated rate. 

Nevertheless, at the full leading-order (LO) limit of the photon rate, Ref. <cit.>, it was shown that in a thermal setting, photon rates are dominated by near-collinear bremsstrahlung for photon energies or p≲ 2 T, while at 2 T≲ p ≲ 10  T, the two-to-two terms are of the same order to the near-collinear contributions. The modification for the rate is applied then by changing the constant under the log

    ℒ→ν_LO(x)

where x =E/T in the thermal case, and ν_LO(x) is given in Eq. (<ref>). We expect a similar behavior to the pre-equilibrium stage, with one difference. During this stage, the characteristic momentum scale is taken to be the saturation scale Q_s, making the near-collinear contributions during the early stages dominant at p≲ 2 Q_s which for the center of mass energy at ALICE is most of the kinematic window at which direct photons are observed. We therefore also use the modification of Eq. (<ref>) in Eq. (<ref>), for x→ x'=E/Q_s.



 §.§  Critical enhancement at late times near T_c


To account for the missing photons one could naively push the initial time to smaller values. Nevertheless, doing so hardens the spectrum, which creates tension with the experimental results <cit.>. If one has to increase the thermal rate, it has to be done increasing the weight of photons coming from later times.This is in line with the idea suggested in Refs. <cit.>, where it is conjectured that the thermal rates are enhanced near a pseudo-critical temperature T_c, 

    E d N_enh/^4 x ^3 p≡  h(T)  E d N_thermal/^4 x ^3 p

by the fact that close to the transition to hadronic degrees of freedom, one has to account for interactions related to confinement. This means that the partonic cross-sections will see a rise which cannot be accounted for by perturbative physics <cit.>. 
For the purpose of this paper, however, we choose to model the enhancement factor, h(T), as follows

    h(T)    =1+h_0  e^-(T-T_c)^2/d^2

where the pseudo-critical temperature is set to be T_c=155. The enhancement parameters are set to be h_0= 3 and d = 50. The enhancement factor is tuned such that the enhancement matches the experimental results from the ALICE collaboration, see Fig. <ref>.



§   RESULTS




We compute the total yield for the thermal baseline, and include as well the both enhancement scenarios, which can be seen in Fig. <ref>.  The pre-equilibrium photon spectrum shows a structure around p_⊥∼ 2.5. This shoulder comes directly from the parametrization of the quark function. Nevertheless, the specific value at which we can find the peak is given by averaging the space dependence of Q_s(x_⊥). The pre-equilibrium spectrum is found to be dominant for 2 < p_⊥∼ 3, while being relatively small in the IR sector. Summing over the prompt, pre-equilibrium and thermal contributions we find good agreement with ALICE data for central collisions, 0-20% (Fig. <ref>, left).
On the other hand, applying the enhancement to the thermal rates, Eqs. (<ref>) and  (<ref>),  just as expected, we see an overall increase of the spectrum, particularly strong for low-p_⊥, photons. It can be seen that both scenarios are compatible with the errorbars, which means that distinguishing such cases experimentally is not possible using only the invariant yield. 



The full HBT correlator, Eq. (<ref>), was computed for midrapidity pairs, K_z=0, along the three diagonals, i.e. q_i with q_j=q_k=0 for i≠ j≠ k. We focus on 0-20% central collisions in ALICE, with √(s_NN)=2.76, where the average saturation scale is ⟨ Q_s^2⟩ = 2.9^2. 
As expected, the longitudinal curves are the most sensitive to the inclusion of both enhancements which are presented in Fig. <ref> for different values of . Although the correlator around the side- and outward diagonals show a difference with the inclusion of both enhancements, the effect is noticeably small. This can be seen better for the diagonal radii,  ,  and  (see Fig. <ref>), which were computed using the characteristic scale method and the aforementioned correlators. 



 Just as expected from the correlators, the change in the longitudinal radius, , is the largest one. The change induced in  by the inclusion of the scenarios was found it to be in the 10-20% range for the outward direction, and 0-10% for the sideward direction. The small change in the transverse radii will make using them to discriminate models difficult. Nevertheless, this gives an interesting case for predictions. Take, for example, the pre-equilibrium case: If pre-equilibrium photons are relevant at the yield level, and the assumption that the pre-equilibrium stage does not create enough pressure gradients is correct, thermal models will be able to reproduce the  but may undershoot significantly . On the other hand, a consistent increase with K_⊥ on the three radii may indicate that photons come from the late stages. 



We  also computed the normalized excess kurtosis, Eq. (<ref>),  for the three diagonals. A clear hierarchy is found, where  breaks Gaussianity the most, followed by  and . We find that the sidewards direction is to good approximation Gaussian (see Fig. <ref>). The non-Gaussianities, as was explained above and in Ref. <cit.> arise from the longitudinal expansion of the fireball. In the case of massless particles these effects will be considerable more important than for e.g. pions. Additionally volume emission will further enhance these effects, opposed to Cooper-Fry surface emission. Non-Gaussianities are quite intuitive to understand in the case of the  direction, since the boosting from longitudinal expansion is largest for the  variable. However, the easiest way to see how the outward direction gets contributions from the expansion is the definition R^ij≡β^i β^j R^00 +2 β^i R^0j +R^ij. From this  formula we see that for the outward direction,  gets a non zero contribution from β_o t =t  K_⊥/K^0, while the sideward direction, by the definition, will not. This means that the outward homogeneity radius not only depends on the spatial size of the source, but also on the lifetime of emissions <cit.>. As it can be seen in Fig.<ref>, the normalized excess kurtosis can be used as an observable complementary to the radii. This is particularly true for K_⊥ < 0.5, where the big difference in Δ_l could be used to differentiate the scenarios.








§  EXPERIMENTAL FEASIBILITY


Measuring direct-photon Hanbury Brown-Twiss correlation is a challenging task. At the LHC, the ALICE experiment measures photons at low transverse momentum (≲ 3GeV) <cit.>. Significantly improved data-taking rates in the upcoming LHC runs 3 and 4 make it possible to collect a sample of Pb–Pb collisions corresponding to an integrated luminosity of 10 nb^-1, or 𝒪 (10^10) collision events. In this section we estimate up to what photon pair transverse momentum  a direct-photon HBT measurement might be possible. 

We concentrate on the longitudinal momentum difference . Statistical uncertainties for measurements of  and  are very similar. For a Gaussian parameterization the correlation function C of direct photons for =  = 0 is given by

    C = 1 + λexp(- ^2 ·^2)
 
with λ = 1/2. The total number of photons, however, is dominated by photons from neutral pion and eta meson decays. Owing to the long lifetime of the neutral pion and the eta mesons the decay photons are not correlated with the direct photons and dilute the measured correlation function, resulting in

    λ = 1/2 r_γ^2,    r_γ = N_dir/N_inc
 
for the correlation strength of pairs of inclusive photons. Here N_dir denotes the number of direct photons and N_inc the number of inclusive photons, i.e., the sum N_inc = N_dir + N_dec of the number of direct and decay photons. We assume a p_⊥-independent fraction of direct photons of N_dir/N_inc≈ 0.1 corresponding to λ = 0.005 <cit.>.

The basis for our estimate is the direct-photon spectrum in 0–20% Pb–Pb collisions at √(s_NN)=2.76 TeV measured by ALICE <cit.>. We parameterize the spectrum by 

    . 1/2 π p_⊥ N_evtd^2 N_dir/d p_⊥dy|_y=0 
    = A exp(- p_T/T)

where the inverse slope parameter is set to T=0.3GeV, see Fig. <ref>.


From this simple parametrization of the measured direct-photon spectrum we calculate the number N_p,u^dir of uncorrelated pairs of direct photons per event in a given  bin. We consider a measurement of C() in 10 MeV wide bins for || < 30MeV and || < 30MeV in various intervals of pair transverse momentum . 

The statistical uncertainty of the total number C · N_evt· N_p,u^inc of pairs of inclusive photons should be much smaller than the number of pairs (C-1) · N_evt· N_p,u^inc above the uncorrelated background. Here N_evt denotes the number of considered Pb–Pb collisions. This corresponds to

    √(C · N_evt· N_p,u^inc)≪ 
    N_evt·(C-1) · N_p,u^inc.

Neglecting the small √(C) term on the left hand-side, the criterion for a significant measurement in the considered bin reads

    σ_rel^inc≪ C - 1 = 1/2 r_γ^2

where 

    σ_rel^inc = 
    1/√(N_evt N_p,u^inc) =
    r_γ/√(N_evt N_p,u^dir).

Results for the statistical uncertainty σ_rel^inc of the measured correlation C() for inclusive photons for N_evt = 10^10 are given in Table <ref>. This table also shows the ratio s = 2 σ_rel^inc / r_γ^2. A value s ≪ 1 corresponds to a significant measurement. We consider the case of a full photon detection efficiency (1) and the case of a limited detection efficiency (2).

From Table <ref> one can conclude that with N_evt = 10^10 Pb–Pb collisions there is enough statistics to measure direct-photon HBT correlations up to a pair transverse momentum of ≈ 1 GeV. For this value of  we illustrate the projected statistical uncertainties of C measured for pairs of inclusive photons in black in Fig. <ref>. For comparison the distribution is also shown in red for ≈ 0.5 GeV, which has much smaller projected statistical uncertainties. This provides a motivation to experimentally explore photon HBT correlation in the upcoming high-luminosity LHC runs <cit.> and to study in detail all sources of systematic uncertainties which might affect the measurement.






§ SUMMARY AND CONCLUSION


In this work, we present a case study of photon interferometry 
 exploring the space-time evolution of the fireball to
investigate possible new photon sources.
In addition to standard thermal and prompt
photons, 
we consider two different scenarios, one in which
 additional photons
are produced from the early pre-equilibrium stage, and one in which
the thermal rates are enhanced close to the transition.
In both cases the mid-rapidity direct photon yields agree with ALICE results in
central (0-20%) Pb-Pb collisions at √(s_NN) = 2.76 TeV.

We then
compute the HBT correlators in the diagonal directions, 
, and  for different  transverse pair momenta. 
In general, including photon emission from the pre-equilibrium stage 
widens the correlation because of a more compact emission source at early times.
Conversely, the late-time enhancement makes the two-photon correlation narrower.
From these correlators we extract the HBT radii    , , and  .
The longitudinal radius exhibits the largest 
difference between the thermal and the other two scenarios, namely ∼ 80% and  ∼ 20% for early and late time enhancements.
In comparison, the  and  radii
are only mildly affected, with ∼ 20% and ∼ 5% changes respectively. 


Direct photons see the entire space-time evolution of the expanding
fireball, which results in pronounced 
non-Gaussianities in the photon
 HBT signal.
 To quantify these effects, we compute the normalized excess kurtosis, which we find to be 
 largest for the longitudinal direction and significantly smaller in 
 the outward and sideward directions. At small transverse momentum, the significant
 differences in the observed longitudinal non-Gaussianities 
  provide a striking new signature sensitive to
 the different photon emission sources.



In view of the potential of two-photon correlation measurements, we perform an experimental
feasibility study.
With the projected count of
∼ 10^10 heavy-ion events at the upcoming LHC Runs 3 and 4, we determine the statistical uncertainties
of the experimental signal.
Owing 
to the photons from neutral meson decays,
the HBT signal is attenuated to a percent level. For 
transverse momenta ≲0.5 GeV
statistics will be sufficient for the measurement
of the correlation function. 
However, the differences between the early and late
time scenarios are most pronounced at higher photon-pair momenta, where statistical uncertainties
are large. Therefore, it is unlikely that
the photon interferometry alone can be used to identify the correct photon emission scenario. 
Nevertheless, we show that photon HBT signal is an experimentally accessible observable with
sensitivity to photon production physics. In conjunction with other
observables, e.g. elliptic flow, HBT correlations could be used to rule out certain
models and, therefore, 
motivate further theoretical studies and experimental estimates of systematic errors.



§ ACKNOWLEDGEMENTS

Authors thank Eduardo Grossi, Jean-François Paquet, and Johanna Stachel for valuable discussions. 
We thank Prithwish Tribedy for sharing IP-Glasma simulation results and Werner Vogelsang for sharing
the pp prompt photon data.
OGM is funded by HGS-HIRe. This work is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)". 



§ THERMAL RATES 
 

After the thermalization of the colored medium, photons can be emitted from either a thermalized quark-gluon-plasma or can be produced by hadronic processes in the hadron resonance gas phase. In the following we will summarize the rates used in this work to compute the radiation from the thermal phases. 



 §.§ Thermal rate for Quark Gluon Plasma



As indicated above, to emit photons from the thermal QGP phase we will use the full LO rate of Ref. <cit.>, which was computed using weak-coupling expansion in a thermal QFT. The rate used is,    


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


with the leading-log coefficient A(p), which is given by 

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


The remaining part of this rate is given by 

    ν_LO(x) ≡   1/2ln(2x) + C_2 ↔ 2(x) 
       + C_brem(x)  + C_annih(x)

with the Fermi distribution function n_f(k)=[exp(k/T)+1]^-1. The dimension of the quark representation is d_F, which is 3 in our case. Summing over the charges of quarks,  q_s, one gets  d_F∑_s q_s^2 = 3 × (1 · (2/3)^2 + 2 · (1/3)^2 ) = 3 × 6/9. The leading-order asymptotic thermal quark-mass m_∞ is given by <cit.> to be


    m^2_∞ = C_F g_s^2T^2/4


with the quadratic Casimir of the quark representation C_F, which is C_F=4/3 for QCD, and the strong coupling g_s = √(4π α_s). Using the running coupling prescription, 

    α_s(Q) = 12π/(11N_c -2N_f)log(Q^2/Λ^2_QCD)

where the cutoff scale, Λ_QCD=0.2. For SU(3), with N_c=3 and three flavours, N_f=3 we get that for ALICE energies, α_s≈ 0.3. The functions that describe the two-to-two particle processes (C_2 ↔ 2) and the in-medium bremsstrahlung (C_brem) and annihilation (C_annih) processes are,

    C_2 ↔ 2   = 0.041 x^-1 - 0.3615 + 1.01e^-1.35x
    
    C_brem   + C_annih≃√(1 + 1/6N_f)
       ×[0.548 log(12.28+1/x)/x^2/3 + 0.133x/√(1+x/16.27)]
 
with x=p/T for three flavours, N_f= 3.
These functions were obtained by approximating the full kinetic kernels. The full logarithm under the log  will also be used to enhance the non-equilibrium rate, with the substitution x=E/t → x'=E/Q. 


 §.§ Photon emission from the hadron resonance gas


For from the hadron resonance gas (HRG) phase, we use the the thermal photon emission rate  the parametrization Ref. <cit.>. These  parametrizations have an underlying error of no more than 20% with the microscopic calculated values <cit.>. We use this parametrization since the inclusion of the full cross section into a phenomenological model is not practical, and very computationally expensive <cit.>. Two different contributions are included, one from the meson channel ππ→ππγ and another one including  the emission from  in-medium ρ mesons. These parametrizations can be applied to photons with energies q_0 between 0.2 and 5, at temperatures between T=100-180 and baryon chemical potentials of μ_B=0-400. For these investigations we will set μ_B=0.

The contribution from in-medium ρ-mesons, including channels like π N →π N γ and N N → N N γ, are universally given by <cit.>,


    q_0dR_γ^ρ/d^3q(q_0;T) = exp[a(T)q_0 + b(T) + c(T)/q_0+0.2] .

Here, and in the following, q_0 and T are given in units of GeV. We will use the fitted  parameters given in Ref. <cit.>

    a(T) = -31.21 + 353.61 T - 1739.4 T^2 + 3105 T^3
          b(T) = -5.513 -42.2T + 333T^2-570T^3
          c(T)=-6.153+57T-134.61T^2+8.31T^3


Nevertheless, this contribution does not include meson-meson bremsstrahlung, strongly dominated by the ππ→ππγ channel. The contribution from  π K scattering is subleading, and will not be included, since it comprises  at most an increase of 20%. The following fit function is used


    q_0dR_γ^Brems/d^3q(q_0;T) =   exp{α_B(T)+ q_0 β_B(T)
    
     +   γ_B q_0^2 + δ_B(T)(q_0+0.2)^-1}

with the following fitted parameters, 


    α_B(T)    = -16.28 + 62.45 T - 93.4T^2 + 7.5 T^3
    β_B(T)    = -35.54 +414.8T -2054T^2+3718.8T^3
    γ_B(T)   =0.7364-10.72T+56.32^2-103.5T^3
    δ_B(T)   =-2.51+58.152T-318.24T^2+610.7T^3

In the HRG, these two contributions are relevant for different kinematic windows of the photons. For a temperature of 150 , soft photons (q_0 < 0.4) are strongly dominated by ππ scattering. On the other hand, the contribution form ρ-meson decays is an order of magnitude larger for q_0 > 1 <cit.>. 



