The emission of electromagnetic radiation from the early stages of relativistic heavy-ion collisions

We estimate the production of electromagnetic radiation (real and virtual photons) from the early, pre-equilibrium, stage of relativistic heavy-ion collisions. The parton dynamics are obtained as a solution of the Boltzmann equation in the Fokker-Planck diffusion limit. The photon and dilepton rates are integrated and the obtained yields are compared with those from standard sources and with available experimental data. Non-equilibrium photon spectra are predicted for Pb+Pb at $\sqrt{s_{\rm NN}} = 5.02$ TeV.


I. INTRODUCTION
The accepted theory of the nuclear strong interaction is Quantum Chromodynamics (QCD), a local gauge theory which admits a spontaneously broken chiral symmetry. This theory has been very successful in describing static properties of strongly interacting systems, and is also able to interpret and predict the outcome of high-energy scattering experiments involving hadronic particles. In spite of all its remarkable successes, there remains much to be learned about QCD. For example, the behaviour of many-body QCD in temperature and density regions far removed from equilibrium is currently the topic of a vibrant research program. The theoretical nature of the transition between degrees of freedom belonging respectively to the partonic and confined phases has only been recently identified as a rapid crossover, occurring at T c ≈ 150 MeV, for zero baryon density [1]). This region is accessible to heavy-ion experiments performed at both the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). Experiments performed at these facilities have revealed an exotic form of matter: the quark-gluon plasma (QGP) [2].
One of the tantalizing properties of the QGP, is that -contrary to early theoretical expectations -it possesses fluid-like characteristics, and therefore can be modelled using relativistic fluid dynamics [3]. This approach has achieved great empirical success, mainly characterized by a quantitative interpretation of the hadronic flow systematics measured in experiments [4]. Modern relativistic hydrodynamics even enables the extraction of the transport parameters of QCD. For instance, 3D approaches now exist [5], and can be used to extract the shear [6] and bulk [7,8] viscosities of QCD matter. Notwithstanding the recent progress in relativistic fluid dynamics, it is fair to write that relativistic heavy-ion collisions can still not be modelled ab initio: hybrid approaches need to be constructed. Those typically consist of an initial state followed by the hydrodynamics phase which ends with a hadronic cascade and kinetic freeze-out, when inter-particle distances exceed mean-free-paths while the interaction volume expands and cools. The cascade stage usually relies on Monte Carlo packages where hadronic species are allowed to collide and interact with each other. A recent example is smash [9].
The quantum nature of the initial state, especially how it evolves towards "hydrodynamization", is the subject of much current research. The short hydro formation times required by modern phenomenological analyses (τ 0 < ∼ 1 fm/c) is a challenge to perturbative approaches and calculations [10]. However, initial conditions computed within the Color Glass Condensate (CGC) framework obtained using the impact parameter dependent saturation model with the classical Yang-Mills evolution of the classical gluon fields [11] have been shown -when combined with a subsequent viscous hydrodynamic evolution -to yield very successful interpretation of the measured azimutal flow distributions [12] and of other related observables [13]. In alternate approaches, the non-equilibrium nature of the initial stages has been captured in several versions of effective kinetic theories based on either the Boltzmann equation [14][15][16][17], or the Kadanoff-Baym equations [18,19].
This work concentrates on the study and the analysis of the very early stages of a relativistic nuclear collision. We will make use of the Boltzmann transport equation, and use the fact that the partonic interactions are dominated by small-angle scattering and low momentum transfer, which permits a diffusion treatment in terms of a Fokker-Planck equation [20][21][22]. Our theoretical treatment of the non-equilibrium ensemble of quarks, antiquarks, and gluons is summarized in the next section. The experimental variable chosen to characterize the initial state needs to be of a penetrating nature, impervious to final state interactions. Two obvious candidates are QCD jets [23] and electromagnetic radiation [24]: here we focus on the latter. Because the electromagnetic interaction is much weaker than the strong interaction that governs the QGP evolution -α/α s 1 -photons and dileptons produced in heavy-ion collisions travel to the detectors es-sentially unscathed. By now, an impressive body of work has been devoted to the calculation and measurement of real and virtual photons, for conditions at both RHIC and the LHC [24], but less attention has been devoted to the electromagnetic emissivity of the pre-hydrodynamics phase. This is one of the purposes of this work: the production of real and virtual photons will be considered in a medium generally out of statistical equilibrium. This paper is organized as follows: the next section outlines the approach used to model the early-time evolving parton distributions, and to obtain the quark, antiquark, and gluon phase space distribution functions. Section III is devoted to details of the calculation of real photon production rates and yields. Section IV contains the calculation of lepton pair production. We compare with data where appropriate. We then devote Section V to a discussion of the results, and conclude.

II. PARTONIC EVOLUTION AND THE BOLTZMANN EQUATION
At the very early stages of heavy-ion collisions, the composition of the medium is dominated by gluon degrees of freedom released from the colliding nuclei. These gluons, whose dynamics initially follow nonlinear field equations [25], are expected to occupy phase space with a probability inversely proportional to the strong coupling constant, 1/α s , and with typical momentum of the order of the saturation scale Q s . As the system expands, kinetic theory becomes applicable for gluons when approximately τ ∼ Q −1 s , and the evolution of gluons is then characterized by a time-dependent phase-space distribution function. Once gluons become on-shell particles and can be treated using kinetic theory, quark degrees of freedom can be introduced accordingly through QCD interactions, e.g., gg ↔ qq. Quark production is ignored when τ Q −1 s . The ingredients of kinetic theory are the phase-space distribution functions, denoted as f g (t, x, p) for gluons and f q (t, x, p) for quarks. The normalization of these functions gives rise to the number density of gluons and quarks respectively. Owing to the charge conjugation symmetry of QCD, distribution functions of quarks and anti-quarks are not distinguished. Similarly, energy density and entropy density are well-defined quantities in kinetic theory, in terms of integrals of f g (t, x, p) and f q (t, x, p) [26]. For later convenience, we define the longitudinal and transverse pressures as, Note that we ignore the quark mass so that E p = | p| and the energy density is related these pressures as = P L + 2P T .

A. The Diffusion Approximation
The evolution of the phase-space distribution function is described by the Boltzmann equation. It is written as where C g and C q are the collision integrals, which in our current work are determined by 2 ↔ 2 scattering processes from QCD. One may further assume that scatterings among quarks and gluons are dominated by those with small angles, so that the collision integral is simplified as a Fokker-Planck diffusion term [20]. In the presence of quarks, an extra source term contributes as well, thus in total one has [21] C where are the effective currents, and are the sources. The constant integrals in the current and source are Note that I c effectively describes the conversion of a quark/anti-quark to a gluon, due to the exchange of a quark/anti-quark with the medium with small momentum. In the above equations, N c and N f denote the number of colors and number of flavors respectively, C F is the square of the Casimir operator of the colour SU (N c ) group in the fundamental representation and is given by . The Coulomb logarithm L is a divergent integral that is related to the strong coupling constant, L ∼ log α −1 s . Alternatively, in realistic simulations the logarithm L can be taken dynamically, if one explicitly quantifies the UV (q max ) and IR (q min ) cutoffs and writes L = log(q max /q min ) [21].
One may verify that the collision integrals conserve energy density 1 . The total number density of quarks and gluons, i.e., n = n g + n q + nq, is also conserved provided that IR gluon modes do not lead to divergence, otherwise there would be a δ(p) in the gluon distribution function corresponding to a gluon Bose-Einstein condensate (BEC) [27]. In the case of elastic 2-to-2 scatterings and Bjorken expansion, the gluon BEC presents as long as the initial gluon occupation exceeds some critical value. However, if the initial gluon occupation is not sufficiently large, the produced gluon BEC is transient and the system eventually approaches local thermal equilibrium. Quarks and gluons are then described using equilibrium distribution functions:the Bose-Einstein distribution for gluons and the Fermi-Dirac distribution for quarks, with temperature T and a finite effective chemical potential µ. In our calculations of pre-equilibrium photons and dileptons, we initialize the system according to the realistic collisions in experiments at RHIC and the LHC (see discussions later in Section II C), for which a gluon BEC is always presents during the pre-equilibrium system evolution.
The obtained equations can be made dimensionless with the help of the momentum saturation scale Q s , by scaling momenta: p → p/Q s . For the evolution time, we take the combined factor 4πα 2 s L as a constant which can be also absorbed to define the dimensionless time parameter t → (4πα 2 s L)tQ s . Note that the effective strength of strong coupling constant is determined as long as this constant factor is specified. In this work we will consider 4πα 2 s L ∼ 1 which implies α s ∼ 0.234.

B. Bjorken expansion of the QGP
The early time evolution of QGP in high energy heavyion collisions is dominated by a longitudinal expansion along the collision beam axis. We describe this system using the Bjorken model, meaning that the system is boost invariant along the collision beam (z-axis) and translationally invariant in directions transverse to the collision beam (x and y). The symmetry becomes apparent if one writes in terms of the proper time τ = √ t 2 − z 2 and the space-time rapidity y = tanh −1 (z/t), such that all physical quantities depend only on τ . Accordingly, this Bjorken symmetry simplifies the Boltzmann equation. The z = 0 slice is of particular interest as τ → t and the distribution function reduces to a function of 1 In a system with one-dimensional Bjorken expansion which we consider throughout this work, both energy density and number density decay as a function of time. The conservation of energy density and number density is reflected in the evolution equation, transverse momentum p ⊥ and longitudinal momentum p z such that The term proportional to 1/t reflects the nature of expansion along z, whose contribution is strong at early times. We initialize the system at t 0 Q s = 1 solely with gluons, as described by the distribution function [28] while f q (t 0 , p) = fq(t 0 , p) = 0. This gluon dominance is inspired by the colour glass picture [29]. The parameter ξ is used to introduce an initial momentum anisotropy which leads to an initial pressure anisotropy. Given Eq. (15) and the definition of pressures in Eq. (1a), it is not difficult to find that, P L /P T < 1 when ξ > 1. Pressure anisotropy is a quantity that characterizes how far an expanding system deviates from an ideal hydrodynamic description. For instance, in dissipative hydrodynamics, the difference of pressures is proportional to the viscous corrections to the stress tensor, P L − P T ∼ η/t [30], where η in this context is the shear viscosity. In our calculations, we consider two types of initial gluon distribution functions with ξ = 1.0 and ξ = 1.5, and we initialize the system evolution with a pressure isotropy P L /P T = 1 and a pressure anisotropy P L /P T ≈ 0.5 respectively [30]. It should be emphasized that, owing to the fast longitudinal expansion, the early-stage evolution tends to drive the system further away from equilibrium, until collisions among quarks and gluons become dominant.
Although, the overall out-of-equilibiurm effect is stronger in the system evolution with ξ = 1.5 comparing to the ξ = 1.0 case, the system evolution rapidly becomes universal. This effect of an attractor solution has been observed in the solution of kinetic theory within the relaxation time aproximation, QCD effective kinetic theory, as well as out-of-equilibrium hydrodynamics [31][32][33]. As a consequence of the attractor solution, out-ofequilibrium system evolution tends to merge into one single, universal path, irrespective of initial conditions. In Fig. 1, the pressure anisotropy P L /P T is plotted as a function of τ T with respect to the nucleus-nucleus collisions carried out at RHIC and the LHC. The effective temperature T is estimated via the Landau's matching condition, i.e., = eq ∝ T 4 . In all these calculations, with initial pressure anisotropy ξ = 1.0 and 1.5, attractor behavior of the system evolution present with a universal curve realized after a short period of time, τ T > ∼ 10. Because of this attractor behavior, the dependence of final results on the switching time to hydrodynamics from kinetic theory can be suppressed. In our calculations, we use τ hydro = 0.4 fm/c as the switching time to hydrodynamics and solve the pre-equiibrium stage evolution for Unlike the initial state pressure anisotropy that varies as a result of quantum fluctuations in the classical gluon field evolution, the initial state gluon occupation is a fixed quantity according to the multiplicity yield in heavy-ion collisions. In Eq. (15), given a specified value of the ξ parameter, the initial gluon occupation is mostly determined by the constant f 0 and the saturation scale Q s . We shall empirically take Q s = 1 GeV at the top RHIC energy, and allow Q s to vary between 1 GeV and 2 GeV at TeV and 5.02 TeV [34] respectively. This freedom is used to explore the sensitivity of results obtained herein to a specific values of the saturation scale.
To determine the value of f 0 , we will use the relation between the initial state entropy per space-time rapidity, dS/dy, deduced from hydrodynamical simulations and the empirical charge particle multiplicity per pseudo- as seen in Fig. 2. The initial state entropy density is calculated using details of the time evolution of hydro calculations that correctly reproduce final state observables [13,35]. Note that the initial time of hydrodynamical evolution would be the final time of the preequilibrium evolution solved from the Boltzmann equation, i.e., τ hydro = 0.4 fm/c. The manifestly linear relation simply follows from the idea that entropy per particle yield in heavy-ion collisions depends little on rapidity [36,37], and the fact that entropy production during the hydrodynamical evolution is subdominant. The constant 7.14 contains the information on the effective microscopic degrees of freedom in the QCD equation of state [38]. It can be extracted from event-by-event hydrodynamical simulations. Results are shown in Fig. 2 for Au+Au collisions at RHIC and Pb+Pb collisions at the LHC in the centrality class 0-20%. These results are obtained using by now standard hydrodynamical simulations of heavy-ion collisions [13], with τ hydro = 0.4 fm/c, and η/s =0. 12 TeV. A temperature dependent ζ/s is used [7]. The observed experimental flow harmonics can be well reproduced, together with other key observables. As shown in Fig. 2, the linear relation between the charged particle multiplicity and initial entropy is indeed apparent, although the slopes from Au+Au and Pb+Pb are very slightly different. Note that the extracted constant is smaller comparing to the simple estimate using a conformal EoS, which is approximately 7.5 [37]. The Eq. (16) allows one to work out the appropriate entropy in the beginning of hydrodynamical simulations, given experimental results for charged particle multiplicities and the access to the details of the hydrodynamical simulations. Values corresponding to the measured charged particle multiplicity at RHIC [39] and the LHC [40,41] are shown in Table I.
The dominant entropy production is from the preequilibrium stage of the system evolution. We solve the Boltzmann equation with respect to initial condition Eq. (15), up to τ hydro = 0.4 fm/c. Entropy density is a well-defined quantity in the kinetic theory. For quarks and gluons, one has which gives the total entropy per space-time rapidity as In realistic calculations, the transverse overlapping area A T can be determined effectively, although there is not a "standard" way to calculate the nuclear overlap area. Since our study pertains to early time dynamics, the geometry associated with the Glauber model is appropriate. Overlap areas calculated using Glauber Monte-Carlo were calculated for different systems colliding at different energies, and binned in different centrality classes. These values were tabulated in Ref. [42], and they are written symbolically here as A T .
In principle, the pre-equilibrium entropy production monotonously depends on the values of f 0 , hence by tuning f 0 , one is able to match the desired entropy per spacetime rapidity, with respect to the realistic colliding systems. Fig. 3 illustrates the evolution of entropy in the pre-equilibrium stage and the matching procedure using various values of f 0 in our calculations. We now have a partonic evolution model which can be used to calculate the emission of electromagnetic radiation for the duration of the non-equilibrium phase.  Table I; see the main text for data references.

III. PHOTON PRODUCTION
To reiterate, electromagnetic probes are penetrating and relay information complementary to that contained in strongly interacting ones. For studies of preequilibrium dynamics their value is matched only by jets in the hard sector, and unparalleled in the soft sector. First, we investigate the production of real photons. In a system comprised of only quarks and gluons, both dileptons and photons can be produced through the annihilation of a quark with an anti-quark. However, unlike dileptons, photons canal be so generated through the Compton scattering process, in which a quark or antiquark scatters with a gluon. Feynman diagrams for both those channels are shown in Figure 4. As was the case for the virtual photon emission treated earlier [43], the real photons attributed to the LPM effect are not considered explicitly here. That contribution can be comparable in magnitude to the sum of Compton and quark antiquark annihilations, depending on the energy of the emitted photon [44].

A. Out-of-Equilibrium Photon Emission Rate
The production rate of photons can be derived starting with the expression for the production of on-shell photons, where the degeneracy factor has been absorbed into the amplitude |M i | 2 . Summation over i denotes contributions from the Compton and quark-antiquark annihilation channels. The distribution functions f 1 , f 2 and f 3 are for quarks and gluons, corresponding to the scattering processes respectively. Following the procedure outlined, for example, in Ref. [45], the expression for the off-equilibrium photon production rate can be derived using the small-angle approximation which assumes the dominance of low momentum transfer between scattering particles. One may thus expand the kinematic variables in terms of the exchanged momentum q = p − p 1 , and perform the kinematic integrals.
Writing explicitly the gluon and fermion distributions, f g and f q (= fq), one arrives at expressions for the photon production rate [45]. The rate from the annihilation process is given by For the Compton scattering contribution, a similar derivation can be performed which yields the expression where the logarithmic divergence is given by The IR cutoff is given by the the Debye mass scale m D ∼ gT and the UV is regulated by the temperature T . Summing the Compton and annihilation contributions gives the expression where I c is defined in Eq. (12). The expansion of the production rate in terms of the exchanged momentum described previously simplifies the analytic expressions, but avoids the details of the HTL (Hard Thermal Loops) regulation of IR divergences [46]. To fix the overall scale of the net rates, the quark and gluon distribution functions, f q and f g , in Eq. (25) are replaced by thermal distribution functions. The factor L is treated as an adjustable constant fixed by matching the result of this expression to that of the equivalent analytical expression from [46]. This procedure yields a constant of ∼ O(1). Importantly, this will neglect the contribution associated with the Landau-Pomeranchuk-Migdal (LPM) effect, which has been shown to contribute an approximate additional factor of 2 [44,47] in equilibrium. However, owing to phase space considerations, the non-equilibrium dynamics may well have a different effect on the LPM contribution than it does on 2 → 2 processes. We consider this possible factor to be part of the systematic theoretical uncertainties in this study, and we argue that leaving it out in fact provides a conservative yield estimate. An appropriate assessment of complete leading order non-equilibrium photon production in the context of the dynamical approach used here requires the numerical implementation of a field-theoretical analysis [48] which we leave for future work. This would include a non-equilibrium assessment of Debye screening and LPM effects.

B. Out-of-Equilibrium Photon Yield
To calculate the off-equilibrium photon yield, the photon production rate is converted using This means that an integral over y and τ is needed in order to obtain an expression of the form dN/dy p d 2 p ⊥ . In the photon production rate calculation, it was assumed that y = 0. Therefore, the y dependence needs to be restored for non-zero values of y. To do this, a change of variables from p z (y) =p z = p ⊥ sinh(y p − y) is performed. This can be further rewritten knowing that p z = p ⊥ sinh(y p ), so that change of variables becomes such that p z (0) = p z returns the original equation. Thus, upon performing this change of variables and noting that the distribution functions are also a function of τ , an integration over y and τ gives The integration over x ⊥ yields A T , the overlapping transverse area of the two colliding nuclei. The expression for the off-equilibrium photon yield can therefore be written as dN In the expression for the yield, the parton distribution functions are obtained from the time-dependent numerical solution of the Boltzmann equation, using the procedures and methods discussed in Section II. Eq. (28) is then integrated from an initial time of 1/Q s , where Q s is either 1 or 2 GeV. For the choice of a final time, we rely on analyses of heavy-ion phenomenology, as commonly practiced. Specifically, many fluid-dynamical simulations of relativistic nuclear collisions adopt a starting time of τ hydro = 0.4 fm/c [3]. Therefore, in this study the pre-equilibrium phase exists for a time τ , with 1/Q s ≤ τ ≤ 0.4 fm/c. As a first step, we study the effect of the anisotropy of the initial gluon distribution on the pre-equilibrium photon spectra. Figure 5 shows the result of integrating the non-equilibrium photon rates (rates evaluated with distributions from the transport calculation) from a time of 1/Q s (∼ 0.2 fm/c for RHIC energies, and ∼ 0.1 fm/c for LHC) to 0.4 fm/c, the end of the pre-equilibrium phase. The effect of using the different values of the asymmetry parameter ξ is shown, and its values are chosen to span the parameter space and sample a of ratios of longitudinal to transverse pressure, as defined by Eqs. (1a).
Recall that the values f 0 follow from requiring the multiplicity per unit pseudo-rapidity to match experimentally measured values, via a fluid-dynamical analysis. One observes that the effect of the initial gluon asymmetry on the final photon spectra is modest. Quantitatively, going from ξ = 1 to ξ = 1.5 we report an increase of ∼35% at RHIC ( √ s N N = 200 GeV and Q s = 1 GeV), and of ∼5% and ∼15% at the LHC (for Q s = 2 GeV, and √ s N N = 2.76 TeV and 5.02 TeV, respectively). These numbers are for a photon transverse momentum of p T = 2.5 GeV.

C. Other sources and experimental data
Electromagnetic radiation is emitted throughout the entire space-time history of the collision process. Therefore, the radiation from the pre-equilibrium sources will compete with others, and will constitute only one of the contributions to the total yield. In what concerns real photons, those other contributors are primordial nucleonnucleon collisions. The photon production there can be calculated using next-to-leading order (NLO) perturbative QCD [49]. These photons are often referred to as "prompt photons". In addition, the strongly interacting medium which is modeled by viscous hydrodynamics will shine throughout its existence [50]. The photons are often referred to as "thermal photons". Finally, the late stages, where matter falls out of thermal equilibrium can also generate photons. Those contributions must be evaluated using a transport approach [51,52]. Fig. 6 shows the net pre-equilibrium photon yield, with the Compton and qq annihilation channels shown separately, for conditions prevalent at RHIC and at the LHC. In both plots, a striking feature is the dominance of the Compton channel over the fermion annihilation channel. This fact is easily understood in terms of the parton dynamics at work here. The Compton channel is linear in the fermionic density, whereas annihilation is quadratic. Initially, the fermions are absent, as the initial state is gluon-dominated; the quark and anti-quark populations then proceed to grow dynamically. This is illustrated in Fig. 7 which shows the time evolution of the gluonic and fermionic parton density. The difference in intensity between the Compton and annihilation channels is therefore a direct consequence of this asymmetry in partonic content. To illustrate this point even more vividly, recall that the photon-producing Compton and quark/anti-quark annihilation rates are identical in equilibrium [53]. The time-evolution of an equilibrated medium would therefore populate the photon final state spectrum with an equal number of "Compton photons" and of "qq photons". Thus, pre-equilibrium photons offer unique insight into the dynamics which control the early-time chemistry of the parton population.
We now turn to other sources and also consider experimental data, in order to set the scale of the early-time photon radiation. Figure 8 shows a variety of sources, together with direct 3 photon data from two different experimental RHIC collaborations. The dotted line represents the photons originating from primordial nucleon-nucleon collisions, the calculation of which considers NLO pQCD contributions, together with corrections accounting for isospin and nuclear in-medium effects [54]. Adding to those the photons generated in the hydrodynamic phase [50] (shown by the dot-dashed line) of the nuclear reac- tion yields the short-dashed contribution. It is important to specify that the "hydro photons" are corrected for viscous effects, in both the shear and bulk sector, as completely as is currently known [50]. Note that, since the late-stage electromagnetic emission models are still under development [52], those photons are simulated by letting the hydrodynamic modeling operate until T ∼ 105 MeV [50]. The solid line is the sum of all photon sources considered in this study. The difference between that line and the short-dashed one represents the pre-equilibrium contribution. At a transverse momentum of p T = 2 GeV, the pre-equilibrium photons represent ∼ 6% of the total yield, and less than 3% at p T = 2.5 GeV. The signature of the pre-equilibrium dynamics go down with increasing momentum, to approximately disappear at p T ∼ 3 GeV. The pre-equilibrium photons are also calculated for LHC conditions, and are plotted in Fig. 9. The different sources here are as described for RHIC. In addition, a study of the effect of varying Q s is performed: the left panel uses Q s = 1 GeV, whereas the right panel uses Q s = 2 GeV. There, one observes a more signifi- cant ∼ 38% contribution from the pre-equilibrium photons to the total signal, at p T = 2.5 GeV. The current data does not have the resolution to exclude either Q s value, and is statistically consistent with both. However, this simple model does generate a pre-equilibrium photon yield which could be within reach of contemporary experiments, provided an upgraded low-p T resolution.
The study performed in this work can be extended to include photon spectra to be measured and analyzed at the LHC, featuring Pb+Pb collisions at an energy of √ s = 5.02 A TeV. This prediction appears in Fig. 10. At p T = 2.5 GeV, pre-equilibrium photons represent ∼ 25% of the total photon yield, and its effect can be seen to persist up to p T ∼ 4 GeV. The pre-equilibrium component shines about as brightly at both LHC energies, but the contribution from the hydro phase increases with increasing energy, therefore outshining the pre-equilibrium contribution more than at the lower colliding energy.

IV. DILEPTON PRODUCTION
The dense system of quarks and gluons formed immediately after relativistic heavy-ion collisions can also be studied using virtual photons, or dileptons produced through quark/anti-quark annihilation. This section contains the details of the derivation.
From relativistic kinetic theory, the rate of production of dileptons from qq → l + l − (Fig. 11), can be derived [58][59][60][61] dR which is the number of dileptons produced per space-time volume and four dimensional momentum-space volume. In this equation, the relativistic relative velocity is and the total cross section is given by where Only u, d, and s massless quarks are used and it is assumed that the rest mass of the leptons is much less than M, the centre-of-mass energy and the dilepton invariant mass. As this is the off-equilibrium case, thermal distribution functions cannot be used. Therefore, the distribution functions must come from the out-of-equilibrium solution to the Boltzmann equation [21], where, as mentioned, the resulting distribution functions are given in terms of p ⊥ , p z , and τ . After integrating over p 2 and φ 1 , the equation becomes , such that the integral is no longer dependent on p 2 . The known identities were used for the boost-invariant case, where the rapidity may be set to y = 0 as in the Bjorken model such that The off-equilibrium dilepton production rate could also be written in terms of mass distribution knowing that where M dM = 1 2 dM 2 and, as before, y = 0. Therefore, by integrating over Q ⊥ , an alternate expression for the rate is given by

B. Out-of-Equilibrium Dilepton Yield
The above expression can be converted into the dilepton yield using the equality Thus, the off-equilibrium dilepton yield can be determined using where d 2 Q ⊥ = 2πQ ⊥ dQ ⊥ and the integration over x ⊥ is simply taken to be the overlapping area of the two colliding nuclei. As for the calculations of real photon production, the overlap area is estimated using Glauber Monte-Carlo results [42]. Finally, the off-equilibrium dilepton yield can be determined from the expression the increase due to the anisotropy is ∼ 15%. Those numbers are for an invariant mass M = 2.5 GeV and for a centrality class 0 − 20%, as reported in Fig. 12.

C. Other sources and experimental data
In the low invariant mass region considered in this study, other contributing dilepton sources are Drell-Yan production from primordial nucleon-nucleon interactions [62]. The pairs radiated from in-medium reactions involving mesons and baryons have traditionally been considered prime sources of information on in-medium properties [63]. The hadrons which freeze-out at the end of the strong interaction era will also emit lepton pair via radiative decay channels. This last source is commonly referred to as "the cocktail". In addition, the leptons coming from semi-leptonic open-charm meson decay can combine and constitute an irreducible background [64] as far as the other sources are concerned. However, if the detector suite has the capability to recognize and analyze displaced decay vertices [65], those sources can be subtracted.
In the case of dileptons, it turns out their revealing potential is substantially less promising, in what concerns the nature of the pre-equilibrium phase. This is situation is due to the fact that the Born graph describing lepton pair production, Figure 11, has an initial state which consists of quarks and anti-quarks only; no gluons at leading order. In a gluon-dominated initial state, the fermions are only produced in secondary processes, as described in Section II. To get a sense of scale, it is useful to compare with the cocktail contribution which defines the threshold for new physics in the dilepton channel. This is shown on Fig. 13 corresponding to pre-equilibrium emission sits orders of magnitude below the cocktail.

V. DISCUSSION AND CONCLUSION
This study has shown that, in the end, the potential for electromagnetic radiation to resolve the anisotropy in the initial gluon distribution is not large. This is illustrated by the photons results of Fig. 5 and by the dileptons results of Fig. 12. We find that this statement holds at both RHIC and the LHC, at least for the systems studied in this work. This is a consequence of the fact that our model assumes an initial parton population which is gluon-dominated, with quarks and anti-quarks appearing dynamically as time evolves, as shown in Fig.   7. The fermions hardly carry a visible imprint of the initial gluon spectrum asymmetry: this is especially visible in the dilepton channel. In addition the global pressure anisotropy is also rapidly quenched, as reported in Fig.  1.
Our estimate of the non-equilibrium photon production has shown that the net yield can contain a measurable signature of the early time dynamics. This is to be contrasted with the findings for lepton pairs, where the pre-equilibrium contribution is suppressed. As mentioned previously, dileptons suffer from the fact that the duration of the pre-equilibrium phase is not sufficiently long to build up an appreciable population of qq pairs. The photons will also suffer from this paucity of fermions, but this more than compensated by the relatively large (in comparison) number of gluons. See Fig. 6. This model calculation therefore suggests that the early, elusive, gluon saturation distribution can inform the final photon spectrum. The influence of a gluon-dominated early hydrodynamic medium on the photon spectrum has been studied in Ref. [67].
One obvious missing element in this work is the contribution of the pre-equilibrium electromagnetic radiation to the so called "direct photon flow puzzle" [50,68]. Naively, an extra early-time component should reduce the overall elliptic flow, as the net v 2 is a weighted average. It is the elliptic flow of each source weighted by the associated photon multiplicity. However, simply ascribing a vanishing elliptic flow to non-equilibrium photons may well neglect the possibly complicated -or even chaotic -dynamics which can influence photons and, to a much lesser extent, hadrons. For instance, the earlytime behavior of the colored plasma is still not well understood, and consequently its modeling is the subject of much debate. The effect on electromagnetic observables of plasma instabilities triggered by the chromoelectric field is still unclear, but recent theoretical developments carrying promising results leave hope for an answer in the near future [69]. Early magnetic fields can also influence electromagnetic spectra, depending on their duration and on their magnitude [70][71][72][73]. In that context, the simple dynamics used in this work do not allow for initial coordinate-space inhomogeneities in the transverse plane. For one-body observables like single particle spectra, this may be acceptable but becomes more questionable for empirical variables that depend on geometry like v 2 , or geometry and fluctuations, like v 3 . These elements all point out the need for a quantitative and realistic spacetime picture of very early relativistic heavy-ion collisions like, for example, 3D IP-Glasma [74], KøMPøST [75], or even magnetohydrodynamics.
Even if the investigation of early time photon spectra, real and virtual, is not yet as popular as that of later photon production, some studies have been devoted to that topic. The pre-equilibrium photon contribution found in this work can be compared with results of the 3D Boltzmann equation simulation approach, BAMPS [76]: an approach close in spirit to what is done here. At RHIC, the photon yields reported here are almost one order of magnitude above those of BAMPS, at low transverse momentum (p T ∼ 1 GeV), they cross at p T ∼ 2.5 GeV. However, BAMPS utilizes pythia, an initial state very different from the simple CGC-like picture used in this work. However, this difference may just be used to prove the point that the photon spectrum is indeed sensitive to details of the initial stages. Work which models the initial state with the Abelian Flux Tube model [77] generates a partonic distribution with the Schwinger mechanism that evolves via 2 → 2 scattering comes close to what is performed here, even though some details differ. That reference also finds imprints of the initial conditions on the final photon spectrum, and the magnitude of the signal is comparable to the effects observed in this work. Similarly, calculations realized in the parton-hadron-string dynamics (PHSD) scenario conclude that the photons (and also dileptons, in that approach) do show characteristic features that can be attributed to specific initial state configurations [78]. It appears that the results reported here for the pre-equilibrium photon spectra are lower than those obtained in a realization of the bottomup thermalization scenario calculated in Ref. [79], by a factor ∼ 2 (LHC), ∼ 3 (RHIC). The bottom-up estimate for the photon yield and the one performed herein cross at p T ∼ 3 GeV. Finally, recent preliminary estimates using KøMPøST to model the pre-equilibrium stage [75] point to a somewhat smaller photon contribution than the one reported here, together with a slight increase in photon elliptic flow. That work is being completed and will appear soon. The partial conclusion which follows from our comparison with other results in the literature is that the evaluation of pre-equilibrium electromagnetic radiation is still a developing effort, and that results obtained so far do not rule out the exciting possibility of an experimental identification.
To conclude, the time-dependent parton dynamics of early-time heavy-ion collisions were modeled with the relativistic Boltzmann equation, solved in the diffusion approximation. It was argued that the dynamics and nonequilibrium emission rates used here make a plausible case for an observable signature of a pre-equilibrium component in the net photon spectra. The magnitude of the real photon signal makes a quantitative statement about kinetic and chemical equilibrium. Indeed, it appears that the fermion suppression observed at early times is largely compensated by the gluon-rich environment in the evaluation of QCD Compton scattering. One early goal of our investigation was to identify a possible electromagnetic signature of a BEC; however this has not been seen in the one-body observables considered here. Work is ongoing to continue these investigations at higher order in α s , to include photon-hadron correlations, and to involve more elaborate simulation approaches.