MARTINI: An event generator for relativistic heavy-ion collisions
aa r X i v : . [ h e p - ph ] O c t MARTINI: An event generator for relativistic heavy-ion collisions
Bj¨orn Schenke, Charles Gale, and Sangyong Jeon
Department of Physics, McGill University, Montreal, Quebec, H3A 2T8, Canada
We introduce the Modular Algorithm for Relativistic Treatment of heavy IoN Interactions(MARTINI), a comprehensive event generator for the hard and penetrating probes in high energynucleus-nucleus collisions. Its main components are a time evolution model for the soft background,PYTHIA 8.1 and the McGill-AMY parton evolution scheme including radiative as well as elasticprocesses. This allows us to generate full event configurations in the high p T region that take intoaccount thermal QCD and QED effects as well as effects of the evolving medium. We present resultsfor the neutral pion nuclear modification factor in Au+Au collisions at RHIC as a function of p T for different centralities, and also as a function of the angle with respect to the reaction plane fornon-central collisions. Furthermore, we study the production of high transverse momentum photonsincorporating a complete set of photon-production channels. I. INTRODUCTION
High transverse momentum jets emerging from thecentral rapidity region in heavy-ion collisions provide im-portant information on the produced hot quark-gluonplasma (QGP). To extract this information from the ex-perimental data, it is important to develop a good theo-retical understanding of the interactions of hard partonswith the medium. Therefore the energy loss of hard par-tons traversing the medium has been under extensive the-oretical investigation. Gluon bremsstrahlung includingthe Landau-Pomeranchuk-Migdal (LPM) [1] effect hasbeen proposed as the dominant mechanism for energy lossand different theoretical formalisms have been developedand applied to describe it [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].Also, binary elastic scattering off thermal partons is po-tentially important for the energy loss and momentumbroadening of high- p T partons. In [13], elastic energyloss was combined with the AMY [10, 11, 12] radiativeenergy loss within the McGill-AMY evolution formalism,and in [14] the description of the collisional processes wasfurther improved.The main goal of these and other calculations is tocreate a quantitative basis for the “tomography” of thehot and dense nuclear medium created in heavy-ion col-lisions. Among the “tomographic variables” is the nu-clear modification factor R AA , which is defined as theratio of the hadron yield in A+A collisions to that inbinary-scaled p+p interactions. A variety of computa-tions that differ significantly in the applied energy lossmechanism can reproduce the measured R AA in centralAu+Au collisions at √ s = 200 AGeV at the RelativisticHeavy-Ion Collider (RHIC), given the present experimen-tal errors (see [15] for a systematic comparison of threedifferent formalisms). For this reason the tomographicusefulness of R AA in central collisions has become ques-tionable [16]. More differential observables, such as R AA as a function of both p T and azimuth in non-central col-lisions [17, 18, 19, 20], can help to improve this situation.Other important observables in heavy-ion collisions areelectromagnetic probes, such as photons and dileptons.Because of the small electromagnetic coupling, once pro- duced, they usually escape the medium without furtherinteraction and thus provide undistorted information onthe early stages of a heavy-ion collision. Photon pro-duction from nuclear collisions at RHIC has been cal-culated in [21, 22, 23] using 1+1 dimensional Bjorkenevolution, 2+1 dimensional, or 3+1 dimensional relativis-tic hydrodynamic evolution for the background, respec-tively. Photons in nuclear collisions come from a varietyof sources, namely direct photons, fragmentation pho-tons, jet-plasma photons and thermal photons. Directphotons are predominantly produced from hard collisionsin the early stage of the heavy-ion collision via annihila-tion and Compton scattering processes. Fragmentationphotons are produced from the surviving high energy jetsafter their passing through the thermal medium. Ther-mal photons have a negligible contribution at high p T and thus are not relevant in this range. On the otherhand, jet-plasma photons from photon radiation q → qγ (bremsstrahlung photons) and jet-photon conversion viaCompton and annihilation processes have been shown tobe very important for the understanding of experimentaldata for photon production in Au+Au collisions at RHIC[21, 22, 23, 24].For a best possible comparison of the theoretical de-scription with experimental data, we incorporate theMcGill-AMY formalism [10, 11, 12, 13, 25] for radiativeenergy loss as well as elastic processes [14] into a newevent generator, dubbed Modular Algorithm for Rela-tivistic Treatment of heavy IoN Interactions (MARTINI).Its main ingredients are PYTHIA 8.1 [26, 27] to gener-ate the hard partons from the individual nucleon-nucleoncollisions and to take care of the final fragmentation intohadrons, the evolution of the background medium (e.g.from hydrodynamic models), and the McGill-AMY par-ton evolution formalism. This way it is possible to studyhard observables in heavy-ion collisions theoretically onan event-by-event basis, keeping information on correla-tions.Several Monte Carlo simulations for heavy-ion colli-sions have been or are being developed [28, 29, 30, 31,32, 33, 34, 35, 36]. The implementation of medium ef-fects and use of approximations varies significantly be-tween the different models. MARTINI is the first MonteCarlo simulation to include AMY bremsstrahlung com-bined with elastic processes and in-medium photon pro-duction. It is very flexible due to its ability to incor-porate any soft background evolution that provides in-formation on the temperature and flow velocities of themedium. Currently, hydrodynamic evolution calcula-tions from four different groups [37, 38, 39, 40, 41, 42]have been implemented.In this work we present first results on hard observablesobtained with MARTINI. These include the azimuth av-eraged R AA as a function of p T , R AA as a function ofthe azimuth and p T , as well as photon yields and photon R AA .This paper is organized as follows. We review theMcGill-AMY evolution formalism and present transitionrates for radiative and elastic processes in Section II. TheMonte Carlo simulation is introduced in Section III, re-sults are presented in Section IV, and we conclude andpresent an outlook in Section V. II. THE MCGILL-AMY EVOLUTIONFORMALISM
At the core of MARTINI lies the McGill-AMY formal-ism for jet evolution in a dynamical thermal medium.Here, we briefly review this formalism and discuss its im-plementation into MARTINI in the following section.In previous works [13, 23, 25] the evolution of the jetmomentum distribution in a hydrodynamic backgroundwas computed within the McGill-AMY scheme. This evo-lution is governed by a set of coupled Fokker-Planck typerate equations of the form dP ( p ) dt = Z ∞−∞ dk (cid:18) P ( p + k ) d Γ( p + k, k ) dk − P ( p ) d Γ( p, k ) dk (cid:19) , (1)where d Γ( p, k ) /dk is the transition rate for processeswhere partons of energy p lose energy k . The k < T ≫ gT ≫ g T emerges. This allows radiativetransition rates to be calculated by means of integralequations [12], which correctly reproduce both the Bethe-Heitler and the LPM results in their respective limits. They are given by [21, 43]: d Γ dk ( p, k ) = C s g πp ± e − k/T ± e − ( p − k ) /T × − x ) x (1 − x ) q → qgN f x +(1 − x ) x (1 − x ) g → q ¯ q x +(1 − x ) x (1 − x ) g → gg × Z d h (2 π ) h · Re F ( h , p, k ) , (2)where g is the strong coupling constant ( α s = g / (4 π )), N f is the number of flavors, C s is the quadratic Casimirrelevant for the process, and x ≡ k/p is the momentumfraction of the radiated gluon (or the quark, for the case g → q ¯ q ). h ≡ ( k × p ) × e k determines how non-collinearthe final state is, where e k is the unit vector in a directionnearly collinear with k , p and p + k that can be fixedby convention (see [12]). Parametrically h is of order gT and hence small compared to p · k . Therefore h can be taken as a two-dimensional vector in transversespace. F ( h , p, k ) is the solution of the following integralequation [21, 43]:2 h = iδE ( h , p, k ) F ( h ) + g s Z d q ⊥ (2 π ) C ( q ⊥ ) × n ( C s − C A / F ( h ) − F ( h − k q ⊥ )]+( C A / F ( h ) − F ( h + p q ⊥ )]+( C A / F ( h ) − F ( h − ( p − k ) q ⊥ )] o . (3)Here δE ( h , p, k ) is the energy difference between the finaland the initial states: δE ( h , p, k ) = h pk ( p − k ) + m k k + m p − k p − k ) − m p p , (4)and m are the medium induced thermal masses. C ( q ⊥ )is the differential rate to exchange transverse (to the k direction) momentum q ⊥ . In a hot thermal medium, itsvalue at leading order in α s is [44] C ( q ⊥ ) = m D q ⊥ ( q ⊥ + m D ) , m D = g s T N c + N f ) . (5)For the case of g → q ¯ q , ( C s − C A /
2) should appear asthe prefactor on the term containing F ( h − p q ⊥ ) ratherthan F ( h − k q ⊥ ).Transition rates for elastic processes can also be com-puted in perturbative QCD. They are given by [14] d Γ dω ( p, ω, T ) = d q Z d q (2 π ) Z d q ′ (2 π ) π pp ′ qq ′ × δ ( p − p ′ − ω ) δ ( q ′ − q − ω ) × |M| f ( q, T )(1 ± f ( q ′ , T )) , (6)where p = | p | and p ′ = | p ′ | are the absolute val-ues of the three-momenta of the incoming and outgo-ing hard parton, respectively, and q = | q | and q ′ = | q ′ | are those of the incoming and outgoing thermal parton. ω = p − p ′ = q ′ − q is the transferred energy. The distribu-tion functions f are either Fermi-Dirac or Bose-Einsteindistributions depending on the nature of the thermal par-ton involved. The + or − sign appears accordingly, with − for Pauli blocking and + for Bose enhancement, and d q describes the degeneracy of the thermal parton. In-tegrations over q and q ′ in Eq. (6) can be rewritten interms of the transferred momentum k . Then, employingan appropriate separation of the regimes of soft and hardmomentum exchange, and doing a hard thermal loop re-summation in the soft regime to cure the infrared diver-gence allows for a numerical evaluation of Eq. (6). It iseven possible to extend the expression for the soft regimeto all momenta [14, 45], which leads to a good approxi-mation of the total result and avoids the introduction ofan intermediate scale. The transition rate as a functionof both the transferred energy and momentum d Γ /dωdk is obtained from Eq. (6) after rewriting and omittingthe integration over the transferred momentum k . Us-ing d Γ /dωdk will allow us to sample both the transferredenergy and the transferred three-momentum in an elasticprocess.Both radiative and elastic transition rates were com-puted and tabulated as functions of the parton energy p ,as well as the radiated energy k or transferred energy ω and momentum k , respectively, and are sampled withinMARTINI. Note that since the radiated partons are on-shell ω = k , while the transferred parton in an elasticcollision can be far off-shell, so that ω and k are treatedseparately. As opposed to the radiative transition rates,the dependence on the coupling of the elastic rates isnon-trivial and the rates are also tabulated as functionsof the coupling α s .Apart from the processes described above we includegluon-quark and quark-gluon conversion due to Comptonand annihilation processes, as well as the QED processesof photon radiation q → qγ and jet-photon conversion.The transition rate for the photon radiation process iscalculated analogously to the gluon radiation [10, 12],and in the limit E ≫ T the transition rate for quark-photon conversion is given by d Γ conv q → γ dω ( E, ω ) = (cid:16) e f e (cid:17) πα e α s T E (cid:18)
12 ln
ETm q − . (cid:19) δ ( ω ) , (7)where m q is the thermal quark mass, m q = g T / e f is the charge of a quark with flavor f , α e is the finestructure constant, and the delta function indicates thatwe neglect the energy loss of the quark during the processof conversion. The conversion rate d Γ conv q → g /dω is equal to that in Eq. (7) times a color factor C F = 4 /
3, and d Γ conv g → q dω = N f N c N c − d Γ conv q → g dω , (8)which follows from interchanging the initial with the finaljet. III. MONTE CARLO SIMULATION
MARTINI solves the rate equations (1) using MonteCarlo methods. Instead of evolving a probability distri-bution P ( p ), every event and every hard parton is treatedindividually. This way we obtain information on the fullmicroscopic event configuration in the high momentumregime, including correlations, which allows for a very de-tailed analysis and offers a direct interface between the-ory and experiment. The average over a large number ofevents will correspond to the solution found by solving(1) for the probability distribution.We follow the evolution of one event to describe thebasic functionality of the simulation. First, the numberof individual nucleon-nucleon collisions that produce par-tons with a certain minimal transverse momentum p min T is determined from the total inelastic cross-section, pro-vided by PYTHIA. The initial transverse positions r ⊥ ofthese collisions are determined by the initial jet densitydistribution P AB ( b, r ⊥ ), which for A+B collisions withimpact parameter b is given by P AB ( b, r ⊥ ) = T A ( r ⊥ + b / T B ( r ⊥ − b / T AB ( b ) . (9)Here we use a Woods-Saxon form for the nuclear den-sity function, ρ ( r ⊥ , z ) = ρ / [1 + exp( r − Rd )], to evaluatethe nuclear thickness function T A ( r ⊥ ) = R dzρ A ( r ⊥ , z )and the overlap function of two nuclei T AB ( b ) = R d r ⊥ T A ( r ⊥ ) T B ( r ⊥ + b ). The values of the parame-ters R = 6 .
38 fm and d = 0 .
535 fm are taken from [46].The initial parton distribution functions can be selectedwith the help of the Les Houches Accord PDF Interface(LHAPDF) [47]. We include nuclear effects on the partondistribution functions using the EKS98 [48] or EPS08 [49]parametrization, by user choice. In practise we samplethe initial positions of nucleons in nucleus A and B, su-perimpose the transverse areas depending on the impactparameter and then divide the overlap region into circleswith area σ inel . In three dimensions these are tubes andwe determine how many jet events with given p min T occurwithin such a tube using the number of nucleons fromA and B in the tube, and the probability for a jet event σ jet ( p min T ) /σ inel for each of their combinations.The soft medium is described by hydrodynamics orother models, which provide information on the sys-tem’s local temperature and flow velocity β . Before thismedium has formed, i.e., before the hydrodynamic evo-lution begins ( τ < τ ), the partons shower as in vac-uum. Currently we include the complete vacuum shower,because there is no apparent reason why the vacuumsplittings should end immediately once the medium hasformed. Since most of the vacuum shower occurs be-fore the medium has formed, this is a reasonable approx-imation. We have also tested the other extreme casewhere we stop the vacuum evolution at the virtualityscale Q min = p p T /τ , determined by the time τ atwhich the medium evolution begins. This modifies theparton distribution such that the strong coupling con-stant has to be chosen approximately 10% larger to de-scribe the pion R AA . At this stage, we do not includefinal state radiation (FSR) of the partons that have leftthe medium - all vacuum showers take place before themedium evolution and there is no interference betweenthe medium and vacuum showers. The improvement ofthis point is the subject of future work. Once τ > τ for acertain parton, its in-medium evolution begins. The par-tons move through the background according to their ve-locity. To compute probabilities for any one of the abovedescribed in-medium processes, we perform a Lorentz-boost into the rest-frame of the fluid cell at the parton’sposition and determine the transition rates according tothe local temperature and the parton’s energy in this lo-cal rest-frame. The probability for a parton to undergoany process during a time step of length ∆ t is given by∆ t Γ( p, T ) total , where Γ( p, T ) total = R dkd Γ total /dk forthe radiative processes, and the integral over both ω and k for the elastic processes. d Γ total /dk is the sum over allpossible transition rates, which include those for photonproduction.In case that some process occurs, we decide which oneit is according to the relative weights of the different pro-cesses at the given temperature and parton energy. Wesample the radiated or transferred energy from the tran-sition rate of the occurring process using the rejectionmethod. In case of an elastic process, we also samplethe transferred transverse momentum, while for radiativeprocesses we assume collinear emission for now, which is agood approximation in the weak coupling limit since theemission angle is of order g [12]. After energy and mo-mentum are transferred, we boost back into the labora-tory frame, where the parton continues to move along itstrajectory. Radiated partons are also further evolved iftheir momentum is above a certain threshold p min , whichcan be set (typically we choose p min ≃ T , where T is the local temperature. As both radiative and elas-tic transition rates were computed in the limit E ≫ T ,it would not make sense to attempt to learn anythingabout the evolution of partons with E ≃ T within thisformalism. For partons that stay above that threshold,the evolution ends once they enter the hadronic phase ofthe background medium. In the mixed phase, processesoccur only for the QGP fraction.When all partons have left the QGP phase, hadroniza-tion is performed by PYTHIA, to which the complete -10 -8 -6 -4 -2
2 4 6 8 10 12 14 / ( π p T ) d N π / ( d y dp T ) [ G e V - ] p T [GeV] PHENIX Au+Au 200 GeV 0-10%PHENIX p+p 200 GeVAu+Au AMY+elastic + 3+1D hydro α s =0.3Au+Au AMY+elastic + 2+1D hydro α s =0.33p+p FIG. 1: (Color online) Neutral pion spectrum in p+p and0 −
10% central Au+Au collisions at RHIC energy comparedto data by PHENIX [18, 50]. information on all final partons is passed. BecausePYTHIA uses the Lund string fragmentation model[51, 52], it is essential to keep track of all color stringsduring the in-medium evolution. This includes gener-ating new strings when a gluon is emitted or during aconversion process. In the latter case, new string endshave to be attached to thermal partons, whose momentaare sampled from Fermi or Bose distributions. Being pro-vided with such consistent information on the color stringstructure, PYTHIA is able to perform the fragmentation.This concludes the evolution of one event.The concept of MARTINI is modular, such that we canturn on and off different processes independently, and usedifferent hydrodynamic or other data inputs. In principleit is also possible to extend MARTINI to use differentenergy loss formalisms. The parameters are set using thesame interface as PYTHIA and all options for PYTHIAcan still be modified by the user.
IV. RESULTSA. Pion production and nuclear modification factor π R AA ( p T ) We begin by showing the MARTINI results for thespectrum of neutral pions in p+p collisions at RHICenergy ( √ s = 200 GeV) as well as in central (0-10%)Au+Au collisions compared to data by PHENIX [18,50] in Fig. 1. The calculations were performed usingCTEQ5L parton distribution functions [53] including nu-clear shadowing effects using the EKS98 parametrization[48]. In the shown results for pion production, isospineffects were ignored. We checked that they have no bigeffect (less than 5%), however, we included them in calcu-lations of photon production where they become impor-tant (15-20%). Au+Au calculations take into account T [GeV]PHENIX Au+Au 200 GeV, 0-10% centralAMY + 3+1D Hydro α s =0.39 b=2.4 fmAMY + elastic + 3+1D Hydro α s =0.3 b=2.4 fm 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 π R AA PHENIX Au+Au 200 GeV, 20-30% centralAMY + 3+1D Hydro α s =0.39 b=7.5 fmAMY + elastic + 3+1D Hydro α s =0.3 b=7.5 fm FIG. 2: (Color online) The neutral pion nuclear modificationfactor for mid-central (upper panel) and central (lower panel)collisions at RHIC with √ s = 200 GeV. MARTINI results( b = 2 . b = 7 . α s compared to PHENIX data from [54]. both radiative and elastic processes in the medium de-scribed by either the 2+1 dimensional hydrodynamics of[38, 39, 40] or the 3+1 dimensional hydrodynamics of[41], using a coupling constant α s = 0 .
33 or α s = 0 . α s was adjusted to describe the experimen-tal measurement of the neutral pion nuclear modificationfactor R AA in most central collisions (see below). Forboth hydro evolutions τ = 0 . /c . PYTHIA param-eters were tuned to fit the p+p data, however there isstill freedom and possibly an even better set of param-eters can be found. Using the same PYTHIA settings,we find very good agreement for both p+p and Au+Auspectra. Note that in order to achieve higher statisticswe used certain cutoffs for the minimal p T produced in anucleon-nucleon collision in PYTHIA for the Au+Au cal-culation. The shown result is a combination of runs withminimal p min T = 1 . p T ≤ . p min T = 2 GeVfor 5 . ≤ p T ≤ . p min T = 4 GeV abovethat. Note that the shown results are averages over typ-ically ∼ · simulated events.In Fig. 2 we present the results for the neutral pionnuclear modification factor, defined by R AA = 1 N coll ( b ) dN AA ( b ) /d p T dydN pp /d p T dy , (10)in Au+Au collisions at RHIC measured at mid-rapidityin two different centrality classes (0-10%) and (20-30%),employing the corresponding average impact parameters,2 . . T [GeV]PHENIX Au+Au 200 GeV, 0-10% centralAMY + elastic + 2+1D Hydro α s =0.33 b=2.4 fm 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 π R AA PHENIX Au+Au 200 GeV, 20-30% centralAMY + elastic + 2+1D Hydro α s =0.33 b=7.5 fm FIG. 3: (Color online) The neutral pion nuclear modificationfactor for mid-central (upper panel) and central (lower panel)collisions at RHIC with √ s = 200 GeV. MARTINI results( b = 2 . b = 7 . measurements by PHENIX [54]. All parameters are thesame as in the calculation shown in Fig. 1, i.e., hydro-dynamic background evolution from [41] and α s = 0 . α s is used in all following calculations.(There are no additional parameters for the later calcula-tion of the azimuthal dependence of R AA or the high- p T photon production.) We find very good agreement withthe data for both centrality classes.The value for α s for which the data is described if weonly include radiative processes is α s ≃ .
39, higher thanfor the calculations in [13], where it was 0 .
33. Differencesin the treatment of vacuum showers and the fragmenta-tion scheme are responsible for this difference. We haveobserved that when varying parameters in PYTHIA, suchas the factorization or renormalization scale, or param-eters in the fragmentation function, we need slighty dif-ferent α s to describe the data. This freedom leads to anapproximately 10 - 20 % uncertainty in α s . The valueof α s is again closer to that used in [13] ( α s = 0 . R AA .Importantly, our approach has the potential to revealthe effect of different soft background models on hardobservables. Fig. 3 shows the result for R AA using the2+1 dimensional hydrodynamic evolution from [38, 39,40]. As initial temperatures are typically lower than inthe 3+1 dimensional hydro evolution, we need to increase α s to 0 .
33 to describe the neutral pion R AA . Furthereffects are demonstrated in the following section. B. Angular dependence of the nuclear modificationfactor π R AA ( p T , ∆ φ ) In Figs. 2 and 3 the neutral pion R AA at midrapid-ity is averaged over the azimuthal angle with respect tothe reaction plane φ . To learn more about the producedmedium, e.g. to disentangle the effects of the densityof the medium and the pathlength traversed on the en-ergy loss, R AA has been studied at midrapidity in non-central collisions not only as a function of p T but alsoas a function of the azimuth φ in [18, 19, 20]. We per-form the same analysis within MARTINI and separatethe azimuth in 6 bins of 15 ◦ each, reaching from 0 − ◦ (in-plane) to 75 − ◦ (out-of-plane), as shown in Fig. 4. R AA is then determined in each bin separately. We em-ploy both 2+1 and 3+1 dimensional hydrodynamic back-grounds and α s = 0 .
33 and α s = 0 .
3, respectively. Fig. 5shows the most extreme cases (0 − ◦ and 75 − ◦ ) at b = 7 . in planeout of plane FIG. 4: (Color online) Binning in azimuth ∆ φ for the calcu-lation of R AA ( p T , ∆ φ ) in non-central collisions. Particularly at lower p T , we find a less pronounceddifference between the in-plane and out-of-plane resultthan the data, as do previous theoretical calculations[19, 55]. Furthermore, one can see that the result iscloser to the experimental data for the 3+1 dimensionalhydrodynamic background. This is due to a larger ini-tial eccentricity in this case. Generally, we do not expectto describe the low p T part correctly, because we do notincorporate soft physics such as flow effects or recombi-nation. The too small difference for p T > ∼ T [GeV] AMY+elastic + 3+1D hydro, b=7.5 fm, ∆φ =0-15 ° AMY+elastic + 3+1D hydro, b=7.5 fm, ∆φ =75-90 ° π R AA ( p T , ∆ φ ) AMY+elastic + 2+1D hydro, b=7.5 fm, ∆φ =0-15 ° AMY+elastic + 2+1D hydro, b=7.5 fm, ∆φ =75-90 ° PHENIX 20-30% ∆φ =0-15 ° PHENIX 20-30% ∆φ =75-90 ° FIG. 5: (Color online) Nuclear modification factor as a func-tion of p T and ∆ φ for mid-central collisions at √ s = 200 GeV.MARTINI results ( b = 7 . α s = 0 .
33; lower panel 3+1D hydro, α s = 0 .
3) with bothradiative and elastic processes for most in-plane and mostout-of-plane angular bins compared to PHENIX data from[20]. The gray band indicates the experimental error in R AA . C. Photon production
Photons in the high p T region produced in nuclearcollisions are dominantly direct photons, fragmentationphotons, and jet-plasma photons. Direct photons areincluded in PYTHIA. Apart from leading order directphotons, PYTHIA produces additional photons emittedduring the vacuum showers. Parts of these overlap witheffects that are found in NLO calculations [61, 62], butthere is no simple theoretical way to identify the amountof overlap between the two. Also fragmentation pho-tons are part of the photons produced in the showersin PYTHIA. We computed photon production withinPYTHIA with and without photons from showers andfound that the shower contribution leads to an effective K -factor of approximately 1.8 in the regarded p T range.For heavy-ion reactions MARTINI adds the very rele-vant jet-medium photons from photon radiation q → qγ (bremsstrahlung photons) and jet-photon conversion viaCompton and annihilation processes. As above, we in-clude the full vacuum shower in the shown calculation -most of the shower photons will be emitted before themedium has formed and the parton has realized that ithas formed. To improve on this approximation and in-clude interference between vacuum and in-medium radi-ation is the subject of future work.We present results for photon production in p+p andAu+Au collisions at √ s = 200 GeV compared to data by -12 -10 -8 -6 -4 -2
0 2 4 6 8 10 12 14 16 / ( π p T ) d N γ / ( d y dp T ) [ G e V - ] p T [GeV] PHENIX Au+Au 200 GeV 0-10%Au+Au dir.+ISR+FSR+jet-medium α s =0.3PHENIX p+p 200 GeVp+p direct+ISR+FSR FIG. 6: (Color online) Photon yield in p+p and Au+Au col-lisions at RHIC energy √ s = 200 GeV. MARTINI results( b = 2 . PHENIX [58, 59, 60] in Fig. 6.Another observable that provides information on theeffect of the nuclear medium on photon production is thephoton nuclear modification factor R γAA = 1 N coll ( b ) dN γAA ( b ) /d p T dydN γpp /d p T dy . (11)Fig. 7 shows R γAA as a function of p T for most centralAu+Au collisions ( b = 2 . −
10% central PHENIX data. We find that in the ap-proximation of including all photons from the vacuumshower leads to good agreement with the data within theerror bars. We checked that not including any photonsfrom the final state vacuum shower before τ , which cor-respond to a pre-equilibium contribution, reduces R AA by approximately 20%. V. CONCLUSIONS AND OUTLOOK
We presented first results obtained with the newly de-veloped Modular Algorithm for Relativistic Treatmentof heavy IoN Interactions (MARTINI). This hybrid ap-proach describes the soft background medium using hy-drodynamics or other medium models and simulates thehard event microscopically, using PYTHIA 8.1 to gener-ate individual hard nucleon-nucleon collisions. Hard par-tons are evolved through the medium using the McGill-AMY evolution scheme including radiative and elas-tic processes. Fragmentation is performed employingPYTHIA 8.1 which uses the Lund string fragmentationmodel. Apart from parameters in PYTHIA which werefixed by matching the neutral pion and photon spectrain p+p collisons to experimental data, α s is the only freeparameter. Employing the 3+1 dimensional hydrody-namic evolution from [41], it was set to α s = 0 . R AA γ p T [GeV] PHENIX Au+Au 200 GeV, 0-10% centralMARTINI AMY+el. + 3+1D hydro α s =0.3 FIG. 7: (Color online) Photon nuclear modification factorfor central collisions at √ s = 200 GeV. MARTINI resultscompared to data from [60]. See main text for details. match the neutral pion R AA measurement for central col-lisions. Using the same value for all other calculations(there was no additional freedom in any of the calcula-tions), we were able to describe the neutral pion R AA inmid-central collisions as well as photon yields and photon R AA . For photon production, we added photons from jet-medium interactions to the prompt and shower photonsthat PYTHIA produces.We have demonstrated that MARTINI is able to re-veal the effect of using different evolution models for thesoft background on results for hard observables. For ex-ample, when using a 2+1 dimensional hydrodynamicalevolution model [38, 39, 40], a slightly different value for α s ( α s = 0 .
33) had to be employed to describe the ex-perimental data. The calculated azimuthal dependenceof the neutral pion R AA in mid-central collisions turnedout to be weaker than that found by PHENIX, for bothshown hydrodynamic models, but more so for the 2+1dimensional one. This is mainly due to a larger initialeccentricity in the 3+1 dimensional model. The result forthe 3+1 dimensional hydrodynamic background is alongthe line with previous calculations using different energyloss schemes [19, 55]. The weak azimuthal dependence,particularly for p T < ∼ Acknowledgments
We are happy to thank Steffen Bass, EvanFrodermann, Ulrich Heinz, Scott Moreland, ChihoNonaka, and Huichao Song for very useful correspon-dence. We thank Torbj¨orn Sj¨ostrand for very help-ful clarifications regarding PYTHIA 8.1, and grate- fully acknowledge fruitful discussions with Guy Moore,Guang-You Qin, Thorsten Renk, Michael Strickland, andVasile Topor-Pop. This work was supported in part bythe Natural Sciences and Engineering Research Coun-cil of Canada. B.S. gratefully acknowledges a RichardH. Tomlinson Fellowship awarded by McGill University. [1] A. B. Migdal, Phys. Rev. , 1811 (1956).[2] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne,and D. Schiff, Nucl. Phys.
B483 , 291 (1997), hep-ph/9607355.[3] A. Kovner and U. A. Wiedemann (2003), hep-ph/0304151.[4] B. G. Zakharov, JETP Lett. , 952 (1996), hep-ph/9607440.[5] M. Gyulassy, P. Levai, and I. Vitev, Nucl. Phys. B594 ,371 (2001), nucl-th/0006010.[6] X.-N. Wang and X.-f. Guo, Nucl. Phys.
A696 , 788(2001), hep-ph/0102230.[7] B.-W. Zhang and X.-N. Wang, Nucl. Phys.
A720 , 429(2003), hep-ph/0301195.[8] A. Majumder, E. Wang, and X.-N. Wang, Phys. Rev.Lett. , 152301 (2007), nucl-th/0412061.[9] A. Majumder and B. Muller, Phys. Rev. C77 , 054903(2008), arXiv:0705.1147.[10] P. Arnold, G. D. Moore, and L. G. Yaffe, JHEP , 009(2001), hep-ph/0111107.[11] P. Arnold, G. D. Moore, and L. G. Yaffe, JHEP , 057(2001), hep-ph/0109064.[12] P. Arnold, G. D. Moore, and L. G. Yaffe, JHEP , 030(2002), hep-ph/0204343.[13] G.-Y. Qin et al., Phys. Rev. Lett. , 072301 (2008),arXiv:0710.0605.[14] B. Schenke, C. Gale, and G.-Y. Qin, Phys. Rev. C79 ,054908 (2009), arXiv:0901.3498.[15] S. A. Bass et al., Phys. Rev.
C79 , 024901 (2009),arXiv:0808.0908.[16] T. Renk, Phys. Rev.
C74 , 034906 (2006), hep-ph/0607166.[17] A. Majumder, Phys. Rev.
C75 , 021901 (2007), nucl-th/0608043.[18] S. S. Adler et al. (PHENIX), Phys. Rev.
C76 , 034904(2007), nucl-ex/0611007.[19] R. Wei (PHENIX) (2009), arXiv:0907.0024.[20] S. Afanasiev et al. (PHENIX) (2009), arXiv:0903.4886.[21] S. Turbide, C. Gale, S. Jeon, and G. D. Moore, Phys.Rev.
C72 , 014906 (2005), hep-ph/0502248.[22] S. Turbide, C. Gale, E. Frodermann, and U. Heinz, Phys.Rev.
C77 , 024909 (2008), arXiv:0712.0732.[23] G.-Y. Qin, J. Ruppert, C. Gale, S. Jeon, and G. D. Moore(2009), arXiv:0906.3280.[24] R. J. Fries, B. Muller, and D. K. Srivastava, Phys. Rev.Lett. , 132301 (2003), nucl-th/0208001.[25] G.-Y. Qin et al., Phys. Rev. C76 , 064907 (2007),arXiv:0705.2575.[26] T. Sj¨ostrand, S. Mrenna, and P. Skands, Comput. Phys.Commun. , 852 (2008), arXiv:0710.3820.[27] T. Sj¨ostrand, S. Mrenna, and P. Skands, JHEP , 026(2006), hep-ph/0603175. [28] X.-N. Wang and M. Gyulassy, Phys. Rev. D44 , 3501(1991).[29] A. Dainese, C. Loizides, and G. Paic, Eur. Phys. J.
C38 ,461 (2005), hep-ph/0406201.[30] I. P. Lokhtin and A. M. Snigirev, Eur. Phys. J.
C45 , 211(2006), hep-ph/0506189.[31] I. P. Lokhtin et al., Comput. Phys. Commun. , 779(2009), arXiv:0809.2708.[32] S. Wicks (2008), 0804.4704.[33] T. Renk, Phys. Rev.
C78 , 034908 (2008),arXiv:0806.0305.[34] T. Renk (2008), arXiv:0808.1803.[35] K. Zapp, G. Ingelman, J. Rathsman, J. Stachel, andU. A. Wiedemann, Eur. Phys. J.
C60 , 617 (2009),arXiv:0804.3568.[36] N. Armesto, L. Cunqueiro, and C. A. Salgado (2009),arXiv:0907.1014.[37] K. J. Eskola, H. Honkanen, H. Niemi, P. V. Ruuskanen,and S. S. Rasanen, Phys. Rev.
C72 , 044904 (2005), hep-ph/0506049.[38] P. F. Kolb, J. Sollfrank, and U. W. Heinz, Phys. Rev.
C62 , 054909 (2000), hep-ph/0006129.[39] P. F. Kolb and R. Rapp, Phys. Rev.
C67 , 044903 (2003),hep-ph/0210222.[40] P. F. Kolb and U. W. Heinz (2003), nucl-th/0305084.[41] C. Nonaka and S. A. Bass, Phys. Rev.
C75 , 014902(2007), nucl-th/0607018.[42] U. W. Heinz, J. S. Moreland, and H. Song (2009),arXiv:0908.2617.[43] S. Jeon and G. D. Moore, Phys. Rev.
C71 , 034901 (2005),hep-ph/0309332.[44] P. Aurenche, F. Gelis, and H. Zaraket, JHEP , 043(2002), hep-ph/0204146.[45] M. Djordjevic, Phys. Rev. C74 , 064907 (2006), nucl-th/0603066.[46] C. W. De Jager, H. De Vries, and C. De Vries, Atom.Data Nucl. Data Tabl. , 479 (1974).[47] M. R. Whalley, D. Bourilkov, and R. C. Group (2005),hep-ph/0508110.[48] K. J. Eskola, V. J. Kolhinen, and C. A. Salgado, Eur.Phys. J. C9 , 61 (1999), hep-ph/9807297.[49] K. J. Eskola, H. Paukkunen, and C. A. Salgado, JHEP , 31 (1983).[52] T. Sjostrand, Nucl. Phys. B248 , 469 (1984).[53] H. L. Lai et al. (CTEQ), Eur. Phys. J.
C12 , 375 (2000),hep-ph/9903282.[54] A. Adare et al. (PHENIX), Phys. Rev. Lett. , 232301 (2008), arXiv:0801.4020.[55] G.-Y. Qin et al., Phys. Rev.
C76 , 064907 (2007),arXiv:0705.2575.[56] T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey, andY. Nara, Phys. Lett.
B636 , 299 (2006), nucl-th/0511046.[57] A. Adil, H.-J. Drescher, A. Dumitru, A. Hayashigaki,and Y. Nara, Phys. Rev.
C74 , 044905 (2006), nucl-th/0605012.[58] S. S. Adler et al. (PHENIX), Phys. Rev. Lett. , 012002 (2007), hep-ex/0609031.[59] A. Adare et al. (PHENIX) (2008), 0804.4168.[60] T. Isobe (PHENIX), J. Phys. G34 , S1015 (2007), nucl-ex/0701040.[61] P. Aurenche, R. Baier, A. Douiri, M. Fontannaz, andD. Schiff, Nucl. Phys.
B286 , 553 (1987).[62] P. Aurenche, R. Baier, M. Fontannaz, and D. Schiff,Nucl. Phys.