Multistage Monte-Carlo simulation of jet modification in a static medium

TL;DR

This paper introduces a multistage Monte Carlo simulation framework for jet modification in a static medium, combining different models to capture the full evolution of jets over various scales, and compares their effects on jet quenching observables.

Contribution

It develops a novel multistage modeling approach where jets transition between different Monte Carlo generators, enhancing the simulation of jet evolution across multiple epochs.

Findings

Sequential generator transitions affect jet observables.

Varying generator boundaries impacts jet quenching results.

Proposed criterion improves modeling consistency.

Abstract

The modification of hard jets in an extended static medium held at a fixed temperature is studied using three different Monte-Carlo event generators (LBT, MATTER, MARTINI). Each event generator contains a different set of assumptions regarding the energy and virtuality of the partons within a jet versus the energy scale of the medium, and hence, applies to a different epoch in the space-time history of the jet evolution. For the first time, modeling is developed where a jet may sequentially transition from one generator to the next, on a parton-by-parton level, providing a detailed simulation of the space-time evolution of medium modified jets over a much broader dynamic range than has been attempted previously in a single calculation. Comparisons are carried out for different observables sensitive to jet quenching, including the parton fragmentation function and the azimuthal distribution of jet energy around the jet axis. The effect of varying the boundary between different generators is studied and a theoretically motivated criterion for the location of this boundary is proposed. The importance of such an approach with coupled generators to the modeling of jet quenching is discussed.

Equations

$$\hat{q}=\frac{⟨p_{\perp }^2{⟩}_L}{L}.$$

$$\Delta (t,t_0)$$

$$\frac{d{P}_q(p)}{dt}$$

$$+2P_g(p+k)\frac{d{\Gamma }_{q\bar{q}}^g(p+k,k)}{dkdt},$$

$$\frac{d{P}_g(p)}{dt}$$

$$+P_g(p+k)\frac{d{\Gamma }_{gg}^g(p+k,p)}{dkdt}$$

$$-P_g(p)\left(\frac{d{\Gamma }_{q\bar{q}}^g(p,k)}{dkdt}+\frac{d{\Gamma }_{gg}^g(p,k)}{dkdt}\theta (2k-p)\right).$$

$$\sqrt{⟨p_{\perp }^2⟩}=(\hat{q}p{)}^{1/4},$$

$$p_1\cdot \partial f_1(x_1,p_1)=E_1({\mathcal{C}}_{\mathrm{el}}+{\mathcal{C}}_{\mathrm{inel}}),$$

$${\Gamma }_{\mathrm{el}}$$

$$\times f_2({\vec{p}}_2)\left[1\pm f_3({\vec{p}}_3)\right]\left[1\pm f_4({\vec{p}}_4)\right]S_2(s,t,u)$$

$$\times (2\pi {)}^4{\delta }^{(4)}(p_1+p_2-p_3-p_4)|{\mathcal{M}}_{12\to 34}{|}^2,$$

$$⟨N_g⟩(E,T,t,\Delta t)=\Delta t\int dxd{k}_{\perp }^2\frac{d{N}_g}{dxd{k}_{\perp }^{2}dt},$$

$$\frac{d{N}_g}{dxd{k}_{\perp }^{2}dt}=\frac{2{\alpha }_s{C}_A\hat{q}P(x)k_{\perp }^4}{\pi {\left(k_{\perp }^2+x^2{m}^2\right)}^4}{\sin }^2\left(\frac{t-t_i}{2{\tau }_f}\right),$$

$$P(n)=\frac{⟨N_g{⟩}^n}{n!}{e}^{-⟨N_g{⟩}^n}.$$

$$Q_0^2=\sqrt{2E\hat{q}},$$

$$\hat{q}=C_R{\alpha }_s{\mu }_D^{2}T\log \left(\frac{6ET}{{\mu }_D^2}\right).$$

Full text

I Introduction

The extensive data sets and highly precise measurements that are now available for heavy-ion collisions at RHIC and the LHC, combined with theoretical developments, enable an increasingly quantitative understanding of the properties of the Quark-Gluon Plasma (QGP) Heinz et al. (2015); Akiba et al. (2015). The bulk behavior (related to low transverse momentum particles) of the plasma is described with increasing precision using relativistic viscous fluid dynamics Romatschke and Romatschke (2007); Song et al. (2011); Gale et al. (2013); Heinz and Snellings (2013); however, a detailed quantitative picture of the microscopic dynamics that is the cause of the evidently hydrodynamic behavior remains elusive.

Jet modification Gyulassy and Wang (1994); Wang and Gyulassy (1992); Baier et al. (1998, 1997a, 1997b); Zakharov (1997, 1996); Gyulassy et al. (2003); Kovner and Wiedemann (2003); Casalderrey-Solana and Salgado (2007); Majumder and Van Leeuwen (2011), and its associated transport coefficients Baier (2003); Majumder (2009a); Burke et al. (2014) provide multi-scale hard probes of the QGP. Jets start with an off-shellness or virtuality $Q^{2}$ that is typically of the order of the hard scale, orders of magnitude higher than any scale in the medium. As the partons in a jet split and radiate more partons, they lose virtuality faster than they lose energy, and may spend a portion of their path in the medium at a scale comparable to the medium scale $Q^{2}\sim\sqrt{\hat{q}E}$. Here, $E$ is the energy of the parton interacting with the medium, and $\hat{q}$ is the quark/gluon transport coefficient that denotes the broadening of the transverse (to the direction of the parton) momentum distribution, per unit length:

[TABLE]

In the equation above, $\langle\cdots\rangle_{L}$ denotes event average over a length $L$. At even lower (or near thermal) momentum scales the jets are expected to be strongly coupled with the medium.

The interaction of the jet with the medium and its splitting process differ at these various scales. To date, no calculation or simulation of jet modification in medium has accounted in a unified way for these different eras; all previous efforts have used a single formalism and applied it to the entire space-time history of the jet Burke et al. (2014). In this work, we provide the first treatment that accounts for different and complementary theoretical approaches, applied in succession to the space-time history of the jet in the medium.

The energy of a reconstructed jet can be theoretically traced back to a single hard parton that undergoes subsequent splits and turns into a parton shower. The energy and virtuality of each parton determines its interaction with the medium. Hence, accurate event generators should allow for each parton to be ascribed the interaction formalism valid for that particular energy and virtuality at any given time. For instance, in the domain of high energy and high virtuality $Q^{2}\gg\sqrt{\hat{q}E}$, one expects to apply the medium modified Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution Majumder (2009b); Deng and Wang (2010) based on the higher-twist formalism Wang and Guo (2001); Majumder (2012). As the virtuality approaches the medium induced scale, one expects to simulate a rate equation Jeon and Moore (2005); Turbide et al. (2005); Qin et al. (2009a) based on the BDMPS/AMY Baier et al. (1997a, b); Arnold et al. (2001, 2002) formalisms or a derivative of the higher twist formalism Wang and Zhu (2013); Cao et al. (2016, 2017); Chen et al. (2017).

To carry out such a program requires detailed analysis of the effect of merging formalisms at the parton level. A full simulation of jet modification in realistic heavy-ion collisions requires several different components: an initial state calculation that samples hard partons from the parton distribution functions of the incoming nuclei and generates hard (high $Q^{2}$) collisions, leading to the formation of hard final state partons that interact with the medium, a viscous fluid dynamical simulation followed by hadronization of the soft sector, the modification of the jet as it propagates through this evolving medium, and hadronization of the jet partons.

In this paper by the JETSCAPE Collaboration JET (a), we perform the first such calculation of jet quenching, using coupled event generators that are applicable to different eras of the in-medium evolution of jet. We take a simplified approach in this study and consider the case of a single mono-energetic parton generated at the surface of a uniform, static medium of fixed length $L$ in the direction of propagation, with all chemical potentials set to zero and fixed temperature $T$. The medium is assumed to have an infinite transverse extent (a brick) Armesto et al. (2012). This simplification of the problem has the advantage of eliminating possibly confounding effects from the dynamical evolution of the medium in a realistic heavy-ion collision, thereby focussing our attention on essential differences between the different energy loss formalisms applied to different jet evolution stages. There may also be effects that manifest themselves clearly in a static medium, e.g. the wake of a jet Neufeld and Muller (2009); Qin et al. (2009b), but are blurred by a dynamical medium that undergoes phase transitions He et al. (2015). Comparison of simulations in the brick and in a dynamical medium will disentangle these different aspects of jet modification.

In Ref. Armesto et al. (2012), the TECHQM tec and the JET collaboration JET (b) studied four different energy loss formalisms within the brick setup. In that effort, calculations were carried out in the event-averaged formalism of energy loss of the leading parton, and all formalisms were treated separately. The goal of the current effort by the JETSCAPE Collaboration is a first comparison between different event generators for a full parton shower evolution within a brick, and to explore the effect of combining different generators at the partonic level. These results will serve as benchmarks for future event generators to be developed within the JETSCAPE framework.

This paper is organized as follows: In Sec. II, we provide a brief overview to the various simulation formalisms and describe how they are coupled together. In Sec. III, we present the results of comparisons between the different approaches when used separately, along with those of the coupled approach. Concluding discussions and an outlook are provided in Sec. IV.

II Methodology

In this section, we first briefly summarize the three approaches for medium-modified parton showers to be employed here – MATTER, MARTINI and LBT – and then discuss how to combine the virtuality-ordered scheme for high virtuality partons and the time-ordered scheme for low virtuality partons into a unified theoretical approach. While MATTER represents the high virtuality scheme, MARTINI and LBT will be used for the low virtuality phase.

It should be pointed out that the single gluon emission kernel, which is sampled by MARTINI and LBT, is quite different in terms of the physics assumptions used to derive it: MARTINI is based on the BDMPS/AMY formalism and LBT is based on the higher-twist formalism. Also, the LBT generator includes the dynamics of soft partons in the medium, an aspect missing in MARTINI (and MATTER) so far. Briefly stated, energy lost by hard partons in LBT reappears as an enhancement of the soft spectrum of partons, while energy lost by hard partons in MATTER and MARTINI is considered to be irrecoverably lost to the medium. As one would expect, the presence of a multitude of such soft particles leads to noticeable difference both in the soft spectrum and its angular distribution away from the jet axis.

II.1 MATTER

The Modular All Twist Transverse-scattering Elastic-drag and Radiation (MATTER) event generator simulates the splitting of high energy high virtuality jets, i.e. jets whose virtuality $Q^{2}\gg\sqrt{\hat{q}E}$ where $E$ is the energy of the parton. At these high virtualities, the dominant mechanism of splitting is described by a medium-modified virtuality-ordered shower Majumder (2012, 2009b); Wang and Guo (2001); Majumder and Shen (2012). The underlying physical picture at this stage is that scattering in the medium produces a small variation in the vacuum shower process. The setup of the formalism ensures that the number of splittings dominates over the number of scatterings.

A virtuality-ordered shower is initiated by a single hard parton produced at a point $r$, whose forward light-cone momentum $p^{+}=(p^{0}+\hat{n}\cdot\vec{p})/\sqrt{2}$ has been specified ($\hat{n}=\vec{p}/|\vec{p}|$ represents the direction of the jet). One then samples a Sudakov form factor to determine its virtuality $t=Q^{2}$ Majumder (2013); Majumder and Putschke (2016),

[TABLE]

In the equation above, $\Phi$ represents a sum over phase factors that depend on $\zeta^{+}$, $p^{+}$, and $Q$. The transport coefficient $\hat{q}$ is evaluated at the location of scattering $\vec{r}+\hat{n}\zeta^{+}$. The function $P(z)$ is the vacuum splitting function. The maximum length sampled $\zeta^{+}_{\mathrm{MAX}}$ corresponds to $1.3\tau_{f}^{+}$, where $\tau_{f}^{+}$ is the mean light-cone formation time $\tau^{+}_{f}=2p^{+}/Q^{2}$ Kordell and Majumder (2017).

Once $Q^{2}$ is determined, $z$ can be determined by sampling the splitting function $P(z)$. The transverse momentum of the produced pair (transverse to $\hat{n}$) is fixed once the shower is determined, by inspecting the difference in invariant mass between the parent and siblings of a given split. To this one may add the transverse momentum generated by the propagation through the medium. This process is continued until the $Q^{2}$ reaches a predetermined value of $Q_{0}^{2}$. The parton splitting process stops at this point. The final partons at this stage may then be passed to a hadronization routine, or to another formalism.

II.2 MARTINI

MARTINI (Modular Algorithm for Relativistic Treatment of Heavy IoN Interactions) is a Monte Carlo event generator for simulating jets in heavy ion collisions Schenke et al. (2009a). The nucleon-nucleon collisions and the vacuum shower are generated by Pythia 8. Then MARTINI deals with the parton evolution in a QGP medium according to the AMY formalism for the radiative energy loss rates Arnold et al. (2002, 2003) combined with collisional processes Schenke et al. (2009b). Those energy loss rates depend on the thermal background, whose information is provided by hydrodynamic calculations such as MUSIC Schenke et al. (2010). Finally, hadronization of the evolved partons can be performed by Pythia based on the Lund string model Andersson et al. (1983). In this work, the initial parton shower is generated by MATTER, a static brick is adopted to study medium modification and all final spectra are analyzed at the parton level.

In MARTINI, the time evolution of the jet momentum distribution is governed by a set of coupled rate equations, which take the following forms:

[TABLE]

$d\Gamma^{a}_{bc}(p,k)/dkdt$ is the transition rate for a process where a parton $a$ of energy $p$ emits a parton $c$ of energy $k$ and becomes a parton $b$. The factor of 2 in front of $d\Gamma^{g}_{q\bar{q}}$ takes into account the fact that $q$ and $\bar{q}$ are distinguishable. For the $g\rightarrow gg$ process, the $\theta$ function is there to avoid double counting of final states. Here $P_{q}(p)$ and $P_{g}(p)$ are the energy distribution of quarks and gluons, respectively. The integration range with $k<0$ indicates energy gain from the thermal medium; the range with $k>p$ for $q\rightarrow qg$ process represents annihilation against anti-quark of energy $k-p$ from the medium.

The AMY formalism describes energy loss of hard jets in heavy ion collisions as parton bremsstrahlung in the evolving QGP medium. The effective kinetic theory described in Arnold et al. (2003) assumes that quarks and gluons in the medium are well defined (hard) quasi-particles and have typical momentum of the order of temperature $T$ and thermal mass of order $gT$. Under this assumption, the radiation rate can be calculated by means of integral equations Arnold et al. (2002).

In the current version of MARTINI, the radiative energy loss mechanism is improved by implementing the effects of finite formation time and running coupling. The formation time of the radiation process increases with $\sqrt{p}$ and a hard parton and an emitted parton are coherent within that time. This interference effect suppresses the radiation rate at early times after the original radiation.

For the renormalization scale of running coupling constant $\alpha_{s}(\mu)$, we use the root mean square of the momentum transfer $\sqrt{\langle p^{2}_{\perp}\rangle}$ between the two particles, parameterized as

[TABLE]

where $\hat{q}$ is the averaged momentum transfer squared per scattering and $p$ the energy of the mother parton Young et al. (2013).

II.3 LBT

A Linear Boltzmann Transport (LBT) model has been devised to describe the in-medium parton showers at low virtuality scale Li et al. (2011); Wang and Zhu (2013); He et al. (2015); Cao et al. (2016, 2017); Chen et al. (2017). In the absence of a mean field, the evolution of the phase space distribution of a hard parton “1” with $p_{1}^{\mu}=(E_{1},\vec{p}_{1})$ can be described using the Boltzmann equation

[TABLE]

in which $\mathcal{C}_{\mathrm{el}}$ and $\mathcal{C}_{\mathrm{inel}}$ are collision integrals for elastic and inelastic scatterings.

For elastic scattering, the collision term $\mathcal{C}_{\mathrm{el}}$ is evaluated with the leading-order matrix elements for all possible “$12\rightarrow 34$” processes between the given jet parton “1” and a thermal parton “2” present in the medium background. To regulate the collinear ($u,t\rightarrow 0$) divergence of the matrix element, $S_{2}(s,t,u)=\theta(s\geq 2\mu_{\mathrm{D}}^{2})\theta(-s+\mu_{\mathrm{D}}^{2}\leq t\leq-\mu_{\mathrm{D}}^{2})$ is imposed in which $\mu_{\mathrm{D}}^{2}=g^{2}T^{2}(N_{c}+N_{f}/2)/3$ is the Debye screening mass. The elastic scattering rate of parton “1” can be evaluated as follows:

[TABLE]

in which $\gamma_{2}$ represents the spin-color degeneracy of parton “2”. And therefore, the probability of elastic scattering of parton “1” in each time step $\Delta t$ is $P_{\mathrm{el}}=\Gamma_{\mathrm{el}}\Delta t$.

For inelastic scattering, or medium-induced gluon radiation, the average number of emitted gluons from a hard parton in each time step $\Delta t$ is evaluated as Cao et al. (2013, 2015, 2016)

[TABLE]

in which the differential spectrum of radiated gluon is taken from the higher-twist energy loss formalism Guo and Wang (2000); Majumder (2012); Zhang et al. (2004):

[TABLE]

where $x$ and $k_{\perp}$ are the fractional energy and transverse momentum of the emitted gluon with respect to its parent parton, $\alpha_{s}$ is the strong coupling constant, $C_{A}=N_{c}$ is the gluon color factor, $P(x)$ is the splitting function, and $\hat{q}$ is the transport coefficient due to elastic scattering and can be obtained by evaluating Eq. (II.3) weighted by the transverse momentum broadening of parton “1”. The mass dependence of gluon emission from heavy quark is included in Eq. (8). In addition, $t_{i}$ denotes an “initial time” or the production time of the parent parton from which the gluon is emitted, and $\tau_{f}={2Ex(1-x)}/{(k_{\perp}^{2}+x^{2}m^{2})}$ is the formation time of the radiated gluon. To avoid possible divergence as $x\rightarrow 0$, a lower cut-off $x_{\mathrm{min}}=\mu_{D}/E$ is implemented for the energy of emitted gluon. Multiple gluon radiation is allowed in each time step. Different emitted gluons are assumed independent with each other and therefore their number $n$ obeys a Poisson distribution with the mean as $\langle N_{g}\rangle$:

[TABLE]

Thus, the probability for the total inelastic scattering process is $P_{\mathrm{inel}}=1-e^{-\langle N_{g}\rangle}$. Note that for the $g\rightarrow gg$ process, $\langle N_{g}\rangle/2$ is taken as the mean instead to avoid double counting.

To combine the above elastic and inelastic processes, the total scattering probability is divided into two regions: pure elastic scattering with probability $P_{\mathrm{el}}(1-P_{\mathrm{inel}})$ and inelastic scattering with probability $P_{\mathrm{inel}}$. Thus the total scattering probability is $P_{\mathrm{tot}}=P_{\mathrm{el}}+P_{\mathrm{inel}}-P_{\mathrm{el}}\cdot P_{\mathrm{inel}}$. Based on these probabilities, the Monte Carlo method is applied to determine whether a given jet parton is scattered inside the thermal medium and whether the scattering is pure elastic or inelastic. With a selected scattering channel, the energies and momenta of the outgoing partons are then sampled based on the corresponding differential spectra given by Eq. (II.3) and (8). The only parameter in this LBT model is the strong coupling constant $\alpha_{s}$ that quantifies the jet-medium interaction, which is determined by comparing model calculation to experimental data of single heavy and light flavor hadron production, single inclusive jet production and $\gamma$-jet production in heavy-ion collisions Cao et al. (2017, 2016); Li et al. (2011); Wang and Zhu (2013); Chen et al. (2017).

II.4 Combining MATTER and MARTINI/LBT

To establish a unified framework for parton showers, we apply MATTER to partons with large virtuality and MARTINI or LBT to partons with low virtuality. MARTINI and LBT, based on AMY and higher-twist energy loss formalisms respectively, will be compared to each other and shown to be consistent within the kinematic region we investigate. To begin with, each hard parton directly produced in hard scatterings is placed in MATTER, in which its virtuality-ordered splitting process is simulated as described in Sec. II.1. In each splitting, the virtualities of the daughter partons are much smaller than that of the parent. If the virtuality of a given parton in the shower drops below a certain scale $Q_{0}$, it is passed to MARTINI/LBT for the subsequent time-ordered evolution inside the medium.

One crucial quantity in this unified framework is the separation scale $Q_{0}$ between MATTER and MARTINI/LBT evolutions. Two different schemes, fixed $Q_{0}$ and dynamical $Q_{0}$, are applied and compared in this study. For the former, fixed values of $Q_{0}$ (1, 2 or 3 GeV) are used as the scale below which both vacuum and medium-modified showers in MATTER cease. For the latter, we define $Q_{0}^{2}=\hat{q}\tau_{f}$ for each parton. This quantifies the average virtuality gain of each parton from its scattering with the medium and serves as a reasonable separation scale between the virtuality-loss splitting process in MATTER and the near-constant-virtuality scattering process in MARTINI/LBT. With $\tau_{f}=2E/Q_{0}^{2}$, we have

[TABLE]

in which the quark/gluon transport coefficient is obtained from a finite temperature field theory calculation as Gyulassy and Wang (1994); Caron-Huot and Gale (2010),

[TABLE]

For small $E$, $Q_{0}=1$ GeV is used if the estimate above yields a $Q_{0}$ smaller than 1 GeV. Thus $Q_{0}$ is dynamically determined based on the energy of each parton and the local temperature of its surrounding medium.

Note that for this dynamical scenario, Eqs. (10,11) are only well defined in a thermal medium, and thus fixed $Q_{0}=1$ GeV is used when partons travel outside the color deconfined nuclear matter (or in vacuum). The value $Q_{0}=1$ GeV is from the standard practice to regard 1 GeV to be the lowest scale where pQCD is expected to be valid. Partons with a virtuality below this scale are assumed to be strongly coupled with the medium, or undergo hadronization if outside the QGP. Since neither of these effects are incorporated in this first attempt at a multi-stage event generator, the $Q_{0}$ is held at a minimum of 1 GeV until exit from the brick.

In this work, we start with a single quark and let it evolve either in vacuum or through a static brick with fixed temperature 250 MeV. With a virtuality greater than $Q_{0}$, each parton in the shower evolves in MATTER; but after hitting $Q_{0}$, it starts evolving in either MARTINI or LBT. We will investigate how such combined theoretical approach affects the energy distribution and jet broadening of the finally produced partons, compared to traditional methods in the literature where a single energy loss approach is applied. The initial energy of the single quark and the length of the brick medium will be varied so that one may study in which region the virtuality ordered splitting in the parton shower is more prevalent and in which region the time ordered splitting dominates. Note that it is possible that a highly virtual parton still has not hit $Q_{0}$ after it traverses the entire brick medium. If so, the vacuum shower is attached in MATTER after its in-medium evolution until each daughter parton reaches $Q_{0}$; and MARTINI/LBT does not have any effect on the shower in such case. Therefore, all plots describe the distribution of partons at the exit from the brick or the moment hitting $Q_{0}$, whichever comes later. Throughout our calculation, a fixed value of $\alpha_{s}=0.3$ is used in the MARTINI/LBT portion and a running $\alpha_{s}$ used in the higher virtuality MATTER portion.

Before discussing the energy distribution and jet broadening of the final partons, we first investigate the time $t_{0}$ at which a parton in the shower hits $Q_{0}$ as a function of its energy in Fig. 1. This $t_{0}$ is obtained by summing the formation times of all previous splittings in MATTER before the produced parton hit $Q_{0}$ as discussed in Sec. II.1. In Fig. 1, the left panel corresponds to 50 GeV energy for the initial quark, and the right panel corresponds to 200 GeV. In each panel, 3 different fixed $Q_{0}$ and the dynamical $Q_{0}$ scenarios are compared. The two main figures show results in a brick with length of 4 fm, and the two subfigures inside show the ratio between 4 fm and 8 fm.

From these figures, one may observe the time it takes to evolve an energetic parton down to $Q_{0}$ can be long (compared to the formation time of the QGP $\tau_{0}\sim 0.6$ fm in realistic heavy-ion collisions). The switching time $t_{0}$ increases if the separation scale $Q_{0}$ decreases or the initial energy $E_{\mathrm{init}}$ (i.e. the possible maximum virtuality) of the parton increases. For $E_{\mathrm{init}}=50$ GeV, $t_{0}$ for the dynamical $Q_{0}$ is consistent with that for fixed $Q_{0}=2$ GeV at the high energy end of the final parton spectrum. This is a natural result from Eq. (10). And if $E_{\mathrm{init}}$ is increased to 200 GeV, $t_{0}$ for the dynamical $Q_{0}$ is then consistent with that for fixed $Q_{0}=3$ GeV at the high energy end. On the other hand, for final partons with lower energies, $t_{0}$ for the dynamical $Q_{0}$ approaches that for fixed $Q_{0}=1$ GeV since at low $E$ and after partons travel outside the brick (back into vacuum) $Q_{0}$ is set as 1 GeV.

Changing the size of the brick from $L=4$ to 8 fm also affects the value of $t_{0}$. For the scenarios of fixed $Q_{0}$, extending the length of the brick increases scattering of partons inside the medium. This is a virtuality gain process and thus may delay the time $t_{0}$ for each parton to hit $Q_{0}$. For $E_{\mathrm{init}}=50$ GeV, this only affects the $Q_{0}=1$ GeV scenario, since for $Q_{0}=2$ and 3 GeV, most partons hit $Q_{0}$ before 4 fm and thus adding another 4 fm of length/time has slight effect. On the contrary, for $E_{\mathrm{init}}=200$ GeV, extending $L$ from 4 to 8 fm clearly increases $t_{0}$ for both fixed $Q_{0}=1$ and 2 GeV. The opposite effect of varying the brick size is observed for the scenario of dynamical $Q_{0}$. Unlike the fixed $Q_{0}$ scenario where the same $Q_{0}$ is adopted for both in-medium and vacuum parton shower, the dynamical scenario uses $Q_{0}$ from Eq. (10) (usually larger than 1 GeV) inside the brick but 1 GeV outside. And therefore, extending $L$ from 4 to 8 fm increases the range where larger $Q_{0}$ is applied and thus shortens $t_{0}$. These discussions on the separation time $t_{0}$ between MATTER and MARTINI/LBT evolutions will be helpful in understanding the final parton spectra within our unified theoretical approach as will be shown in the next section.

III Energy distribution and jet broadening in MATTER+MARTINI/LBT

In this section, we will present the energy distribution and jet angular distribution of the parton shower and compare our unified approach to traditional approaches in which a single energy loss scheme is applied. In the simplified scheme of the brick, these two distributions form the underlying basis of almost all jet observables. In the limit of near onshell particles, only three components of the four-momentum are relevant, which we have chosen as the energy $E_{k}$, and the components transverse to the jet axis, $k_{\perp}$ [$k=\left(E_{k}^{2}-Q_{0}^{2}-k_{\perp}^{2}\right)$, where $Q_{0}=1$ GeV].

In Fig. 2, we show the energy distribution of final partons with the brick size $L=4$ fm. The parton energy of 50 GeV is used for the initial quark in Fig. 2 and 200 GeV for Fig. 2. In each figure, the four panels correspond to different choices of the separation scale $Q_{0}$: three fixed cases and one dynamical case. For the three panels with fixed $Q_{0}$, we compare results between vacuum shower, medium modified shower through MATTER alone, MARTINI alone and MATTER+MARTINI. For the panel with dynamical $Q_{0}$, we compare results between medium modified shower through MATTER alone, MATTER+MARTINI and MATTER+LBT. Note that for the vacuum shower or the medium modified shower in MATTER alone, we let all partons evolve down to the given $Q_{0}$ in the virtuality ordered scheme. For the pure MARTINI evolution, we let each parton produced by the vacuum shower evolve through MARTINI for the entire 4 fm as the traditional implementation in the literature. MATTER+MARTINI/LBT corresponds to our unified approach as discussed in Sec. II.4 in which each parton evolves though MATTER till $t_{0}$ (i.e., with virtuality larger than $Q_{0}$) after which it evolves through MARTINI or LBT until 4 fm.

From Fig. 2, we observe the energy spectra of final partons from MARTINI and LBT are consistent with each other, within our setup, except for the very low energy region, since in this calculation MARTINI is set up to regard emitted partons with an energy less than 2 GeV as a part of the medium. The similarity between MARTINI and LBT may be due to the fact that the results shown in Fig. 2 are consequence of multiple emissions. Even though the single gluon emission spectra differ in MARTINI and LBT in detail, the multiple emission spectra obtained either by solving Eq. (II.2) or via the Poisson ansatz Eq. (9) may wash away some of the differences.

For $E_{\mathrm{init}}=50$ GeV, the maximum scale of the medium $\sqrt{\hat{q}\tau_{f}}$ probed by the hard parton is around 2 GeV, and therefore pure MATTER evolution (compared to the vacuum shower) leads to a minimal suppression of the spectra if one sets $Q_{0}\gtrsim 2$ GeV. Within our unified approach of MATTER+MARTINI, we see that the effect of just applying MARTINI for the entire length of the medium on a vacuum shower leads to a considerably larger suppression of the spectrum of final state partons, specifically for high energy partons, compared to the suppression of first allowing the partons to split using MATTER and then with MARTINI once the virtuality reaches $Q_{0}$. The cause of this, clearly illustrated in Fig. 1, is that it takes a longer time for more energetic partons to evolve down to $Q_{0}$ and therefore leaves shorter time for them to evolve inside MARTINI. This is more apparent when $Q_{0}$ is smaller. In the end, results for the dynamical $Q_{0}$ scenario is close to the fixed $Q_{0}=2$ GeV scenario which can be naturally understood from Eq. (10).

If the energy of the initial quark is increased by a factor of 4, one may observe in Fig. 2 that the amount of suppression obtained by pure MATTER evolution (vs. vacuum) is considerable up to a higher $Q_{0}$ cut (around 3 GeV), since the scale of the medium probed by the jet increases with the jet energy. In addition, the additional amount of suppression obtained from the MARTINI portion in our MATTER+MARTINI approach is reduced compared to the $E_{\mathrm{init}}=50$ GeV scenario, since the switching time $t_{0}$ between MATTER and MARTINI is larger (see Fig. 1). For $E_{\mathrm{init}}=200$ GeV, results with dynamical $Q_{0}$ are closer in form to those with the fixed $Q_{0}=3$ GeV.

In Fig. 3 we investigate the path length dependence of the energy distribution of the final partons within our MATTER+LBT framework. In Fig. 3, the energy of the initial quark is set as 50 GeV. We observe for the case of pure MATTER evolution with $Q_{0}=1$ GeV, the energy loss of the incoming parton is greater, i.e., the spectra is more suppressed, when the medium size $L$ is extended from 4 to 8 fm. However, for $Q_{0}\gtrsim 2$ GeV, there is no apparent difference between $L=4$ and 8 fm since most partons evolve down to the virtuality of $Q_{0}$ before 4 fm (see Fig. 1) and these partons are propagated without effect through the remaining 4 fm of medium. On the other hand, extending $L$ from 4 to 8 fm leaves longer time for parton evolution through LBT. In other words, the LBT (and MARTINI) evolution has a much larger effect on the suppression if the path length inside the medium is longer, as is readily observed for $Q_{0}\gtrsim 2$ GeV. For $Q_{0}=1$ GeV, the separation point between the pure MATTER curve and the MATTER+LBT curve also shifts to the right (from around 10 to around 20 GeV) when $L$ is extended from 4 to 8 fm, indicating a wider range of partons affected by the LBT evolution when $L$ is larger.

A similar investigation of the length dependence is presented in Fig. 3 where the energy of the initial quark is increased to 200 GeV. Compared to Fig. 3, the difference between $L=4$ and 8 fm for the pure MATTER scenario is not only obvious for $Q_{0}=1$ GeV, but also for $Q_{0}=2$ GeV now, because with $E_{\mathrm{init}}=200$ GeV, it takes much longer than 4 fm for MATTER to evolve partons down to 2 GeV (see Fig. 1) and thus evolution between 4 and 8 fm becomes important. However, since now it takes longer for MATTER to evolve partons down to a given $Q_{0}$, shorter time is left for the subsequent LBT evolution and thus extending $L$ from 4 to 8 fm results in a lower effect of the LBT contribution to the suppression of the final parton spectra compared to the previous results for $E_{\mathrm{init}}=50$ GeV.

Finally, we also investigate the energy flow inside the jet within our unified approach of MATTER+MARTINI/LBT. In Fig. 4, we present the energy distribution with respect to the jet cone angle, compared between vacuum shower, in-medium shower from pure MATTER, pure LBT and MATTER+LBT/MARTINI. The four panels are for the four different choices of $Q_{0}$, and in order to clearly display the effect of medium modification, we show the ratios between the medium modified spectra and their corresponding vacuum spectra in the four subfigures (For the dynamical scenario, $Q_{0}=1$ GeV is used for the vacuum baseline in calculating the ratio).

In most cases, medium induced emission leads to a depletion of the energy at small angles and a milder enhancement at larger angles. However, there are several caveats in this effect. If the switching virtuality is set at $Q_{0}=1$ GeV, or MATTER evolution is allowed to run down to the minimum possible value, then there is a narrowing of the cone (enhancement) at very small angles ($\theta<0.05$). This is mostly caused by MATTER evolution. Both MARTINI and LBT lead to a suppression of the energy in the most collinear bins.

In MATTER, the narrowing is caused by the enhanced Sudakov factor, which leads to an enhanced probability of collinear emission. The broadening and enhancement at larger angles is mostly present in LBT. The particles shifted to larger angles tend to be low energy partons, which are explicitly excised from the MATTER and MARTINI showers. Also, LBT is more effective in shifting the energy distribution into larger angles than MATTER Majumder and Wang (2006), since elastic scattering is also included in LBT. The contributions from LBT and MARTINI are also compared: MARTINI is consistent with LBT at small angle, or for high-energy partons, but they deviate at large angles since partons with very low energy are regarded as part of the background medium in the current MARTINI simulation.

The results presented in Fig. 4 are somewhat surprising, specifically those in the top left panel, where $Q_{0}=1$ GeV. The multi-stage simulation involving MATTER and LBT qualitatively reproduces the feature seen in experimental data of an enhancement at small angles, a depletion at intermediate, and an enhancement at larger angles (see inset). This is a feature absent in any single generator applied to a static medium. It has been proposed that a back reaction from a dynamical medium, enhanced by radial flow is responsible for the enhancement at large angles Casalderrey-Solana et al. (2017). The simulations presented in this work, do not involve a dynamical medium, however, they do involve elastic scattering of soft partons off the medium. Due to the switch between MATTER and LBT at $Q_{0}=1$ GeV, hard partons at the lowest angles spend most of their lifetime in the MATTER phase and are not affected by the rescattering. Whereas, soft wide angle partons that reach a virtuality of $Q_{0}$ deep in the medium, are scattered out to larger angles, this leads to the depletion at intermediate angles and the enhancement at larger angles. As $Q_{0}$ is increased, even the hard partons reach the transition between MATTER and LBT some distance from exit and are affected by the re-scattering in LBT. As a result, the enhancement at small angles is lost.

IV Summary

In this work, different event generators for medium modified parton showers are coupled for the first time at the parton level. The virtuality-ordered event generator MATTER, based on the higher twist energy loss formalism, is adopted for shower partons at large virtuality, while two time-ordered transport models, MARTINI and LBT, are applied at low virtuality. MARTINI, based on the AMY energy loss formalism, and LBT, based on the higher twist energy loss formalism, are shown to be consistent with each other. Both fixed and dynamical separation scales between MATTER and MARTINI/LBT have been explored in this study. Varying the parton energy, medium length and virtuality separation scale $Q_{0}$, we studied the changing relative weights of the different energy loss schemes on the medium modifications experienced by the jet.

It was shown that the time $t_{0}$ it takes to evolve each parton down to a given $Q_{0}$ can be large. The switching time $t_{0}$ grows with increasing initial and final parton energies and decreases with larger $Q_{0}$. Its dependence on the size of the thermal medium varies with the parton energy and $Q_{0}$. When combining these two schemes – MATTER and MARTINI/LBT – into a unified approach, we observe that MATTER plays a larger role in the evolution of high energy partons whereas MARTINI/LBT play a larger role for lower energy ones since the latter reach $Q_{0}$ earlier and thus are subject to MARTINI/LBT evolution for a larger fraction of time. A larger value of $Q_{0}$ typically suppresses the medium modification by MATTER and thus increases the relative contribution from MARTINI/LBT. The importance of the MARTINI/LBT contribution also increases with the path length in the medium, but this increase may be suppressed if the initial parton starts with a higher energy. Finally, it was shown that jet can be broadened more effectively, i.e. stronger energy flow towards larger jet cone angles, if elastic scattering is included in addition to parton splitting. And combining MATTER and LBT may provide the non-monotonic nuclear modification of the angular distribution of jet energy within a small jet cone – a feature that is hard to obtain when a single energy loss mechanism is applied.

This work contributes a crucial step towards establishing a complete theoretical picture of parton shower evolution in relativistic heavy-ion collisions, applying within a single unified framework different energy loss formalisms based on complementary approximations, consistently to different kinematic regions during the parton evolution history. Our study will be extended in the near future in two directions. The framework will be coupled to a hydrodynamic background in order to study jet evolution in a realistically evolving dynamical medium. Secondly, additional theoretical schemes for parton evolution will be added to this framework, such as the AdS/CFT formalism for almost thermalized partons. All these insights will be incorporated in the upcoming JETSCAPE framework.

Acknowledgments

This work was supported in part by the National Science Foundation (NSF) within the framework of the JETSCAPE collaboration, under grant numbers ACI-1550172 (R.A.B. and D.V.), ACI-1550221(S.R. and R.J.F.), ACI-1550223 (D.B. and U.H.), ACI-1550225 (J.B. and J.C.), ACI-1550228 (L.-G.P. and X.-N.W.), and ACI-1550300 (S.C., A.K., D.P., A.M., C.Si. and R.A.S.). It was also supported in part by the NSF under grant numbers PHY-1207918 (M.K.), 1516590 (S.R. and R.J.F.), and by the US Department of Energy, Office of Science, Office of Nuclear Physics under grant numbers DE-AC02-05CH11231 (Y.H. and X.-N.W.), DE-AC52-07NA27344 (R.A.S.), DE-SC0012704 (B.S. and C.S.), DE-SC0013460 (S.C., A.K., A.M. and C.Si.), DE-SC0004286 (U.H.), and DE-FG02-05ER41367 (S.A.B.). The work was also supported in part by the National Science Foundation of China (NSFC) under grant number 11521064, Ministry of Science and Technology (MOST) of China under Projects number 2014CB845404 (Y.H., T.L. and X.-N.W), and in part by the Natural Sciences and Engineering Research Council of Canada (C.G., S.J. and C.P.) and the Fonds de recherche du Québec - Nature et technologies (FRQ-NT) (G.V.). Computations were made in part on the supercomputer Guillimin from McGill University, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), NanoQuébec, Réseau de Médicine Génétique Appliquée (RMGA) and FRQ-NT. C.G. gratefully acknowledges support from the Canada Council for the Arts through its Killam Research Fellowship program. C.S. gratefully acknowledges a Goldhaber Distinguished Fellowship from Brookhaven Science Associates.