Probing the evolution of heavy-ion collisions using direct photon interferometry

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 \textit{bottom-up} thermalization scenario. In the second, the thermal rates are enhanced close to the pseudo-critical temperature $T_c\approx 155\,\text{MeV}$ using a phenomenological ansatz. We compute the correlators for relative momenta $q_o, \,q_s$ and $q_l$ for different transverse pair momenta, $K_\perp$, 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 $\sim 10^{10}$ heavy ion events a measurement of the HBT correlations for $K_\perp<1\, \text{GeV}$ is statistically significant.

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 Tc ≈ 155 MeV using a phenomenological ansatz. We compute the correlators for relative momenta qo, qs and q l for different transverse pair momenta, K ⊥ , 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 K ⊥ < 1 GeV is statistically significant.

I. INTRODUCTION
The relativistic nuclear collision experiments explore the physics of dense and hot QCD matter, also known as the quark-gluon plasma (QGP) [1]. 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 multiplerescatterings 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 [2]. 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 [3][4][5][6][7][8].
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 [9]. HBT correlations, originally introduced to measure the radii of stars from the incoming photons [10,11], have been used extensively across physics, from atomic gas correlations in cold atom experiments [12,13], to pion interferometry in heavy ion collisions experiments [14][15][16]. 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 [17][18][19][20][21][22][23][24] and references therein. So far only one measurement in Pb-Pb collisions at √ s NN = 17.3 GeV at the CERN SPS was reported [25]. In view of the upcoming high-luminosity runs at the LHC [26], 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 [27]. Motivated by the idea that thermal rates might be enhanced around the pseudo-critical temperature by confining modes during hadronization [28,29], 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.

II. HBT CORRELATIONS
Quantum statistical effects can be used to understand the spacetime distribution of particle sources [9,[30][31][32] finding the spatial extension of the photon source in the fireball. For this, we use the normalized HBT correlator, 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 Here, p 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, We can further simplify Eq. (1) by splitting the fourpoint 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 † p1,λ1 a † p2,λ2 a p2,λ2 a p1,λ1 a † p1,λ1 a p1,λ1 a † p2,λ2 a p2,λ2 + a † p1,λ1 a p2,λ2 a † p2,λ2 a p1,λ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. [9] 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 where S(q, K) is the Fourier transform of the emission function, 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(q, K) and the photon invariant yield, S(0, p), will be calculated by associating the Wigner function with photon emission rates, that is A. 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) (8) 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 i.e. K 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 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 [9,23]. For the former, the pair momenta itself is taken to be on-shell, K 0 ≈ |K| + O(q 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 q 2 /K 2 , they are controlled by the group velocity β = K/K 0 [9,23]. However, photons are massless, and this expansion will break at |q|/(2|K|) ∼ 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 |q|, with reasonable confidence (see Sec. V). This means that in general |q| and |K| 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 [33], or by computing the curvature of the correlator at q = 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 |q|/(2|K|) 1 cannot be met for all the kinematic regime. Furthermore, the Wigner function in Eq. (6) is generally given for any combination of momenta. The function can be evaluated off-shell [9,23], 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 q dependence in the denominator of Eq. (5), via S(0, p 1,2 ) → S(0, K). The correlator is given in this limit for which we will introduce the commonly used averaging notation [23]: The smoothness approximation is accurate if the curvature logarithm of the one particle distribution is small [34], which is not true for photons at small momentum (see Fig. 1). In reference [23] it was found that the convergence of the smoothness approximated to the full correlator is restricted for values of |q| < 2 |K|. 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 q o and q s , where the other directions of q are set to zero. In this case, the q o and q s direction look Gaussian, and the explored K ⊥ values will be larger than the inverse half width of the correlator, which makes this approximation safe.

B. Homogeneity Radii
We can get a general form of C(q, K) for an arbitrary Gaussian source around the origin in q space For sources with relatively small non-Gaussianities this approximation is still valid, since the perturbations around C increase only at high q values [35]. The halfwidth tensor,R µν =R µν (K), is a function of the pair momentum. To clean the notation, we avoid writing its K dependence. Using the orthogonality relationship, Eq. (11), we can reduce this expression to 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 In this work we only focus on the diagonal of this matrix. While it has been shown that for longitudinally expanding sources the R os term is relevant [36], it is also true that numerically calculating such cross-terms is more computationally complex.
To compute the radii in Eq. (16), we use the the method of moments, which is stable for correlators with strong non-Gaussianities [30]. We use the moments of the true correlator C(q, K) − 1 in relative momentum space, where R −1 is the inverse matrix of Eq. (16). We have defined the distribution function 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 It is important to clarify that this method requires the correlator to be highly localized around q = 0, to give sensible results for the characteristic scale. In other words, the correlator needs to decay faster than a powerlaw. We can use also this method to quantify the deviations from Gaussianity by computing the normalized excess kurtosis, 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.

III. 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 [37][38][39], 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 TeV in 0-20% centrality class. The initial conditions at τ hydro = 0.6 fm were provided by the twocomponent Monte Carlo Glauber model [40]. 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 GeV/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 fm and dτ = 0.2 fm (x max = y max = 25.2 fm). The final time τ max varied depending on the initial conditions, but at least 100 recorded events had τ max ≥ 15.8 fm. 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 [27,41]. The second source is a phenomenological enhancement of the thermal rates near a pseudo-critical temperature T pc , presented first in Ref. [28]. We discuss these and other photon contributions below.

A. 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 N N → γX process can be calculated using perturbative QCD (pQCD) [42], which is then scaled by the number of binary collisions, N coll , via the relation 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. [3]. The fit function is given by the functional form 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.

B. 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. [43]. 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 [44][45][46]. The parametrization used in this work is given explicitely in Appendix A.

C. Photon emission from the hadron resonance gas
For the thermal photon emission rate from the hadron resonance gas phase the parametrizations of Ref. [47] is used. The given parametrizations agree within 20% with the microscopic calculated values. Microscopic calculations have already been performed [48,49], but, as pointed out in [47], 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 AL-ICE, vanishing chemical potential is assumed.

D. Photon production from pre-equilibrium
Using the "bottom-up" thermalization scenario [41], recent estimates [27,50] 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 −1 s , 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, which combines the geometry of the Glauber model [51] with the IP-Sat model [52,53], 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. [27] to be τ th ∼ 2.4 fm 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 fm 1 . From the field theoretical point of view, in this stage, the gluon medium approaches a nonthermal fixed point [55][56][57], where the gluon occupation is given by Here, α S is the strong coupling, and f S is a scaling function, which can be parametrized from the results of clas-sical statistical simulations [55] as follows 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 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 , 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. (23) and the QGP energy density one can obtain where τ 0 is taken to be the spatially averaged saturation scale, Q s , and can be determined parametrically using the method described in Ref. [27].
For the rate, we will use a kinetic rate, generally given by 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. (28). Expanding in the exchange momentum of the mediator and keeping only the leading term one finds the simplified rate [27,58], where α is the electromagnetic coupling, where N c is the number of colors, and c is the gluon liberation factor described in Ref. [59]. 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 L 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. [60] can be found from the smallangle approximated rate. Nevertheless, at the full leading-order (LO) limit of the photon rate, Ref. [43], 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 where x = E/T in the thermal case, and ν LO (x) is given in Eq. (A3). We expect a similar behavior to the preequilibrium 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. (30) in Eq. (29), for x → x = E/Q s .

E. Critical enhancement at late times near Tc
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 [28,61]. 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. [28,29,62,63], where it is conjectured that the thermal rates are enhanced near a pseudo-critical temperature T c , 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 [64]. For the purpose of this paper, however, we choose to model the enhancement factor, h(T ), as follows where the pseudo-critical temperature is set to be T c = 155 MeV. The enhancement parameters are set to be h 0 = 3 and d = 50 MeV. The enhancement factor is tuned such that the enhancement matches the experimental results from the ALICE collaboration, see Fig. 1.

IV. RESULTS
We compute the total yield for the thermal baseline, and include as well the both enhancement scenarios, which can be seen in Fig. 1. The pre-equilibrium photon spectrum shows a structure around p ⊥ ∼ 2.5 GeV. 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 GeV < p ⊥ ∼ 3 GeV, 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. 1, left). On the other hand, applying the enhancement to the thermal rates, Eqs. (31) and (32), 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. (5), 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 TeV, where the average saturation scale is Q 2 s = 2.9 GeV 2 . As expected, the longitudinal curves are the most sensitive to the inclusion of both enhancements which are presented in Fig. 2 for different values of K ⊥ . 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, R l , R o and R s (see Fig. 4), which were computed using the characteristic scale method and the aforementioned correlators.
Just as expected from the correlators, the change in the longitudinal radius, R l , is the largest one. The change induced in R os 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 R os but may undershoot significantly R l . 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. (20), for the three diagonals. A clear hierarchy is found, where q l breaks Gaussianity the most, followed by q o and q s . We find that the sidewards direction is to good approximation Gaussian (see Fig. 5). The non-Gaussianities, as was explained above and in Ref. [23] 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 q l di-rection, since the boosting from longitudinal expansion is largest for the q l 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, R o gets a non zero contribution from β o t = tK ⊥ /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 [9]. As it can be seen in Fig.5, the normalized excess kurtosis can be used as an observable complementary to the radii. This is particularly true for K ⊥ < 0.5 GeV, where the big difference in ∆ l could be used to differentiate the scenarios.

V. 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 ( 3 GeV) [6,7]. Significantly improved data-taking rates in the upcoming LHC runs 3 and 4 make it possi- ble to collect a sample of Pb-Pb collisions corresponding to an integrated luminosity of 10 nb −1 , or O(10 10 ) collision events. In this section we estimate up to what photon pair transverse momentum K ⊥ a direct-photon HBT measurement might be possible. We concentrate on the longitudinal momentum difference q l . Statistical uncertainties for measurements of q o and q s are very similar. For a Gaussian parameterization the correlation function C of direct photons for q o = q s = 0 is given by 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 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 [6].
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 [6]. We parameterize the spectrum by where the inverse slope parameter is set to T = 0.3 GeV, see Fig. 6. From this simple parametrization of the measured direct-photon spectrum we calculate the number N dir p,u of uncorrelated pairs of direct photons per event in a given q l bin. We consider a measurement of C(q l ) in 10 MeV wide bins for |q o | < 30 MeV and |q s | < 30 MeV in various intervals of pair transverse momentum K ⊥ .
The statistical uncertainty of the total number C ·N evt · N inc p,u of pairs of inclusive photons should be much smaller than the number of pairs (C − 1) · N evt · N inc p,u above the uncorrelated background. Here N evt denotes the number of considered Pb-Pb collisions. This corresponds to Neglecting the small √ C term on the left hand-side, the criterion for a significant measurement in the considered bin reads where a photon detection efficiency of ε = pconv × εreco = 0.04 is assumed where pconv = 0.08 and εreco = 0.5 roughly correspond to the photon conversion and reconstruction efficiencies in the photon conversion measurements of the ALICE experiment [6]. The table also shows the ratio s = 2σ inc rel /r 2 γ for these two cases. For a significant measurement s needs to be significantly smaller than unity.
Results for the statistical uncertainty σ inc rel of the measured correlation C(q l ) for inclusive photons for N evt = 10 10 are given in Table I. This table also shows the ratio s = 2σ inc rel /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 I 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 K ⊥ ≈ 1 GeV. For this value of K ⊥ we illustrate the projected statistical uncertainties of C measured for pairs of inclusive photons in black in Fig. 7. For comparison the distribution is also shown in red for K ⊥ ≈ 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 [26] and to study in detail all sources of systematic uncertainties which might affect the measurement.

VI. 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, q o , q s and q l 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, Projected statistical uncertainties for a measurement of C(q l ) in 10 MeV wide bins for two pair transverse momentum ranges 0.45 < K ⊥ < 0.55 GeV (in red) and 0.95 < K ⊥ < 1.05 GeV (in black) in 10 10 Pb-Pb collisions at √ sNN = 2.76 TeV in the centrality range 0-20%. The other components of the momentum difference are constrained to |qo| < 30 MeV, |qs| < 30 MeV. The shown correlation function corresponds to a Gaussian parameterization (Eq. 33) with an arbitrarily chosen radius R l = 2 fm. the late-time enhancement makes the two-photon correlation narrower. From these correlators we extract the HBT radii R l , R o , and R s . 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 R o and R s 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 K ⊥ 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 iden-tify 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.
After the thermalization of the colored medium, photons can be emitted from either a thermalized quarkgluon-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. [43], which was computed using weak-coupling expansion in a thermal QFT. The rate used is, with the leading-log coefficient A(p), which is given by The remaining part of this rate is given by 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 2 s = 3 × (1 · (2/3) 2 + 2 · (1/3) 2 ) = 3 × 6/9. The leading-order asymptotic thermal quark-mass m ∞ is given by [65] to be 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, where the cutoff scale, Λ QCD = 0.2 GeV. 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, 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. [47]. These parametrizations have an underlying error of no more than 20% with the microscopic calculated values [48,49]. We use this parametrization since the inclusion of the full cross section into a phenomenological model is not practical, and very computationally expensive [47]. 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 GeV, at temperatures between T = 100 − 180 MeV and baryon chemical potentials of µ B = 0 − 400 MeV. 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 [47], Here, and in the following, q 0 and T are given in units of GeV. We will use the fitted parameters given in Ref. [47] a(T ) = −31.21 + 353.61T − 1739.4T 2 + 3105T 3 b(T ) = −5.513 − 42.2T + 333T 2 − 570T 3 (A8) c(T ) = −6.153 + 57T − 134.61T 2 + 8.31T 3 Nevertheless, this contribution does not include mesonmeson 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 In the HRG, these two contributions are relevant for different kinematic windows of the photons. For a temperature of 150 MeV , soft photons (q 0 < 0.4 GeV) are strongly dominated by ππ scattering. On the other hand, the contribution form ρ-meson decays is an order of magnitude larger for q 0 > 1 GeV [49].