Light and heavy flavor dijet production and dijet mass modification in heavy ion collisions

Back-to-back light and heavy flavor dijet measurements are promising experimental channels to accurately study the physics of jet production and propagation in a dense QCD medium. They can provide new insights into the path length, color charge, and mass dependence of quark and gluon energy loss in the quark-gluon plasma produced in reactions of ultra-relativistic nuclei. To this end, we perform a comprehensive study of both light and heavy flavor dijet production in heavy ion collisions. We propose the modification of dijet invariant mass distributions in such reactions as a novel observable that shows enhanced sensitivity to the QGP transport properties and heavy quark mass effects on in-medium parton showers. This is achieved through the addition of the jet quenching effects on the individual jets as opposed to their subtraction. The latter drives the subtle effects on more conventional observables, such as the dijet momentum imbalance shifts, which we also calculate here. Results are presented in Pb+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV for comparison to data at the Large Hadron Collider and in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV to guide the future sPHENIX program at the Relativistic Heavy Ion Collider.


I. INTRODUCTION
The high energy nuclear physics and particle physics communities are in the planning process for the upcoming proton and heavy ion runs at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). This is an opportune time to reflect on the success of recent runs and explore new opportunities. Important observables in high energy and heavy ion physics are related to hadronic jets. For example, the dominant Higgs decay channel H → bb was only recently observed [1,2] and involved b-jet reconstruction in the final state. There has also been a resurgence in the Quantum Chromodynamics (QCD) theory of jets with emphasis on their substructure, for a recent review of this physics in elementary p+p collisions see Ref. [3]. In collisions of ultra-relativistic nuclei at RHIC and the LHC, the modification of the production cross sections and substructure of jets is more sensitive to the in-medium strong interaction dynamics in comparison to the leading hadron attenuation [4]. As such, jets are excellent diagnostics of the hot and dense medium, the quark-gluon plasma (QGP), that is created in heavy ion collisions (HIC). These jet quenching phenomena have been widely studied at both RHIC and the LHC, for a recent review of jet physics in HIC see Ref. [5].
Heavy flavor physics in reactions of ultra-relativistic nuclei is another incredibly active area of research [6] that predates the study of jets. For the purposes of this paper, we will restrict our discussion to open heavy flavor, where experimental measurements and related phenomenology have traditionally focused on D-meson and B-meson production. There is a great deal of interest in the use of heavy flavor to constrain the transport properties of the QGP [7,8]. Still, the mechanisms of its in-medium modification are not yet fully understood. In light of this, heavy flavor jets have been proposed as a new tool to test the theory of heavy quark production, parton shower formation, and modification in nuclear matter. The first theoretical study of single inclusive b-jet production in HIC [9] has found that the cross section receives a large contribution from prompt gluons, where heavy flavor emerges from gluon splitting only in the late stages of the parton shower evolution. Thus, the suppression of inclusive b-jets at high transverse momenta can be nearly as large as the quenching of light jets, as confirmed by the first CMS measurement [10]. For the same reason, the connection between b-jets suppression and b-quark energy loss can be quite indirect. On the other hand, B-meson-tagged b-jets [11] suppress such a contribution from gluon splitting, and are most effective in ensuring that the dominant fraction of recoiling jets originate from prompt b-quarks. Such a conclusion also applies to the back-to-back b-tagged dijet production, as we will show below. New measurements of heavy flavor jets and their substructure, which is particularly sensitive to mass effects on in-medium parton shower evolution [12], are expected from the upcoming LHC runs and from the future sPHENIX experiment at RHIC [13].
Back-to-back jet pair (or dijet) production is among the most exciting channels used to probe QGP properties, where one typically focuses on the most energetic ("leading") and second most energetic ("subleading") jets. It is instructive to recall that the first definitive measurement of quenching effects on reconstructed jets came from the enhanced dijet asymmetry measurements at the LHC [14,15]. Further studies of this observable have been carried out not only at the LHC [16], but also at RHIC [17]. The origin of the additional imbalance to the dijet transverse momentum distribution in heavy ion collisions in comparison with the elementary p+p collisions has been attributed to the path length and color charge dependence of parton energy loss [18][19][20]. The interplay of Sudakov and in-medium collisional broadening on dijet acoplanarity has also been explored [21]. More recently, efforts have been put forward to understand the dependence of the quenching on the type of parton that initiates the jet. The first measurement of the back-to-back b-jet momentum imbalance [22] has been performed at the LHC and modeled theoretically [23].
The dijet asymmetry and momentum imbalance measure the difference of potentially large attenuation effects on the leading and subleading jets. Thus, those observables show a somewhat reduced sensitivity to the physics of jet quenching and the transport properties of the QGP. It has been pointed out early on that the asymmetry and momentum imbalance shifts in HIC may be influenced by background fluctuations [24] and, more recently, by parton shower fluctuations on an event-by-event basis [25]. To this end, we set out to find an observable where the effects that arise from the in-medium modification of parton showers add rather than subtract, and lead to enhanced sensitivity to the interactions of jets in the QGP, as well as the mass dependence of parton energy loss.
In the current work, we provide an extensive study of dijet production in heavy ion collisions at RHIC (or sPHENIX) kinematics and at LHC energies for both inclusive and b-tagged dijets. We compare the similarities and differences between those channels in A+A collisions to understand the flavor dependence of the quenching effects. Most importantly, we propose to use the dijet invariant mass modification as a novel probe of the QGP. An earlier study of dijet mass in proton-nucleus collisions showed negligible cold nuclear matter effects [26], suggesting that any significant modification of the dijet invariant mass distribution in A+A collisions arises from radiative and collisional energy loss processes in the QGP. At the same time, we include the studies for the more conventional observables such as two-dimensional nuclear modification factor R AA as a function of leading and subleading jet transverse momenta, and the imbalance z J distribution. We present theoretical predictions at √ s N N = 200 GeV for future Au+Au collisions relevant to the sPHENIX kinematics at RHIC and at √ s N N = 5.02 TeV for comparison to Pb+Pb data at the LHC.
The rest of our paper is organized as follows. In Sec. II, we present the evaluation of the differential cross sections for both inclusive and b-tagged dijet production in p+p collisions using the Pythia 8 event generator [27]. We also determine the flavor origin of the dijet production for the proper implementation of the energy loss effects. In Sec. III, we first present the basic formalism used to generate dijet invariant mass distributions and imbalance distributions, starting from the double differential cross section in terms of the transverse momenta of leading and subleading jets. We then provide the information on how we implement the medium effects to obtain the modification of inclusive and b-tagged dijet production in dense QCD matter. In Sec. IV, we present our phenomenological results for both RHIC and LHC kinematics. We give predictions for sPHENIX at RHIC, and provide detailed comparison with the most recent experimental measurements by the CMS collaboration at the LHC. We conclude our paper in Sec. V.

II. LIGHT AND HEAVY FLAVOR DIJET PRODUCTION IN P+P COLLISIONS
In this section, we present the evaluation of the double differential cross sections for inclusive and b-tagged dijet production in p+p collisions using Pythia 8 [27], which is a widely-used high energy phenomenology event generator that describes the main properties of the event structure well. In our simulations, 8 million events are simulated for each of these two processes. We construct jets with the anti-k T jet clustering algorithm [28], where a b-jet is identified if there is at least one b-quark within the jet.
Both inclusive and b-tagged dijet production in p+p collisions have been measured at the LHC. To show the validation of Pythia 8 simulation against experimental measurements on inclusive dijet production in p+p collisions, we present dijet cross section as a function of dijet invariant mass in the left panel of Fig. 1, compared to experimental measurements by the CMS [29] collaboration at center-of-mass energy √ s = 7 TeV at the LHC. Here, the dijet invariant mass m jj is defined as with p L and p S being the four-momenta for the leading and subleading jets, respectively. The jets are constructed with a jet radius R = 0.6, along with the following rapidity cut 0.5 < y * < 1.0, TeV. The left is for inclusive dijets from CMS collaboration [29], while the right is for b-tagged dijets from ATLAS collaboration [30].
where y * = |y L − y S | with y L (y S ) being the rapidity of leading and subleading jets. At the same time, we implement additional cuts on the transverse momentum and rapidity of individual jets, which are matched to those given in the experimental paper [29]. The red histograms are the results from Pythia 8 simulations. As one can see, the Pythia 8 event generator describes the experimental dijet invariant mass data very well. This gives us confidence in extracting information on the parton flavors initiating the dijets in heavy ion collisions.
In the right panel of Fig. 1, we compare our Pythia 8 simulation for b-tagged dijet invariant mass distribution with the ATLAS measurement [30] at √ s = 7 TeV. The jet radius is R = 0.4 and the distance in rapidity and azimuthal angle between b-quark and the b-jet, ∆R = (∆η) 2 + (∆φ) 2 , is required to be smaller than 0.3. Additionally, the transverse momentum of the b-quark is required to satisfy p T > 5 GeV. All other event selection and kinematic cuts are implemented to match the experimental measurements. For details, see Ref. [30]. Again, we obtain satisfactory agreement between our Pythia 8 simulation and the experimental data.  Kinematic cuts are implemented in our simulations as in CMS measurements, see Ref. [22]. The roughness of the b-tagged dijet cross section relative to that for inclusive dijets is due to the inherently lower statistics.
With the confidence that our comparisons to experimental measurements afford us, we now present the detailed baseline information for b-tagged and inclusive dijet production in p+p collisions, at √ s = 5.02 TeV for the LHC and √ s = 200 GeV for RHIC. These are the same center-of-mass energies (per nucleon pair) for the current heavy ion collisions at the LHC and for the planned sPHENIX experiment, respectively.
In Fig. 2, we show the three-dimensional (3D) plots of the cross section (weighted by the transverse momenta p 1T p 2T ) at the LHC energy √ s = 5.02 TeV as a function of the transverse momenta of the two jets (p 1T and p 2T ) in the mid-rapidity region |y| < 2. The jets are reconstructed with a jet radius R = 0.4. Here we label the dijet transverse momenta as p 1T and p 2T (instead of p L T and p S T ), because we do not distinguish which jet is leading or Kinematic cuts implemented in our simulations are the same as those from the sPHENIX collaboration [31]. Here again, the slight bumpiness of the b-tagged dijet cross section is due to its lower statistics relative to the inclusive case.
subleading in making the 3D plots. We will follow such a convention in the rest of the paper: when we need to specify leading and subleading jets, we label them as p L T and p S T . Otherwise, we simply label them as p 1T and p 2T . The left plot is for b-tagged dijet production, while the right is for inclusive dijets. Fig. 3 is the same as Fig. 2, but for sPHENIX energy √ s = 200 GeV. The roughness of the b-tagged dijet cross section relative to that for inclusive dijets is due to the inherently lower statistics. As usual [20], the cross section reaches its maximum for p 1T ≈ p 2T , and is broad and slowly varying as one goes away from this main diagonal. Such features will help us understand the behavior of nuclear modification in heavy ion collisions as we will see below. Let us now turn to the flavor origin of the dijets, which will be of central importance for our simulations of medium effects in heavy ion collisions, presented in the next section. The detailed kinematic constraints are shown in each plot. Pythia 8 utilizes leading order (LO) perturbative QCD matrix elements combined with parton showers. For b-tagged dijet production, there are 7 channels in our simulations: g + g → b +b, q +q → b +b, g + g → g + g, q +q → g + g, q + g → q + g, g + g → q +q, q + q → q + q. We classify these 7 channels to 4 subprocesses according to the flavor information of the final state partons in LO matrix elements: g + g includes the contributions from category (2), with g + g in the final state. In this case, both b-tagged jets are initiated by prompt gluons through g → b +b splitting in the showering process. Thus, the medium modification of these b jets would resemble that of a massive gluon of effective mass 2m b . Similarly, the red curve denotes the process from category (3). Thus, one b-jet is initiated by a gluon g like above. On the other hand, the other b-jet is initiated by a light quark q, for which the medium modification would resemble that of a massive quark. Finally, the black curve denotes the processes in category (4). In this case, both of the b-tagged jets are initiated by light quarks q. As we can see, for a wide kinematic coverage, the subprocesses with b +b in the final state provide the dominant contributions ( > ∼ 50%) to b-tagged dijet production in p+p collisions at the LHC at √ s = 5.02 TeV. On the other hand, the b +b channel dominates b-tagged dijet production across the p T range above 10 GeV, which is the relevant range for the sPHENIX experiment. This indicates that b-tagged dijet production provides an excellent opportunity to study the effects of heavy quark energy loss in heavy ion collisions. On the other hand, for inclusive dijet production, the usual 5 partonic processes will be reclassified into three subprocesses through their final state parton contents: (1) g + g → g + g, q +q → g + g; (2) q + g → q + g; (3) g + g → q +q, q + q → q + q. In category (1), both jets are initiated by gluons, while for category (3), both jets are initiated by quarks. For category (2), the dijets are initiated by a light quark q and a gluon g, respectively. One can clearly see in Fig. 6 that at LHC energy √ s = 5.02 TeV, for a large kinematic region, the process from category (2) is the dominant channel for inclusive dijet production. In other words, inclusive dijets at LHC kinematics are mostly initiated by a quark q on one side and a gluon g on the other end of the azimuthal plane. In addition, we plot such , respectively. We find that at relatively lower jet transverse momenta ( < ∼ 20 GeV), the inclusive dijet cross section is dominated by category (2) with q + g in the final state. At the higher jet transverse momenta, the cross section is dominated by category (3) with q + q in the final state. This is expected since as the jet transverse momenta increase, the parton momentum fractions x in the protons reach the region x ∼ 1, where valence quarks dominate.

III. LIGHT AND HEAVY FLAVOR DIJET PRODUCTION IN HOT QCD MATTER
In this section, we provide the main formula and basic information on how we implement parton energy loss for both inclusive and b-tagged dijet production in heavy ion collisions.

A. Dijet production: main formula
Our starting point for both p+p and A+A collisions is the double differential cross section, dσ/dp 1T dp 2T , in twodimensional transverse momentum bins (p 1T , p 2T ) of the leading and subleading jets. With such a double differential cross section at hand, one can compute the dijet invariant mass distribution, as well as the so-called imbalance distribution as follows.
The dijet invariant mass m 2 12 = (p 1 + p 2 ) 2 can be written in terms of the jet transverse momentum and rapidity as follows where m 2 1 = p 2 1 and m 1T = m 2 1 + p 2 1T are the invariant mass squared and the transverse mass for one of the jets, likewise we have m 2 and m 2T for the other jet. At the same time, we have the difference in the rapidities and the azimuthal angles as where η 1,2 and φ 1,2 are the rapidities and azimuthal angles for the jets. In the relevant kinematic regimes where the transverse momentum is much larger than the jet mass, p T m, we approximate m T ≈ p T and obtain In the actual Pythia 8 simulations for dijet production, we generate the averaged m 2 1 , m 2 2 , and cosh(∆η)−cos(∆φ) for each (p 1T , p 2T ) bin. With this information, we compute the dijet invariant mass distribution through the double differential dijet momentum distribution via the following formula dσ dm 12 = dp 1T dp 2T dσ dp 1T dp 2T δ m 12 where the transverse momenta p 1T and p 2T are integrated over the desired experimental cuts. Kinematic cuts are implemented in our simulations as in CMS measurements [22]. The black histograms represent dσ/dm12 simulated directly from Pythia 8, while the red histograms are dσ/dm12 computed using Eq. (6) and Pythia 8-simulated dσ/dp1T dp2T . Let us now confirm that such a procedure yields correct dijet mass distributions. To do this, we compare the dijet invariant mass distribution indirectly computed using Eq. (6) and Pythia 8-simulated dσ/dp 1T dp 2T , with dσ/dm 12 simulated directly from Pythia 8. We perform such a comparison in Figs. 8 and 9 for b-tagged (left panel), as well as for inclusive (right panel), dijet production at LHC energy √ s = 5.02 TeV and sPHENIX energy √ s = 200 GeV, respectively. The black histograms represent dσ/dm 12 simulated directly from Pythia 8, while the red histograms are dσ/dm 12 computed using Eq. (6) and Pythia 8-simulated dσ/dp 1T dp 2T . We observe a near-perfect matching mass spectra obtained via direct implementation of dijet mass in Pythia 8 and our approximate formula in Eq. (6), as indicated by the fact that the black and red histograms overlay each other. This validates the use of our formula in applications to heavy-ion collisions and subsequent dijet mass modifications.
On the other hand, one of the more conventional observables, the dijet momentum imbalance shift, is based on the cross section as a function of the imbalance variable which can be derived from the double differential cross section dσ/dp 1T dp 2T . The formula is given as follows dσ dz J = dp 1T dp 2T dσ dp 1T dp 2T where again, the limits of integration for p 1T and p 2T are matched to the desired experimental cuts. Comparing Eqs. (6) with (8), one can gain some insights why medium modification of dijet invariant mass distribution leads to enhanced medium effects than that of the dijet momentum imbalance. This is because dijet invariant mass m 12 ∝ p 1T p 2T , i.e. product of two jet momenta, and thus leads to an addition of the jet quenching effects on the individual jets. On the other hand, the momentum imbalance z J = p 2T /p 1T , i.e. quotient of two jet momenta, and thus leads to an subtraction of the jet quenching effects on the individual jets. We elaborate more on this point in the presentation of our numerical results below.

B. Modification of dijet production
In the presence of the hot and dense QCD medium, the vacuum parton shower gets modified due to the radiative [32][33][34][35][36][37][38][39] and collisional [40][41][42][43][44][45] energy losses of the propagating partons that initiate and form the jets. The implementation of energy loss effects in heavy ion collisions is explained in detail in, e.g., Refs. [20,46]. For a given impact parameter |b ⊥ | in the transverse plane of the nucleus-nucleus collisions, we evaluate the inclusive dijet double differential cross sections in (p 1T , p 2T ) as follows where |b ⊥ | is the mean impact parameter for a given collision centrality. In the b-tagged dijet case, we further include the b-quark contributions that initiates prompt b-jets. In Eq. P q,g ( ) is the probability density for the parent parton to redistribute a fraction of its energy through medium-induced soft gluon bremsstrahlung. For reconstructed jets, what matters is the out-of-cone energy loss fraction f loss q,g [46] f loss q,g (R; rad + coll) which includes both radiative and collisional energy loss effects, with ω min being a parameter that controls the energy dissipated by the medium-induced parton shower into the QGP due to collisional processes [45]. On the other hand, dN g q,g (ω,r) dωdr is the medium-induced gluon distribution [48], which is the soft emission limit of the complete in-medium splitting functions [49].
Splitting functions themselves are calculated using the formula derived in the SCET M,G framework [50]. They have been independently obtained in the lightcone wavefunction approach [51] for both massless and massive partons, and are evaluated in the QGP medium simulated by the iEBE-VISHNU code package [52]. The same model of the medium has been recently used to calculate quarkonium suppression at the LHC [53] and soft-drop groomed momentum sharing distributions [12]. Numerical evaluation of the splitting functions requires multi-dimensional integration over the jet production point, the propagation of the jet in matter, and the transverse momentum dependence of the jet-medium cross section. Since the integral dimension is larger than 4, we use a numerical integration based on the Monte Carlo method. In particular, the VEGAS algorithm [54] implemented in the CUBA multidimensional numerical integration library [55] is used because the adaptive importance sampling algorithm is efficient for integrands with localized peaks. The splitting function calculation code is written in C++. The integrals are evaluated on a Xeon cluster with task parallelization for different kinematic variables such as energy, momentum, quark mass, or the splitting channel, utilizing multiple CPU cores. Integration ranges are determined following the study presented in Ref. [49].
Once we obtain the medium-modified differential cross section dσ AA /dp 1T dp 2T , we then use Eqs. (6) and (8) to compute the dijet invariant mass distribution dσ AA /dm 12 and imbalance distribution dσ AA /dz J in heavy ion collisions. Such a procedure is perfectly fine for dσ AA /dz J , but is an approximation for dσ AA /dm 12 , where we assume that the medium modification for the single jet mass distributions m 2 1 and m 2 2 are much smaller than those for the transverse momenta p 1T and p 2T . Thus, starting from Eq. (6), we obtain dσ AA dm 12 = dp 1T dp 2T dσ AA dp 1T dp 2T δ m 12 − m 2 1 pp + m 2 2 pp + 2p 1T p 2T cosh(∆η) − cos(∆φ) pp , where we have used the same values for m 2 1 , m 2 2 , and cosh(∆η) − cos(∆φ) as those in p+p collisions, as denoted by the subscript pp. Such an approximation is well-justified. For example, mass distributions for single inclusive jets are indeed not significantly modified, as observed by ALICE collaboration at the LHC [56].

IV. PHENOMENOLOGICAL RESULTS AT RHIC AND THE LHC
In this section we first present our phenomenological results for both inclusive and b-tagged dijet production in A+A collisions at the LHC, as well as the future sPHENIX experiment at RHIC. To investigate dijet production in heavy ion collisions and quantify its deviation from the baseline results in elementary p+p reactions, we start with the two-dimensional nuclear modification factor dσ AA (|b ⊥ |)/dp 1T dp 2T dσ pp /dp 1T dp 2T , where |b ⊥ | is the corresponding impact parameter and N bin is the average number of nucleon-nucleon scatterings for a given centrality class. In this paper, we focus on the most central collisions. In Fig. 10, we make 3D plots for nuclear modification factor R AA as a function of the jet transverse momenta p 1T and p 2T simultaneously. The calculations are done for the production of dijets with radii R = 0.4 in central (0 − 10%) Pb+Pb collisions at the LHC energy √ s N N = 5.02 TeV. We integrate the rapidities of both jets over the interval |y| < 2. For the medium effects, we choose the coupling between the jet and the medium to be g = 1.8. This is consistent with the value used in our previous studies for single inclusive jets [57], vector-boson-tagged jets [46], jet substructure [58,59], and single inclusive hadrons [50,60,61] in A+A collisions. The left figure is for b-tagged dijet production, while the right is for inclusive dijets. We note that while we plot the full symmetric range in p 1T and p 2T , we do have in mind that the first jet (1) will be the trigger or leading jet and the second jet (2) will be the recoil or subleading jet. Thus, we incorporate on average path length and color charge bias effects in our calculation. As one can clearly see, the largest suppression occurs along the diagonal p 1T = p 2T , consistent with our expectation. In the region away from the diagonal, there is a striking enhancement. As the future sPHENIX [13] experiment will have good sensitivity in measuring both inclusive and b-tagged dijet production, it is an opportune time to make predictions for sPHENIX kinematics. In Fig. 11 we make similar 3D plots of R AA for b-tagged (left) and inclusive (right) dijet production at sPHENIX energy √ s N N = 200 GeV. Kinematic cuts implemented in our simulations are the same as those from the sPHENIX collaboration [31]. Obviously the kinematic coverage for the jet transverse momenta is much smaller than that of the jets at the LHC, due to a much smaller center-of-mass energy. However, the suppression is even stronger along the diagonal p 1T = p 2T . This is simply because the cross sections at RHIC energies fall much faster as functions of jet transverse momenta due to limited phase space, and thus jet quenching effects get amplified [62][63][64][65][66].
If such two-dimensional nuclear modification ratios could be measured in detail, they would provide the most information and insight into jet quenching and heavy flavor dynamics in the medium. However, the statistics necessary to perform such measurements make this, at present, quite difficult. In practice, one usually integrates out one of the differential variables and, thus, achieves a one-dimensional nuclear modification ratio. In this respect, the conventional dijet momentum imbalance z J and asymmetry A J distributions have been extensively studied in the literature. The medium modification on these traditional distributions emphasize the difference in the quenching of the dijet production, which has been observed to be relatively small. We will present such studies toward the end of this section. Here instead, we present the nuclear modification for another observable, the dijet invariant mass distribution, defined as follows Again, the impact parameter |b ⊥ | indicates the centrality class for the A+A collisions. The numerator and denominator are the dijet mass distribution in A+A and p+p collisions, respectively. They are computed through the double differential cross sections dσ/d 1T dp 2T as in Eqs. (11) and (6), respectively. In Eqs. (6) and (11), one can immediately see the advantage of such an observable. First, being only differential in the dijet invariant mass m 12 , it is a one-dimensional observable, hence one should have enough statistics to perform these measurements experimentally. Second, since the dijet invariant mass is proportional to the product of the dijet transverse momenta, as can be clearly seen in Eq. (5), the dijet mass distribution incorporates the medium modification of the dσ/d 1T dp 2T in an additive way, as we have already emphasized in Sec. III A. In other words, compared to the traditional momentum asymmetry observables, the dijet mass distribution adds rather than subtracts the medium modifications of the two jets. Naturally, one would expect the medium modification of dijet mass distributions to be greatly enhanced and thus to be more sensitive to the properties of the medium. In Fig. 12, we plot the nuclear modification factor R AA as a function of dijet invariant mass m 12 for inclusive (left) and b-tagged (right) dijet production in Pb+Pb collisions at √ s N N = 5.02 TeV at the LHC. For inclusive dijet production, the band corresponds to a range of coupling strengths between the jet and the medium: g med = 1.8 − 2.0. On the other hand, for b-tagged dijet production, we fix g med = 1.8, and the band corresponds to a range of masses of the propagating system between m b and 2m b , implemented as detailed in [9]. We make transverse momentum cuts requiring both leading and subleading jets to have p L,S T > 30 GeV. This is why we have a lower limit on the dijet invariant mass m 12 > ∼ 100 GeV in these plots. As one can clearly see from the figures, being an additive effect, R AA can be as small as 0.1, i.e., suppressed by a factor of 10 in the lower end of the invariant mass m 12 ∼ 100 GeV. This is a dramatic suppression, much stronger than the suppression for single inclusive jet production, around a factor of 2 [57]. As one increases the invariant mass m 12 , the suppression gets smaller, but it is still around a factor of 2 or more even at m 12 ∼ 500 GeV. The suppression for b-tagged dijet production is smaller than that of inclusive dijets at smaller dijet mass m 12 ∼ 100 GeV, and becomes similar to inclusive dijet production as m 12 increases. This is to be expected, as heavy quark mass effects on jet quenching are more important at lower transverse momenta, or naturally smaller dijet invariant mass.
In Fig. 13, we present the same plots but for Au+Au collisions at √ s N N = 200 GeV, relevant to the sPHENIX experiment at RHIC. For inclusive dijet production, the band corresponds to a range of coupling strengths between the jet and the medium: g med = 2.0 − 2.2. On the other hand, for b-tagged dijet production, we fix g med = 2.0, and the band again corresponds to a range of masses of the propagating system between m b and 2m b . We choose a slightly larger coupling strength at RHIC compared to that for the above LHC kinematics, which is also consistent with our previous studies and that of the JET collaboration [67]. Since the center-of-mass energy is much lower, we select jets with much lower p T > ∼ 8 GeV, and correspondingly lower dijet invariant mass m 12 > ∼ 20 GeV for RHIC kinematics. Having smaller jet transverse momenta and cross sections that fall off strongly as functions of jet transverse momenta, the suppression for inclusive dijet cross sections is even larger compared with those of LHC energies. We observe a factor of ∼ 10 or more suppression even up to a relatively high invariant mass m 12 ∼ 100 GeV.
On the other hand, the suppression pattern for b-tagged dijet production as a function of m 12 at sPHENIX energy √ s N N = 200 GeV, as shown in right panel of Fig. 13, appears quite different from inclusive dijet production in left panel, and looks nothing like the b-tagged dijet production at the LHC energy in Fig. 12. It is, thus, important to understand why we observe such a behavior. If one recalls the behavior of the suppression pattern for single inclusive heavy meson/heavy quark production as a function of its transverse momentum, see, e.g. Ref. [8,50], one can understand the above behavior of R AA as a function of m 12 . Due to the heavy quark mass effect in the jet quenching formalism, R AA for heavy quark mesons first decreases and then increases when plotted as a function of p T . In other words, there is a dip in R AA as a function of p T . Now one can translate such a behavior into the behavior of R AA as a function of m 12 . For the mass region in Fig. 13, b-tagged dijets mostly fall into the relatively low values of jet transverse momenta, i.e., before the dip of R AA (as a function of p T ). This explains why R AA decreases as a function of m 12 . If one has a larger phase space to explore much higher values of transverse momenta, as is the case at the LHC energy in Fig. 12, once passing the dip of R AA , one should naturally expect R AA to increase as a function of m 12 . This is precisely what is observed in our calculations, see Fig. 12 (right). This comparison informs us that sPHENIX is sitting in a very interesting kinematic regime for testing heavy quark mass effects within the jet quenching formalism. To quantitatively compare the medium modification of b-tagged and inclusive dijet production, we further plot the ratio of nuclear modification factors for b-tagged (R bb AA ) and inclusive dijet (R jj AA ) production, R bb AA /R jj AA , as a function of dijet invariant mass m 12 in Fig. 14 In both kinematic regimes, we see a smaller suppression (thus larger R AA ) for b-tagged dijets compared to inclusive dijets, though the figure also indicates a markedly different effect at low energies than at higher ones. The most pronounced differences occur in the low mass range m 12 ∼ 20 GeV accessible by sPHENIX, where such a ratio reaches up to almost a factor of 10, R bb AA /R jj AA ∼ 10. On the other hand, at LHC energy, one should observe roughly a factor of 2 less suppression for b-tagged dijet at relatively low dijet invariant mass m 12 . For large m 12 ∼ 500 GeV, the difference diminishes and one should expect to see similar suppressions, R bb AA /R jj AA ∼ 1. Let us now turn to the conventional observable, the momentum imbalance distributions, dσ/dz J . In the absence of in-medium interactions, one expects from perturbative QCD that the transverse momenta of the two jets are balanced, p 1T ≈ p 2T . Consequently, dσ/dz J in elementary p+p collisions will be peaked around z J ≈ 1. On the other hand, in heavy ion collisions, jet quenching plays an important role and one jet will lose more energy than the other one. As a result, one expects to see a downshift of the peak in z J distribution because of strong in-medium interactions.
In Fig. 15, we display the normalized dijet imbalance z J distributions for inclusive (left) and b-tagged (right) dijet production at the LHC energy √ s N N = 5.02 TeV. The black histogram is the result for p+p collisions, while the colored curves are the results for central (0 − 10%) Pb+Pb collisions. In the left panel, the band corresponds to a range of coupling strengths between the jet and the medium: g med = 1.8 − 2.0, respectively. In the right panel, we fix g med = 1.8, and the band corresponds to a range of masses of the propagating system between m b and 2m b . The experimental data points are from CMS collaboration [22]. We clearly see a downshift in the peak of z J distribution for both inclusive and b-tagged dijet production. There is an excellent agreement between our calculations for inclusive dijets and the CMS data. On the other hand, our calculations do not describe very well the CMS data for b-tagged dijets. We attribute this to the use of purely LO matrix elements via Pythia 8 and the specific nature of the reweighting procedure carried out by CMS [22]. We do not carry out such a re-weighting procedure in order to maintain consistency with the rest of our simulations. Note that the visual difference between our results in A+A and the experimental data is also largely driven by the p+p baseline. Our calculation with g med = 2.0 appears closer to the Pb+Pb data. However, as shown below, the results with g med = 1.8 already quantitatively capture the downshift of the z J distribution in heavy ion collisions. This again emphasizes the fact that from the momentum imbalance distributions alone it might be difficult to assess whether a theoretical model correctly represents the physics of jet quenching. Fig. 16 contains the dijet imbalance z J distributions for inclusive (left) and b-tagged (right) dijet production with sPHENIX kinematics at √ s N N = 200 GeV. Our results for b-tagged dijets in p+p collisions are consistent with the preliminary simulations carried out by the sPHENIX collaboration [31] (denoted as the blue "data" points). Our calculations show that a larger shift in z J should be observed for inclusive dijets compared with b-tagged dijets.
To further quantify the downshift of the z J distribution, we define the mean value of z J , We further define the difference for z J in p+p and A+A collisions as and the positive values of ∆ z J represents downshifts of the z J distribution in A+A collisions in comparison with that of the p+p collisions. In Table I, we list our theoretical calculations for z J pp , z J AA , and ∆ z J for both inclusive and b-tagged dijet production. The values labelled as "LHC theory" are our theoretical calculations for Pb+Pb collisions at 0 − 10% centrality at √ s N N = 5.02 TeV, and can be compared with the CMS experimental data. For inclusive dijets, we perform the calculations for the coupling between the jet and the medium g med = 1.8 − 2.0, which explains the uncertainties in our theoretical values. For b-tagged dijets, we vary such a coupling in the same range while the mass of the propagating system is held fixed at m b . We find that in general the downshift ∆ z J is slightly larger for inclusive dijet production than that for b-tagged dijets, though the uncertainties are still large. Nevertheless, within the theoretical and experimental uncertainties, our theoretical calculations for all these observables z J pp , z J AA , and ∆ z J , agree well with the CMS experimental data. Finally, we also perform calculations for central Au+Au collisions for sPHENIX kinematics at √ s N N = 200 GeV in Table I, which are labelled as "sPHENIX theory." We expect such measurements will become available once the sPHENIX experiment starts running in the future.

V. CONCLUSIONS
In this paper we present detailed theoretical predictions for inclusive and b-tagged dijet production and modification in heavy ion collisions at RHIC and the LHC. We propose a new observable, the modification of dijet invariant mass, as a novel diagnostic of the quark-gluon plasma (QGP) created in ultra-relativistic heavy ion collisions. Our comprehensive studies conclusively demonstrate that this observable exhibits enhanced sensitivity to the strength of jet-medium interactions, the transport properties of nuclear matter, and to the mass effects on in-medium parton showers. Complete characterization of the quenching of multi-jet events is given by a multi-dimensional nuclear modification ratio, as we also show in this work. The statistics necessary to perform such measurements at present, however, makes them impractical even for the case of two jets. By integrating out one of those dimensions, which is usually accomplished through an auxiliary variable such as the dijet momentum asymmetry or dijet momentum imbalance, the statistics can be greatly improved and experimental results may be obtained. Such traditional momentum imbalance measurements emphasize potentially small differences in the quenching of leading and subleading jets. Hence, the shift in the mean value of the momentum imbalance variable z J is only on the order of 7 − 15%. Differences between the dijet asymmetries of inclusive and b-tagged dijets have been difficult to identify. In contrast, the dijet mass modification combines the suppression of the individual jets and enhances the observable jet quenching effect by up to an order of magnitude.
To obtain reliable predictions, we use radiative and collisional parton energy losses evaluated in a realistic hydrodynamic background. Furthermore, we utilize couplings between the jet and the medium that have successfully predicted or described the nuclear modification of a number of observables related to the reconstructed jets in heavy ion collisions -inclusive light and heavy flavor jet suppression, jet substructure, and vector-boson-tagged jets. We find that the theoretical model gives an excellent description of the most recent measurements of the inclusive dijet momentum imbalance distributions and their modification in heavy ion collisions at the LHC. For the b-tagged dijets the description is not nearly as good, though the theoretical model qualitatively and even quantitatively captures the shift of the momentum imbalance in Pb+Pb relative to p+p collisions. The deviation in shape can likely be attributed to the lowest order Pythia 8 baseline simulation.
For the main result of this paper, the dijet mass distribution modification, we find that the suppression at m 12 ∼ 100 GeV is around a factor of 10. In contrast, the suppression of single inclusive jets is only around a factor of 2. We note that as the dijet mass grows, the suppression is reduced and at m 12 ∼ 500 GeV, it is about a factor of 2. When the nuclear modification due to final-state interactions diminishes, initial-state cold nuclear matter effects may not be negligible [68]. We defer this study until after the first experimental measurements of dijet mass modification appear. In the mass region studied for the LHC, the modification of the b-tagged dijet mass distribution can be twice as small as that of inclusive dijets, though this difference disappears at larger masses.
Since sPHENIX at RHIC is expected to become available in the near future, we also perform calculations of dijet mass distributions and momentum imbalance distribution. We find that jet quenching effects on the dijet mass distribution can be significantly amplified in the kinematic range accessible by the future sPHENIX experiment, because of steeply falling spectra. In the mass region m 12 = 20 − 100 GeV, the QGP-induced suppression is a factor of 10 or larger for inclusive dijet production. A 10% change in the strong coupling constant g that describes the jet-medium interactions can lead to a factor of 2 larger suppression. On the other hand, the suppression for b-tagged dijets shows a different behavior, which can be traced back to the heavy quark mass effects. In other words, at sPHENIX kinematics, there is an enhanced sensitivity to heavy quark mass effects, and we find that in the smaller dijet mass range the suppression for b-tagged dijets can be an order of magnitude smaller.
To conclude, upcoming runs at RHIC and the LHC present compelling opportunities for experiments to explore novel jet quenching observables. The modification of light and heavy flavor dijet mass distributions will be a promising avenue of exploration in this direction.