# Iterative Prompt Refinement for Radiation Oncology Symptom Extraction
Using Teacher-Student Large Language Models

Reza Khanmohammadi  Department of Computer Science and Engineering, Michigan
State University, East Lansing, USA  Ahmed I Ghanem  Department of Radiation
Oncology, Henry Ford Cancer Institute, Detroit, USA  Alexandria Department of
Clinical Oncology, Alexandria University, Egypt  Kyle Verdecchia  Department
of Radiation Oncology, Henry Ford Cancer Institute, Detroit, USA  Ryan Hall
Department of Radiation Oncology, Henry Ford Cancer Institute, Detroit, USA
Mohamed Elshaikh  Department of Radiation Oncology, Henry Ford Cancer
Institute, Detroit, USA
Benjamin Movsas  Department of Radiation Oncology, Henry Ford Cancer
Institute, Detroit, USA  Hassan Bagher-Ebadian  Department of Radiation
Oncology, Henry Ford Cancer Institute, Detroit, USA  Departments of Radiology
and Osteopathic Medicine, Michigan State University, East Lansing, USA
Department of Physics, Oakland University, Rochester, USA  Indrin Chetty
Department of Radiation Oncology, Cedars Sinai Medical Center, Los Angeles,
USA  Department of Physics, Oakland University, Rochester, USA  Mohammad M.
Ghassemi Shared senior author. Department of Computer Science and Engineering,
Michigan State University, East Lansing, USA  Kundan Thind  Department of
Radiation Oncology, Henry Ford Cancer Institute, Detroit, USA  Department of
Medicine, Michigan State University, East Lansing, USA

Abstract This study introduces a novel teacher-student architecture utilizing
Large Language Models (LLMs) to improve prostate cancer radiotherapy symptom
extraction from clinical notes. Mixtral, the student model, initially extracts
symptoms, followed by GPT-4, the teacher model, which refines prompts based on
Mixtral’s performance. This iterative process involved 294 single symptom
clinical notes across 12 symptoms, with up to 16 rounds of refinement per
epoch. Results showed significant improvements in extracting symptoms from
both single and multi-symptom notes. For 59 single symptom notes, accuracy
increased from 0.51 to 0.71, precision from 0.52 to 0.82, recall from 0.52 to
0.72, and F1 score from 0.49 to 0.73. In 375 multi-symptom notes, accuracy
rose from 0.24 to 0.43, precision from 0.6 to 0.76, recall from 0.24 to 0.43,
and F1 score from 0.20 to 0.44. These results demonstrate the effectiveness of
advanced prompt engineering in LLMs for radiation oncology use.
---

##  1 Introduction

The development of Large Language Models (LLMs) has marked a transformative
era in the field of Artificial Intelligence (AI). These advanced models, built
on the foundation of extensive data and sophisticated algorithms, have
demonstrated remarkable capabilities in understanding and generating human-
like text. In the medical domain, specialized versions of these models such as
Med-PaLM 2 [1] and HuatuoGPT [2] have been tailored to understand and process
medical terminology, patient data, and clinical notes. This focus on
healthcare has enabled LLMs to assist in diagnostic processes, treatment
recommendation, and patient care, thereby revolutionizing the intersection of
AI and medicine [3]. The ability of these models to handle vast and complex
datasets has been pivotal in uncovering insights and augmenting the decision-
making process in healthcare [4]. However, the effectiveness of LLMs in
specialized areas like cancer text analysis is often limited by the
availability of relevant training data. Since these models learn and evolve
based on the data they are exposed to [5], their performance in domains with
scarce publicly available data, such as specific cancer-related texts, can be
suboptimal [6, 7]. Further, a major challenge for LLMs in specialized domains
is the lack of effective prompt engineering, which is essential for enhancing
their performance. Prompt engineering is a process that involves iterative
refinement of the input and output formats, as well as the incorporation of
domain knowledge and terminology. This process can facilitate LLMs to
supplement embedded information in specific domains, and thus increase their
applicability and effectiveness in niche fields such as oncology [8]. Recent
studies have demonstrated the efficacy of prompt engineering in improving
domain-specific inference. Wang et al. [9] showed that combining prompt
engineering with paraphrases of rare biomedical words can significantly
enhance the performance of pre-trained models in the biomedical domain.
Similarly, Alsentzer et al. [10] utilized a large language model for zero-shot
interpretable phenotyping of postpartum hemorrhage from clinical notes,
highlighting the model’s strong performance in extracting granular concepts
associated with the condition. These examples underscore the potential of
prompt engineering to refine the capabilities of LLMs in specialized medical
tasks, offering a promising avenue for the systematic abstraction of critical
concepts in fields such as radiation oncology.

Figure 1: Schematic representation of our proposed method.

Traditional methodologies in prompt refinement have varied, with some relying
on manual tuning by experts, while others explore automated techniques.
However, these methods often require a balance between human input,
computational resources, and specific task requirements [11]. This balance is
especially crucial in medical text analysis, where precision and context-
specific accuracy are paramount. Recent advancements have been observed in the
domain of radiation oncology, particularly in the extraction of key concepts
from clinical notes, including the identification of toxicity [12, 13].
Consequently, the development of prompt engineering methodologies for the
extraction of key concepts in radiation oncology offers the potential to
facilitate large-scale, systematic abstraction of critical concepts in this
field.

We introduce a novel, automatic framework for prompt refinement, using a
teacher-student dynamic between two LLMs: GPT-4 as the teacher and Mixtral as
the student. This novel architecture is shown in Figure 1. This approach is
especially pertinent in the clinical domain for key concept extraction, which
is reliant on domain specific knowledge. We demonstrate the application of
this architecture with prostate cancer radiotherapy toxicity symptom
extraction from clinical notes.

##  2 Materials and Methods

The architecture employs a novel teacher-student framework designed to
iteratively refine the prompts used by LLMs for prostate cancer radiotherapy
toxicity symptom extraction from clinical notes. The student model, Mixtral,
is tasked with the initial extraction of symptoms based on provided prompts.
The teacher model, GPT-4, subsequently refines these prompts in a series of
rounds and epochs based on the student’s performance and reasoning.

The student model, Mixtral-8x7B, operates as the primary extractor of symptoms
from the clinical notes. Given a prompt and a target symptom, it classifies
the symptom as present, negated, or unknown within the notes. The model’s
output for each note includes the symptom classification and a reasoning for
its decision.

The teacher model GPT-4, receives inputs from the student model’s performance
across all training notes for a given symptom. These inputs include:

  * •

Optimal Prompt: The prompt that has led to the highest level of accuracy in
the student’s symptom classification.

  * •

Peak Performance: The student model’s best classification accuracy acquired on
the optimal prompt.

  * •

Decision-making Reasoning: Insights into the reasoning behind the student
model’s classification tied to the optimal prompt.

  * •

Failed Prompt History: A list of prompts previously tested within the current
epoch, which have yielded inferior results compared to the optimal prompt.

Using these inputs, GPT-4 generates refined prompts intended to improve the
student’s performance in symptom extraction. This refinement process occurs
through a sequence of rounds and epochs.

Full optimization encompasses a series of epochs, each consisting of 16
rounds. In every epoch, the teacher model iteratively refines prompts to
improve the student model’s symptom extraction performance. At each epoch’s
start, the teacher is updated with the optimal prompt and peak performance
metrics from the previous epoch where it has 16 opportunities (one per round)
to refine the prompt and exceed the student’s current best performance,
measured by accuracy. If, in any of the 16 rounds, the teacher model enhances
the student’s performance beyond its prior peak, the epoch concludes, and the
process advances to the next epoch where the teacher is granted a reset of 16
refinement opportunities. The refinement process terminates if the teacher
model fails to improve the student’s performance in all 16 rounds of an epoch,
signaling the achievement of the student’s maximum capability within the
framework’s parameters.

Symptom | Single-symptom Notes | Multi-symptom Notes |  |  |  |
---|---|---|---|---|---|---
Accuracy | Precision | Recall | F1 | Accuracy | Precision | Recall | F1
| Init. | Post. | Init. | Post. | Init. | Post. | Init. | Post. | Init. | Post. | Init. | Post. | Init. | Post. | Init. | Post.
Urgency | $0.80$ | 1.00 | $0.87$ | 1.00 | $0.80$ | 1.00 | $0.80$ | 1.00 | $0.37$ | 0.52 | $0.72$ | 0.69 | $0.37$ | 0.52 | $0.27$ | 0.49
Cystitis | $0.40$ | 0.60 | $0.80$ | 0.90 | $0.40$ | 0.60 | $0.53$ | 0.67 | $0.04$ | 0.38 | $0.27$ | 0.78 | $0.04$ | 0.38 | $0.06$ | 0.49
Urinary Obs. | $0.40$ | 0.80 | $0.16$ | 0.67 | $0.40$ | 0.80 | $0.23$ | 0.72 | $0.01$ | 0.12 | 0.98 | $0.97$ | $0.01$ | 0.12 | $0.01$ | 0.19
Dysuria | $0.80$ | 0.80 | $0.65$ | 0.70 | $0.80$ | 0.80 | $0.71$ | 0.73 | $0.56$ | 0.62 | $0.43$ | 0.67 | $0.56$ | 0.62 | $0.44$ | 0.61
Erectile Dys. | 1.00 | $0.80$ | $1.00$ | 1.00 | 1.00 | $0.80$ | 1.00 | $0.89$ | 0.26 | $0.21$ | $0.38$ | 0.62 | 0.26 | $0.21$ | 0.19 | $0.12$
Rectal Bleed. | $0.00$ | 0.20 | $0.00$ | 0.20 | $0.00$ | 0.20 | $0.00$ | 0.20 | $0.16$ | 0.57 | $0.16$ | 0.76 | $0.17$ | 0.57 | $0.10$ | 0.58
Stricture | $0.40$ | 1.00 | $0.16$ | 1.00 | $0.40$ | 1.00 | $0.23$ | 1.00 | $0.04$ | 0.15 | 0.96 | $0.83$ | $0.04$ | 0.15 | $0.02$ | 0.20
Nocturia | 1.00 | $0.80$ | 1.00 | $0.64$ | 1.00 | $0.80$ | 1.00 | $0.71$ | $0.30$ | 0.34 | $0.28$ | 0.75 | $0.30$ | 0.34 | $0.25$ | 0.30
Proctitis | $0.00$ | 0.20 | $0.00$ | 0.80 | $0.00$ | 0.20 | $0.00$ | 0.32 | $0.01$ | 0.20 | $0.74$ | 0.77 | $0.01$ | 0.20 | $0.02$ | 0.32
Hematuria | $0.60$ | 0.60 | $0.80$ | 1.00 | $0.60$ | 0.60 | $0.57$ | 0.70 | $0.72$ | 0.76 | $0.63$ | 0.75 | $0.72$ | 0.76 | $0.66$ | 0.71
Urothelial Carc. | $0.00$ | 1.00 | $0.00$ | 1.00 | $0.00$ | 1.00 | $0.00$ | 1.00 | $0.05$ | 0.70 | $0.99$ | 0.99 | $0.05$ | 0.70 | $0.08$ | 0.81
Incontinence | $0.80$ | 0.80 | $0.87$ | 0.87 | $0.80$ | 0.80 | $0.80$ | 0.80 | $0.39$ | 0.53 | $0.64$ | $0.58$ | $0.39$ | 0.53 | $0.30$ | 0.51
Average Score | $0.51$ | 0.71 | $0.52$ | 0.82 | $0.51$ | 0.71 | $0.49$ | 0.73 | $0.24$ | 0.43 | $0.60$ | 0.76 | $0.24$ | 0.43 | $0.20$ | 0.44
Table 1: Symptom extraction test results: Initial (Init.) vs Post-refinement
(Post) Algorithm 1 Iterative Prompt Refinement for Symptom Extraction on a
target symptom

0: Target symptom, Initial prompt, MaxEpochs, MaxRounds

1: Initialize Mixtral-8x7B (student model) with initial prompts

2: Initialize GPT-4 (teacher model)

3: Set best performance metric to 0

4: Set optimal prompt to initial prompt for the symptom

5: for epoch = 1 to MaxEpochs do

6: for round = 1 to MaxRounds do

7: Apply current optimal prompt to the student model across all clinical notes

8: Gather labels and reasonings

9: Calculate performance score

10: if performance improves over best performance then

11: Update best performance and optimal prompt

12: Break and proceed to the next epoch

13: end if

14: Use the teacher model to refine the optimal prompt based on:

  * •

Current optimal prompt

  * •

Best performance metrics

  * •

Reasoning of the student model

  * •

History of failed prompts

15: end for

16: if no improvement in any of the 16 rounds then

17: Terminate refinement for the symptom

18: end if

19: end for

The data consisted of clinical notes from a cohort of prostate cancer patients
treated with a radiotherapy (RT) dose of 78 Gy between 2013 and 2020. Clinical
notes 6 months post-RT were extracted to focus on long term toxicities. 294
notes with presence of single symptoms were selected. Twelve common post-RT
toxicity symptoms: Urgency, Cystitis, Urinary Obstruction, Dysuria, Erectile
Dysfunction, Rectal Bleeding, Stricture, Nocturia, Proctitis, Hematuria,
Urothelial Carcinoma, and Incontinence were selected. For the training set, we
utilized 20 notes per symptom, except for Urothelial Carcinoma, where only 15
notes were available. Correspondingly, the test set comprised 5 notes for each
symptom, with the exception of Urothelial Carcinoma, which had 4. These notes
were labelled with symptom assignment of -1, 0, or 1, signifying its absence,
unknown status, or presence, respectively. A separate cohort of 375 clinical
notes with multiple symptom within each note were also labelled and used
exclusively for validation. Multi-symptom notes ranged from 1 to 10 for the
number of distinct symptom mention in each notes, with an average of 5
symptoms mentions per mutli-symptom note. This comprehensive dataset structure
and labeling approach provided a robust foundation for evaluating the
architecture performance, particularly in demonstrating its scalability and
effectiveness in handling complex, multi-symptom clinical notes. Use of
clinical notes, and the overall study was reviewed and approved by the
institutional review board.

The evaluation methodology utilized accuracy as the primary metric for
assessing the student model’s performance during the prompt refinement
process. To measure the effectiveness of these refined prompts compared to the
initial ones, we employed accuracy along with weighted precision, recall, and
F1 scores on the test set. These metrics provided a comprehensive view of the
refined prompts’ performance enhancements over the initial versions. You can
find the algorithmic implementation of our method in Algorithm 1.

##  3 Results

Figure 2: The evolution of Mixtral’s symptom extraction accuracy per different
symptoms at different epochs.

Table 1, summarizes accuracy, precision, recall and F1 scores on both single
symptom notes, and multi-symptom notes. Initial and post-refinement refinement
results are presented for both categories, with initial performance metrics
based on the the initial prompt, and post-refinement performance based on
optimized prompt. Improvements across all evaluation metrics is observed.
Highlighted entries show performance improvement from initial to post-refined
extraction. For single symptom notes, a strong performance increase is
observed across all metrics. The average improvement in accuracy from 0.51 to
0.71, precision from 0.52 to 0.82, recall from 0.52 to 0.72, and F1 from 0.49
to 0.73 is observed in single symptom notes. Noteworthy improvements include
the symptom "Urgency," where accuracy soared from 0.80 to a perfect 1.00, and
precision mirrored this perfect score. Similarly, "Stricture" and "Urothelial
Carcinoma" witnessed accuracy enhancements from 0.40 and 0.00 to 1.00,
respectively. Performance improvement across all metrics was observed for
multi-symptom notes, with average improvement in accuracy from 0.24 to 0.43,
precision from 0.6 to 0.76, recall from 0.24 to 0.43 and F1 from 0.20 to 0.44.
The improvement was was not as strong as the single-symptom notes. Notable
example of "Rectal Bleeding," showcased an increase in accuracy from 0.16 to
0.57 and an F1 score increase from 0.10 to 0.58.

Figure 2 depicts the progression of Mixtral’s symptom extraction accuracy
across different symptoms over successive epochs, as well as the average and
standard deviation of symptom accuracy per epoch. A clear trend of improvement
is observed in nearly all symptoms as the number of epochs increases,
illustrating the efficacy of our prompt refinement process. Notably, symptoms
such as "Urothelial Carcinoma" and "Proctitis" show a remarkable upward
trajectory, indicating significant gains from the iterative approach. This
demonstrates adaptive learning capabilities of this novel architecture.

##  4 Discussion

This novel student LLM-teacher LLM architecture demonstrates performance
improvement in accuracy, precision, recall and F1 scores for prostate cancer
RT symptoms. This novel architecture makes contribution in the following key
areas:

  * •

Automatic Prompt Engineering This automated method for refining prompts using
a teacher-student combination of two LLMs, streamlines the traditionally
manual process of prompt optimization.

  * •

Zero-Shot Learning This framework enhances LLM performance without additional
model training, showcasing the efficiency of our zero-shot learning approach.

  * •

Local Model Improvement for Privacy The improvement of the student model is
conducted locally, ensuring sensitive clinical data, like Protected health
information (PHI), is processed securely, addressing major data privacy
concerns in healthcare.

The performance improvement across metrics of accuracy, precision, recall and
F1 was observed for both single and multi-symptom notes. The prompt-
optimization was performed on single symptom notes, and applied across test
set of single symptom notes, and multi-symptom notes. This could potentially
have contributed to inferior performance improvement in multi-symptom notes.
Limited sample size of single symptom notes utilized for optimization may
potentially have resulted in over fitting. Higher sample size for
optimization, both single and multi-symptom could potentially further improve
robustness, and the performance. Hyperparameter selection for both student and
teacher models can also have a significant impact on the overall optimization
process. Repetition studies were performed to tune hyperparameters of student
model to ensure absolute consistency in response over three repetitions. As
such, low temperature (0.2) and top_p (0.1) were utilized for the student
model, and high temperature (2) and top_p (0.9) were used for the teacher
model in this architecture. This may have limitations and further
experimentation may lead to better results.

##  5 Conclusion

In conclusion, this novel architecture presents a pathway for self-optimized
local LLM agents that can extract key concepts in medical notes with a zero
shot learning approach. The complete local inference by the student model is
suitable for protecting privacy of healthcare data. This system provides
proof-of-concept in the research setting, and is promising for the potential
of natural language processing methods to support clinical care.

## References

  * [1] Karan Singhal et al.  “Towards Expert-Level Medical Question Answering with Large Language Models”, 2023  arXiv:2305.09617 [cs.CL]
  * [2] Hongbo Zhang et al.  “HuatuoGPT, towards Taming Language Model to Be a Doctor”, 2023  arXiv:2305.15075 [cs.CL]
  * [3] Mert Karabacak and Konstantinos Margetis  “Embracing Large Language Models for Medical Applications: Opportunities and Challenges”  In _Cureus_ 15, 2023, pp. e39305  DOI: 10.7759/cureus.39305
  * [4] Niloufar Eghbali, Daniel Siegal, Chad Klochko and Mohammad Ghassemi  “Automation of Protocoling Advanced MSK Examinations Using Natural Language Processing Techniques”  In _AMIA Joint Summits on Translational Science proceedings. AMIA Joint Summits on Translational Science_ 2023, 2023, pp. 118–127
  * [5] Mehrsa Mardikoraem, Zirui Wang, Nathaniel Pascual and Daniel Woldring  “Generative models for protein sequence modeling: recent advances and future directions”  In _Briefings in Bioinformatics_ 24.6, 2023, pp. bbad358  DOI: 10.1093/bib/bbad358
  * [6] Reza Khanmohammadi et al.  “An Introduction to Natural Language Processing Techniques and Framework for Clinical Implementation in Radiation Oncology”, 2023  arXiv:2311.02205 [cs.CL]
  * [7] Danial Kamali et al.  “Using Persuasive Writing Strategies to Explain and Detect Health Misinformation”, 2023  arXiv:2211.05985 [cs.CL]
  * [8] Bertalan Meskó  “Prompt Engineering as an Important Emerging Skill for Medical Professionals: Tutorial”  In _Journal of Medical Internet Research_ 25, 2023, pp. e50638  DOI: 10.2196/50638
  * [9] Haochun Wang et al.  “Prompt Combines Paraphrase: Teaching Pre-trained Models to Understand Rare Biomedical Words”  In _Proceedings of the 29th International Conference on Computational Linguistics_ Gyeongju, Republic of Korea: International Committee on Computational Linguistics, 2022, pp. 1422–1431  URL: https://aclanthology.org/2022.coling-1.122
  * [10] Emily Alsentzer et al.  “Zero-shot Interpretable Phenotyping of Postpartum Hemorrhage Using Large Language Models”  In _medRxiv_ Cold Spring Harbor Laboratory Press, 2023  DOI: 10.1101/2023.05.31.23290753
  * [11] Banghao Chen, Zhaofeng Zhang, Nicolas Langrené and Shengxin Zhu  “Unleashing the potential of prompt engineering in Large Language Models: a comprehensive review”, 2023  arXiv:2310.14735 [cs.CL]
  * [12] Miller T Savova GK. Bitterman DS Goldner E Finan S Harris D  “An End-to-End Natural Language Processing System for Automatically Extracting Radiation Therapy Events From Clinical Texts”  In _Int J Radiat Oncol Biol Phys._ 117.3, 2023, pp. 262–273  DOI: 10.1016/j.ijrobp.2023.03.055
  * [13] Julian C Hong et al.  “Natural language processing for abstraction of cancer treatment toxicities: accuracy versus human experts”  In _JAMIA Open_ 3.4, 2020, pp. 513–517  DOI: 10.1093/jamiaopen/ooaa064


# Bayesian inference of the path-length dependence of jet energy loss

Jordan Wu<EMAIL_ADDRESS>Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA  Physics Department,
University of California, Berkeley, California 94720, USA  Weiyao Ke
<EMAIL_ADDRESS>Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  Physics Department, University of
California, Berkeley, California 94720, USA  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM 87545, USA  Xin-Nian Wang<EMAIL_ADDRESS>Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley,
California 94720, USA  Physics Department, University of California, Berkeley,
California 94720, USA

###### Abstract

A simple model for medium modification of the jet function can be used to
extract the jet energy loss distribution through a parameterized form. We
carry out a comprehensive Bayesian analysis of the world data on single
inclusive jet spectra in heavy-ion collisions at both RHIC and LHC energies.
We extract the average jet energy loss $\langle\Delta E\rangle$ as a function
of jet transverse momentum $p_{T}$ for each collision system and centrality
independently. Assuming jet energy loss is proportional to the initial parton
density $\rho\sim dN_{\rm ch}/d\eta/\pi R_{\rm eff}^{2}$ as estimated from the
pseudorapidity density of charged hadron multiplicity $dN_{\rm ch}/d\eta$ and
the effective system size $R_{\rm eff}\sim N_{\rm part}^{1/3}$ given by the
number of participant nucleons $N_{\rm part}$, the scaled average jet energy
loss $\langle\Delta E\rangle/\rho\sim R_{\rm eff}^{0.59}p_{T}^{0.13}\ln p_{T}$
for jet cone-size $R=0.4$ is found to have a momentum dependence that is
slightly stronger than a logarithmic form while the system size or length
dependence is slower than a linear one. The fluctuation of jet energy loss is,
however, independent of the initial parton density or the system size.

Jet energy loss, machine learning, Markov Chain Monte Carlo, Bayesian analysis

##  I Introduction

Parton energy loss in dense medium was predicted to lead to the suppression of
large transverse momentum hadrons and jets, known as jet quenching, in high-
energy heavy-ion collisions Gyulassy and Plumer (1990); Wang and Gyulassy
(1992). Jet quenching was indeed observed in experiments at both the
Relativistic Heavy-ion Collider (RHIC) Adcox  _et al._ (2002); Adler  _et al._
(2002); Wang (2005) and Large Hadron Collider (LHC) Aad  _et al._ (2010);
Aamodt  _et al._ (2011a); Chatrchyan  _et al._ (2012, 2012); Khachatryan  _et
al._ (2017); Acharya  _et al._ (2020). Phenomenological studies have extracted
the jet transport coefficient from comparisons between experimental data and
model calculations Chen  _et al._ (2010); Burke  _et al._ (2014); Cao  _et
al._ (2021); Liu  _et al._ (2022); Apolinário  _et al._ (2022); Xie  _et al._
(2023, 2022) whose values point to the formation of the quark-gluon plasma
with extremely high temperatures in high-energy heavy-ion collisions. These
data-model comparisons are all based on perturbative QCD calculations of
collisional and radiative parton energy loss inside QGP. In a static and
uniform QGP medium, the total radiative parton energy loss is predicted to
have a quadratic path-length dependence because of the non-Abelian Landau-
Pomeranchuck-Migdal interference in gluon radiation induced by multiple
scattering Baier  _et al._ (1997); Zakharov (1996); Wiedemann (2000); Gyulassy
_et al._ (2000); Guo and Wang (2000); Wang and Guo (2001); Majumder (2012). It
is also proportional to the jet transport coefficient, which is proportional
to the local blue color charge density. Taking into account the rapid
longitudinal and transverse expansion in high-energy heavy-ion collisions, the
effective total parton energy loss averaged over the azimuthal angle is
approximately linear in the system size or the average escape time Gyulassy
_et al._ (2002). Such a system size dependence is consistent with the
effective parton energy loss extracted from the measured suppression of single
inclusive hadrons Arleo and Falmagne (2022).

One can define the effective energy loss for a reconstructed jet with a given
jet cone size $R$ as the difference between the energy of a jet in
proton+proton ($p$+$p$) collisions and the energy of the final jet in heavy-
ion ($A$+$A$) collisions that originates from the same hard process in $p$+$p$
collisions. The relation between the jet energy loss and that of an individual
parton is not straightforward since some of the radiated gluons can end up
inside the jet cone as part of the final jet. In addition, some of the lost
energy carried by the recoil medium partons as part of the jet-induced medium
response can also contribute to the final jet energy with a given jet cone
size $R$. The momentum and system size dependence of the jet energy loss can
only be studied through calculations that take into account of both the above
effects He  _et al._ (2015); Wang  _et al._ (2017); He  _et al._ (2019a).

Since the jet production cross section can be factorized as the convolution of
the hard parton cross section and parton jet function which can be further
expressed as a convolution of jet function in vacuum (in $p$+$p$ collisions)
and jet energy loss distribution, one can extract the jet energy loss
distribution from the experimental data with Bayesian inference He  _et al._
(2019b). In this study we will follow the same procedure of Bayesian inference
in Ref. He  _et al._ (2019b) that extracted the momentum dependence of the jet
energy loss in the most central Pb+Pb collisions at LHC. We carry out a
systematic analysis of world data on the single inclusive jet cross sections
in $p$+$p$ and $A$+$A$ collisions with all possible centrality selections at
both RHIC and LHC energies. We will focus on the system size dependence of the
extracted jet energy loss as well as its scaling behavior with respect to the
charged hadron pseudorapidity density in the final state.

##  II Jet production cross section

We first briefly summarize here the pQCD parton model for jet production that
we use to extract jet energy loss distributions from experimental data. The
differential cross section for single inclusive jet production in $p$+$p$
collisions can be expressed in a factorized and schematical form Kang  _et
al._ (2016, 2017),

$\frac{d\sigma_{pp}^{\rm jet}}{dp_{T}d\eta}=\sum_{a,b,c}\int f_{a/{\rm
p}}\otimes f_{b/{\rm p}}\otimes H_{ab}^{c}\otimes J_{c}(p_{T},R|p_{Tc}),$ (1)

where $f_{a/p}$ is the parton distribution function of proton, $H_{ab}^{c}$ is
the hard function for parton scattering $a+b\rightarrow c+X$ and
$J_{c}(p_{T},R|p_{Tc})$ is the semi-inclusive jet function that describes the
probability for a parent parton $c$ with initial transverse momentum $p_{Tc}$
to produce a jet with transverse energy $p_{T}$ and jet-cone size $R$.
Similarly, the single inclusive jet production cross section in $A$+$A$
collisions can be written as,

$\displaystyle\frac{d\sigma^{\rm jet}_{AA}}{dp_{T}d\eta}$ $\displaystyle=$
$\displaystyle\sum_{a,b,c}\int
d^{2}\mathbf{r}d^{2}\mathbf{b}t_{A}(\mathbf{r})t_{A}(|\mathbf{b}-\mathbf{r}|)\frac{d\phi_{c}}{2\pi}$
(2) $\displaystyle\times f_{a/{A}}\otimes f_{b/{A}}\otimes
H_{ab}^{c}\otimes\widetilde{J}_{c}(p_{T},R,\mathbf{r},\mathbf{b},\phi_{c}|p_{Tc}),$

where $t_{A}(r)$ is the nuclear thickness function with normalization $\int
d^{2}rt_{A}(r)=A$, $f_{a/A}$ is the parton distribution function per nucleon
inside the nucleus $A$, $\mathbf{r}$ is the transverse coordinate of the hard
production vertex of the initial parton ($c$), $\phi_{c}$ is azimuthal angle
between its transverse momentum $p_{Tc}$ and the impact parameter $\mathbf{b}$
of the nucleus-nucleus collision, and
$\widetilde{J}_{c}(p_{T},R,\mathbf{r},\mathbf{b},\phi_{c}|p_{Tc})$ is the
medium-modified semi-inclusive jet function for a given path of jet
propagation in the QGP. The impact parameter $\mathbf{b}$ is integrated over a
range that is determined by the centrality class of the nucleus-nucleus
collisions according to experimental measurements.

The modified jet functions $\widetilde{J}_{c}$ take into account the jet
energy loss due to induced gluon radiation and collisional energy loss carried
by the medium response outside the jet cone. Assuming the in-medium jet
function can be approximated by a shift in the jet energy of the vacuum jet
function and further considering event-by-event fluctuates of jet energy loss
for a given propagation path, the medium-modified jet function can be
expressed as the convolution,

$\displaystyle\widetilde{J}_{c}(p_{T},R,\mathbf{r},\mathbf{b},\phi_{c}|p_{Tc})$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\Delta EJ_{c}(p_{T}+\Delta
E,R|p_{Tc})$ (3) $\displaystyle\times w_{c}(\Delta E,p_{T}+\Delta
E,R,\mathbf{r},\mathbf{b},\phi_{c}),$

of the jet function in vacuum with transverse energy $p_{T}+\Delta E$ and a
jet energy loss distribution $w_{c}$ for a given path specified by
$\mathbf{r}$, $\mathbf{b}$ and $\phi_{c}$. The jet transverse energy loss is
defined as the difference between the jet transverse energy in $p$+$p$ and
$A$+$A$ collisions originating from the same initial parton $c$.

One can define the jet energy loss distribution averaged over the initial
parton production point and propagation direction,

$\displaystyle W^{c}_{AA}(\Delta E,p_{T},R)$ $\displaystyle=$
$\displaystyle\int
d^{2}\mathbf{r}d^{2}\mathbf{b}\frac{t_{A}(r)t_{A}(|\mathbf{b}-\mathbf{r}|)}{N_{\rm
bin}(b)}$ (4) $\displaystyle\times$
$\displaystyle\int\frac{d\phi_{c}}{2\pi}w_{c}(\Delta
E,p_{T},R,\mathbf{r},\mathbf{b},\phi_{c}),~{}~{}$

for a given centrality class of $A$+$A$ collisions, where $N_{\rm bin}(b)=\int
d^{2}\mathbf{r}d^{2}\mathbf{b}t_{A}(r)t_{A}(|\mathbf{b}-\mathbf{r}|)$ is the
number of binary collisions. Then, the single inclusive jet production cross
section in $A$+$A$ collision can be written as

$\displaystyle\frac{d\sigma^{\rm jet}_{AA}}{dp_{T}d\eta}$ $\displaystyle=$
$\displaystyle N_{\rm bin}(b)\sum_{a,b,c}\int d\Delta EW^{c}_{AA}(\Delta
E,p_{T}+\Delta E,R)$ (5) $\displaystyle\times f_{a/{A}}\otimes
f_{b/{A}}\otimes H_{ab}^{c}\otimes J_{c}(p_{T}+\Delta E,R|p_{Tc}),$

In principle, the nuclear modification of the parton distributions (nPDF)
should be considered in the above jet production cross section Adhya  _et al._
(2022). However, the nuclear modification of PDF is usually limited to small
$x$ or very large $x$ (for EMC effect) and low $Q^{2}$ regions Eskola  _et
al._ (2017). For our analyses in this study, we will limit the jet transverse
momentum to $15<p_{T}<30$ GeV/$c$ at RHIC and $50<p_{T}<800$ GeV/$c$ at LHC
which correspond to $0.15<x\approx 2p_{T}/\sqrt{s}<0.3$ and $0.02<x<0.32$ in
the central rapidity region ($y=0$), respectively. For such large
$Q^{2}\approx p_{T}^{2}$, the nuclear modification of the PDF is mostly
negligible in these regions of $x$ Eskola  _et al._ (2017). This is also
confirmed by NLO pQCD parton model calculations Xie  _et al._ (2021) and
indicated by recent experimental data on the single inclusive jet cross
section in minimum-bias p+Pb collisions at LHC Aad  _et al._ (2015a).
Furthermore, the jet cross sections also do not depend on the isospin of the
initial quark flavors.

Under the assumption that jet cross section in the considered kinematic region
is not sensitive to the nuclear modification of parton distribution, the
single inclusive jet cross section in $A$+$A$ collisions can be expressed as
the convolution of jet cross section in $p$+$p$ collisions and a flavor-
averaged (quarks and gluon) jet energy loss distribution $W_{AA}$. The
suppression factor for single inclusive jet production in $A$+$A$ collisions
can be written as

$\displaystyle R_{AA}(p_{T})$ $\displaystyle\approx$
$\displaystyle\frac{1}{d\sigma^{\rm jet}_{pp}(p_{T})}\int d\Delta
Ed\sigma^{\rm jet}_{pp}(p_{T}+\Delta E)$ (6) $\displaystyle\times$
$\displaystyle W_{AA}(\Delta E,p_{T}+\Delta E,R).$

This formula has been used for the study of jet suppression Spousta and Cole
(2016); Mehtar-Tani and Tywoniuk (2018) and similar approximate expression for
single inclusive hadron spectra has been used in Refs. Arleo (2017); Arleo and
Falmagne (2022); Baier  _et al._ (2001) assuming a constant average momentum
fraction of hadrons $z_{h}=p_{Th}/p_{T}$ in the energy loss distribution. This
is also the expression we use to extract jet energy loss distribution from
experimental data through Bayesian inference in this study. In this study, we
will use Pythia8 Sjostrand  _et al._ (2008) simulations to calculate the jet
production cross sections in $p+p$ collisions which are shown to describe the
experimental data very well at both RHIC and LHC energies.

In the above equation the jet energy loss distribution is averaged over parton
flavors (quark versus gluon) weighted by their respective cross sections. It
is possible to extend this analysis to include the flavor dependence with
given fractions of the total jet cross section by the pQCD parton model
calculations Pablos and Soto-Ontoso (2023); Zhang:2023oid; Zhang:2022rby. This
will double the number of parameters in the jet energy loss distributions.
Investigation of such flavor-dependent Bayesian inference and the
corresponding constraints on the momentum and path-length dependence from
single inclusive jet modification is beyond the scope of this study and will
be left for future studies.

##  III Bayesian analysis of single jet suppression $R_{\rm AA}$

Using Eq. (6) that relates the nuclear modification factor $R_{AA}$ of single
inclusive jet cross section and the jet energy loss distribution $W_{AA}$, we
attempt to reverse engineer $W_{AA}$ from experimental data on $R_{AA}$.
Though an exact inversion of the convolution problem can be ill-defined, one
can determine the probability distribution of $W_{AA}$ using Bayesian
inference. This method has found a broad application in heavy-ion collisions,
including the extraction of the QCD equation of state at high-temperature
Pratt  _et al._ (2015), the QGP shear and bulk viscosity Bernhard  _et al._
(2016, 2019); Everett  _et al._ (2021a, b); Nijs  _et al._ (2021a, b);
Parkkila  _et al._ (2021, 2022), the heavy quark diffusion constant Xu  _et
al._ (2018); Ke  _et al._ (2018); Liu  _et al._ (2022), and the jet transport
coefficient in the QGP Ke and Wang (2021); Cao  _et al._ (2021); Liu  _et al._
(2022); Xie  _et al._ (2023, 2022). We will follow Ref. He  _et al._ (2019b)
and carry out a comprehensive Bayesian analysis of the world data on single
inclusive jet spectra in heavy-ion collisions with different centralities at
both RHIC and LHC energies and extract the colliding energy, jet momentum and
system size dependence of the jet energy loss and its fluctuations.

We assume that the average jet energy loss $\langle\Delta E\rangle_{i}$ in a
given colliding system and centrality (collectively labeled by $i$) is given
by

$\displaystyle\langle\Delta E\rangle_{i}=\beta_{i}\left(\frac{p_{T}}{p_{T\rm
ref}}\right)^{\gamma_{i}}\ln\left(\frac{p_{T}}{p_{T\rm ref}}\right),$ (7)

with $p_{T\rm ref}=1$ GeV/$c$ as a reference momentum. $\gamma_{i}>0$
guarantees that the jet energy loss vanishes when jet momentum goes to zero.
The jet energy loss distribution is assumed to only depend on the self-
normalized energy loss fluctuation $x=\Delta E/\langle\Delta E\rangle$, i.e.,
$W_{AA}(\Delta E,p_{T},R)\approx W_{AA}(x,R)$. Such an approximate feature was
corroborated in the previously study using the linear Boltzmann transport
(LBT) model simulations He  _et al._ (2019a). We parametrize the energy loss
fluctuation as

$W_{AA}(x)=\frac{\alpha_{i}^{\alpha_{i}}x^{\alpha_{i}-1}e^{-\alpha_{i}x}}{\Gamma(\alpha_{i})}.$
(8)

From now on, we will omit the label of the jet cone size $R$ dependence of
$W_{AA}$ as we will only use jet measurements with $R=0.4$ for the rest of
this study.

By systematic comparing to dataset $D_{i}$ with $N_{i}$ data points, the
posterior probability distributions of the parameters
$\theta_{i}\equiv[\alpha_{i},\beta_{i},\gamma_{i}]$ are given by the Bayesian
theorem,

$P(\theta_{i}|D_{i})=\frac{P(D_{i}|\theta_{i})P(\theta_{i})}{P(D_{i})},$ (9)

where $P(\theta_{i}|D_{i})$ is the posterior distribution after model-data
comparison. $P(\theta_{i})$ is the prior distribution of model parameters.
$P(D_{i}|\theta_{i})$ is the likelihood function between experimental data and
model $(M)$ calculations using parameter $\theta_{i}$. The likelihood is
assumed to take a Gaussian form,

$\displaystyle\ln P(D_{i}|\theta_{i})$ $\displaystyle=$
$\displaystyle-\frac{N_{i}}{2}\ln(2\pi)-\frac{1}{2}\ln|\Sigma_{i}|$ (10)
$\displaystyle-$
$\displaystyle\frac{1}{2}[M(\theta_{i})-D_{i}]^{T}\Sigma^{-1}[M(\theta_{i})-D_{i}].~{}~{}$

$M(\theta_{i})-D_{i}$ is the discrepancy vector between model calculation and
data. The covariance matrix $\Sigma_{i}$ contains experimental and
computational uncertainty. Uncertainties of the experimental data are assumed
to be uncorrelated between different data points. The normalization
$P(D_{i})=\int d\theta_{i}P(D_{i}|\theta_{i})P(\theta_{i})$ is called the
evidence. The properties of the posterior distribution can be explored with
importance sampling. The common practice is the Markov Chain Monte Carlo
(MCMC) Andrieu  _et al._ (2003) method that performs importance sampling
$\theta_{i}$ according to $P(D_{i}|\theta_{i})P(\theta_{i})$. Here, we use the
affine-invariant MCMC algorithm Goodman and Weare (2010) as implemented in the
emcee package Foreman-Mackey  _et al._ (2013). Projecting the posterior
samples to lower dimensions is equivalent to the marginalization of the high-
dimensional distribution. From the marginalization procedure one can define
the one-parameter posterior distributions and the pairwise correlations.
However, the physically meaningful quantify is not individual parameter but
the functional form of the average energy loss. We can marginalized all three
parameters to obtain the posterior distribution of the averaged energy loss
functional as

$\displaystyle P[\langle\Delta E_{i}\rangle(p_{T})]=\int\langle\Delta
E\rangle(p_{T};\beta_{i},\gamma_{i})P(\theta_{i}|D_{i})d\theta_{i}.$ (11)

Then, one can define the median, and percentile credible interval (C.I.) of
the posterior energy loss at each jet energy.

It is important to note that we do not impose that the parameters
$[\alpha_{i},\beta_{i},\gamma_{i}]$ to be the same for different colliding
systems (colliding energy and centrality). They will be extracted
independently for each experimental data set “$i$”. We will then analyze these
“piecewise” information to study the momentum dependence extracted from
different colliding systems and check whether they are consistent with each
other. Eventually, we will determine the path-length dependence of the jet
energy loss by correlating the extracted jet energy loss with the averaged
path-length for each centrality class of collisions at each colliding energy.

##  IV Extract energy loss parameters

Using the above Bayesian method we have carried out analyses of experimental
data on the nuclear modification of single inclusive jet spectra and extract
the jet energy loss distributions in Pb+Pb collisions at $\sqrt{s}=5.02$ and
2.76 TeV and Au+Au collisions at $\sqrt{s}=200$ GeV for all available
centrality classes. The experimental data on single inclusive jet spectra are
from ATLAS Aad  _et al._ (2015b); Aaboud  _et al._ (2019) for Pb+Pb collisions
at $\sqrt{s}=5.02$ and 2.76 TeV with 8 centrality classes
(0-10%,10-20%,20-30%,30-40%, 40-50%, 50-60%, 60-70% and 70-80%), ALICE Acharya
_et al._ (2020) for Pb+Pb collisions at $\sqrt{s}=5.02$ TeV with 0-10%
centrality class and STAR Adam  _et al._ (2020) for Au+Au collisions at
$\sqrt{s}=200$ GeV with 2 centrality class (0-10%, 60-80%). These experimental
measurements are all at the central rapidity region with jet cone size
$R=0.4$. In Tab. 1, we list these dataset along with the range of transverse
momentum of jet measurements and the value of $\langle N_{\rm
part}\rangle^{1/3}$ which is proportional to the system size or the averaged
path length for a given centrality class. Values of $\langle N_{\rm
part}\rangle$ are obtained from the Glauber model that are used in
experimental analysis Abelev  _et al._ (2013); Adam  _et al._ (2016); Adare
_et al._ (2016). We have not included the jet data from CMS experiment
Khachatryan  _et al._ (2017). Nevertheless, given the relatively larger
experimental uncertainty of CMS data than that from ATLAS in a similar range
of the jet transverse momentum $p_{T}$, we do not expect a significant change
of accuracy for the present Bayesian analysis.

System | Centrality |  $p_{T}$ range [GeV/$c$] | Refs. | $\langle N_{\rm part}\rangle^{1/3}$
---|---|---|---|---
Pb+Pb | 0-10% | $(100,1000)$ | ATLAS Aaboud  _et al._ (2019) | 7.11
5.02 TeV | 10-20% | $(100,630)$ |  | 6.41
| 20-30% | $(79,630)$ |  | 5.73
| 30-40% | $(79,630)$ |  | 5.08
| 40-50% | $(50,398)$ |  | 4.42
| 50-60% | $(50,398)$ |  | 3.77
| 60-70% | $(50,398)$ |  | 3.12
| 70-80% | $(50,251)$ |  | 2.50
| 0-10% | $(60,140)$ | ALICE Acharya  _et al._ (2020) | 7.11
Pb+Pb | 0-10% | $(50,398)$ | ATLAS Aad  _et al._ (2015b) | 7.09
2.76 TeV | 10-20% | $(50,316)$ |  | 6.38
| 20-30% | $(50,316))$ |  | 5.71
| 30-40% | $(39,316)$ |  | 5.05
| 40-50% | $(39,316)$ |  | 4.39
| 50-60% | $(39,316)$ |  | 3.61
| 60-70% | $(39,251)$ |  | 3.10
| 70-80% | $(39,199)$ |  | 2.48
Au+Au | 0-10% | $(15.57,29.14)$ | STAR Adam  _et al._ (2020) | 6.86
200 GeV | 60-80% | $(14,24.52)$ |  | 2.62
Table 1: List of experimental data used in this analysis.  Figure 1:  A
demonstration of calibration to the ATLAS data on single inclusive jet
$R_{AA}$ in 0-10% Pb+Pb collisions at $\sqrt{s}=5.02$ TeV Aaboud  _et al._
(2019). $R_{AA}$ calculated with the prior samples of the parameters
$[\alpha,\beta,\gamma]$ for the jet energy loss distribution are shown as
green lines. The posterior of $R_{AA}$ at 90% credible interval (C.I.) after
the calibration is shown as the red band. The same procedure has been
performed for each dataset listed in Tab. 1.  Figure 2: The posterior
distribution of $[\alpha,\beta,\gamma]$ after calibrating to the ATLAS data on
single inclusive jet $R_{AA}$ in 0-10% Pb+Pb collisions at $\sqrt{s}=5.02$ TeV
Aaboud  _et al._ (2019). The same procedure has been performed for each
dataset listed in Tab. 1 Figure 3: Posterior predictions with median (colored
lines) and 65% credible intervals (shaded bands) compared to the experimental
data Aaboud  _et al._ (2019); Acharya  _et al._ (2020) on single inclusive jet
($R=0.4$) $R_{AA}$ in Pb+Pb collisions at $\sqrt{s}=5.02$ TeV with different
centralities.  Figure 4: The same as Fig. 3 except for Pb+Pb collisions at
$\sqrt{s}=2.76$ TeV with data from Ref.Aad  _et al._ (2015b).  Figure 5: The
same as Fig. 3 except for Au+Au collisions at $\sqrt{s}=200$ GeV with data
from Ref.Adam  _et al._ (2020).  Figure 6: The posterior of the average jet
energy loss ($R=0.4$) with the median (solid lines) with 65% C.I. (solid
bands) as a function of the initial jet transverse momentum $p_{T}$ extracted
from experimental data for Pb+Pb collisions at $\sqrt{s}=5.02$ TeV with
different centralities. Figure 7:  The same as Fig. 6 except for Pb+Pb
collisions at $\sqrt{s}=2.76$ TeV.  Figure 8:  The same as Fig. 6 except for
Au+Au collisions at $\sqrt{s}=200$ GeV.

In Figs. 1 and 2, we show the prior and posterior observables and posterior
distribution of parameters by calibrating to the measured single jet $R_{AA}$
in $0-10$% Pb+Pb collisions at $\sqrt{s}=5.02$ TeV from the ATLAS
Collaboration Aaboud  _et al._ (2019). The prior ranges of parameters are
$0<\alpha<10$, $0<\beta<3.5$, and $0<\gamma<0.3$, which provide a good prior
coverage (green lines) of the experimental data. Calculations using the
posterior distributions of parameters at 90% credible interval (C.I.) (solid
red band) are in good agreement with the experimental data.

Similarly, independent calibrations have been performed for each dataset in
Tab. 1. Shown in Figs. 3, 4 and 5 are the final fits to the single inclusive
jet suppression factor $R_{AA}$ as functions of the final jet transverse
momentum $p_{T}$ in $A$+$A$ collisions with different centrality classes at
three different colliding energies. The solid lines are the mean averages of
the fit and shaded bands are uncertainties at 65% C.I. We observe that errors
at 65% C.I. are consistent with the experimental uncertainties. Note that the
STAR experiment measures charged jets and uses Pythia8 simulations as the
$p+p$ baseline, while ALICE and ATLAS measures full jets with experimental
data as the $p+p$ baseline. Therefore, we have scaled the jet transverse
momentum by 2/3 when Eq. (6) is used to fit the STAR data on charged jets.
Furthermore, STAR charged jets covers relatively low $p_{T}^{\rm jet}$ region,
a leading charged particle trigger $p_{T}^{h^{\pm}}>5$ GeV/$c$ is used to
suppress “fake jet” contributions. However, the leading particle trigger also
biases the jet suppression. Therefore, we will only calibrate to the highest
three $p_{T}^{\rm jet}$ bins where the bias effects are negligible Adam  _et
al._ (2020).

From Fig. 2, one can notice that the posterior distributions of $\beta$ and
$\gamma$ are highly correlated, therefore it is more unambiguous to study the
posterior of the averaged energy loss $\langle\Delta E\rangle$. We plot
$\langle\Delta E\rangle$ at 65% C.I. as a function of the jet transverse
momentum in Figs. 6, 7, and 8 for $A$+$A$ collisions with different centrality
classes at three different colliding energies. The jet energy loss from each
extraction is plotted only within the jet $p_{T}^{\rm jet}$ range of the
dataset. It increases slowly with the jet transverse momentum. It is bigger in
more central collisions and at higher colliding energies, implying that the
jet energy loss increases with the medium density and system size or averaged
propagation length. As we will show in the next section, the extracted
parameter $\alpha$ for the jet energy loss fluctuation in Eq. (8) is
approximately a constant within the accuracy of this analysis, independent of
the colliding energy and centrality within the uncertainties at 65% C.I.

##  V Momentum, density and system size dependence of jet energy loss

The extraction of the averaged jet energy loss and its fluctuation from
experimental data on $R_{AA}$ in the last section is independently performed
for each dataset with the assumed convolution form of the jet suppression
factor $R_{AA}$ in Eq. (6) without specifying the system size and medium
density dependence of the jet energy loss. In this section, we will interpret
the inferred jet energy loss distribution and discuss its dependence on jet
momentum, initial medium density and system size or path-length of the
collision system.

We first focus on the medium density and jet momentum dependence of the
average jet energy loss $\langle\Delta E\rangle$. The averaged jet energy loss
are extracted in different jet momentum range from systems with different beam
energy and centrality. To factor out the system-size dependence from the
discussion, we can first look at the collection of averaged jet energy loss in
heavy-ion collisions at different beam energies but only in 0-10% centrality
class, where the averaged system sizes are expected to be similar. To be more
precise, the averaged path-lengths in $A$+$A$ collisions at RHIC and LHC are
slightly different due to several reasons, including 1) different mass number
of the colliding nuclei 197Au versus 208Pb, 2) larger inelastic nucleon-
nucleon cross-section $\sigma_{NN}^{\rm inel}\approx$70mb at $\sqrt{s}=5.02$
TeV than $\sigma_{NN}^{\rm inel}\approx 42$mb at $\sqrt{s}=200$ GeV, and 3)
longer lifetime of the QGP fireballs at the LHC energies than at RHIC.
However, these are all sub-leading effects and can be estimated through
dynamic simulations. For our discussion here, we assume the average path
length (including both QGP and hot hadronic matter) is proportional to the
average system size, which is related to the average number of participant
nucleons in each centrality class by $R_{\rm eff}\sim N_{\rm part}^{1/3}$. The
estimates of the system size for 0-10% centrality at three beam energies are
in good agreement as shown in Tab. 1.

In order to exam the colliding energy dependence of the jet energy loss
distributions, we first plot in the left panel of Fig. 9 the experimental data
and posterior predictions at 90% C.I. of the jet suppression factor $R_{AA}$
with jet cone size $R=0.4$ as a function of jet $p_{T}$ in 0-10% central Au+Au
at RHIC and Pb+Pb collisions at LHC energies. Across the range of colliding
energies from RHIC and LHC, the suppression factor $R_{AA}$ seems to follow a
common trend in its transverse momentum dependence in the respective kinematic
ranges, increasing with $p_{T}$. The increase also seems to taper off at the
respective upper bound of the kinematic region, corresponding to
$x_{T}=2p_{T}/\sqrt{s}>0.3$.

Figure 9:  (Left) The posterior single inclusive jet suppression factor
$R_{AA}$ with $R=0.4$ at 90% C.I. in 0-10% central Pb+Pb collisions at
$\sqrt{s}=2.76$ and 5.02 TeV and Au+Au collisions at $\sqrt{s}=200$ GeV, as
compared to experimental data Aaboud  _et al._ (2019); Acharya  _et al._
(2020); Aaboud  _et al._ (2019); Adam  _et al._ (2020). (Right) The
corresponding posterior jet energy loss scaled by the charge multiplicity
density $(dN_{\rm ch}(\sqrt{s})/d\eta)/(dN_{\rm ch}(2.76{\rm TeV})/d\eta)$.
The black solid line $\sim p_{T}^{0.13}\ln p_{T}$ is a fit to the scaled
posterior jet energy loss from all systems shown in this figure.

Since the parton energy loss is proportional to the initial value of the jet
transport coefficient or initial parton density Baier  _et al._ (1997);
Zakharov (1996); Wiedemann (2000); Gyulassy  _et al._ (2000); Guo and Wang
(2000); Wang and Guo (2001); Majumder (2012) which in turn is proportional to
the rapidity density $dN_{\rm ch}/d\eta$ of final charged multiplicity, we
factor out such colliding energy loss dependence when compare the extracted
jet energy loss from experimental data at different colliding energies. Shown
in the right panel of Fig. 9 is the extracted average jet energy loss scaled
by the charged multiplicity density $dN_{\rm ch}(\sqrt{s})/dN_{\rm
ch}(2.76{\rm TeV})$ in central 0-10% Au+Au collisions at $\sqrt{s}=200$ GeV
and Pb+Pb collisions at $\sqrt{s}=2.76$ and 5.02 TeV Adare  _et al._ (2016);
Aamodt  _et al._ (2011b); Adam  _et al._ (2016). The colored solid lines are
the median energy loss and shaded bands are the 90% C.I. We can see the scaled
jet energy loss has a common momentum dependence which increases a little
faster than a simple logarithmic dependence.

Since the initial parton density can be estimated as proportional to
$\rho\sim(dN_{\rm ch}/d\eta)/\pi R_{\rm eff}^{2}$ and the effective system
size can be related to the average number of participant nucleons for a given
centrality class $R_{\rm eff}\sim\langle N_{\rm part}\rangle^{1/3}$, we can
parameterize the jet energy loss as,

$\langle\Delta E\rangle=\frac{dN_{\rm ch}(\sqrt{s})}{d\eta}\frac{1}{\langle
N_{\rm part}\rangle^{2/3}}f\left(\langle N_{\rm
part}\rangle^{1/3}\right)g(p_{T}).$ (12)

where $f\left(\langle N_{\rm part}\rangle^{1/3}\right)$ and $g(p_{T})$ are now
assumed to be universal functions for all collision systems and centralities.
Fitting to the momentum dependence of the extracted jet energy loss in the
right panel of Fig. 9 for central 0-10% Au+Au and Pb+Pb collisions, we find
$g(p_{T})\approx p_{T}^{0.13}\ln p_{T}$ which is slightly stronger than a
logarithmic dependence. This is consistent with the earlier Bayesian analysis
He  _et al._ (2019b). The $p_{T}$ ranges of available experimental data on
single inclusive jet spectra with other centralities are smaller than the most
central collisions in Fig. 9. The extracted $p_{T}$ dependence of the jet
energy loss for these semi-central and semi-peripheral collisions is
consistent with the above functional form.

To find out the system size dependence of the extracted jet energy loss, we
plot in the left panel of Fig. 10 the scaled jet energy loss $\langle N_{\rm
part}\rangle^{2/3}\langle\Delta E\rangle/(dN_{\rm ch}/d\eta)$ as a function of
$\langle N_{\rm part}\rangle^{1/3}$ for final jet $p_{T}=60-100$ GeV/$c$ in
Au+Au collisions at RHIC and Pb+Pb collisions at LHC energies with different
centrality classes. The values of $\langle N_{\rm part}\rangle^{1/3}$ from
Glauber model are listed in Tab. 1. Excluding the extracted values in very
(60-80%) peripheral collisions where the Bayesian fitting does not do well,
the extracted jet energy loss has an approximate $f(\langle N_{\rm
part}\rangle^{1/3})\sim(\langle N_{\rm part}\rangle^{1/3})^{0.59}$ dependence
on the system size. This is quite different from the approximate linear
dependence on the system size for the energy loss of a single parton in an
expanding system Gyulassy  _et al._ (2002), which is needed to explain system
size dependence of the suppression of single inclusive hadron spectra Arleo
and Falmagne (2022). Jet energy loss, however, is not proportional to energy
loss of individual shower partons. Energy carried by radiated gluons and
medium response within the jet cone is recovered by the jet, reducing the jet
energy loss. This will both lead to the weaker system size dependence as shown
by the LBT model simulations He  _et al._ (2019a, b).

We also extract the parameters $\alpha$ for the jet energy loss fluctuations
in Eq. (8) and they are shown in the right panel of Fig. 10 as a function of
$\langle N_{\rm part}\rangle^{1/3}$. Within the larger errors at 65% C.I., the
jet energy loss fluctuation parameter $\alpha\approx 5.0\pm 3$ does not depend
on the system size and the initial parton density (or colliding energy).

Figure 10: The posterior (Left) average jet energy loss scaled by the initial
medium density $(dN_{\rm ch}/d\eta)/\langle N_{\rm part}\rangle^{2/3}$ and
(Right) the parameter $\alpha$ for jet energy loss fluctuations as a function
of $\langle N_{\rm part}\rangle^{1/3}.$ The solid line is a power-law fit
$(\langle N_{\rm part}\rangle^{1/3})^{0.557}$ excluding the peripheral
collisions with $\langle N_{\rm part}\rangle^{1/3}<3.3$.

##  VI Summary

Employing Bayesian inference, we have analyzed the world data on the nuclear
suppression of single inclusive jet spectra in heavy-ion collisions at both
RHIC and LHC energies with a wide selection of centrality classes. We have
extracted the average jet energy loss as a function of the jet transverse
momentum for each centrality class of collisions at each colliding energy and
the jet energy loss fluctuations. We found that the extracted average jet
energy loss scales with the initial parton density $\rho\propto(dN_{\rm
ch}/d\eta)/\langle N_{\rm part}\rangle^{2/3}$. The average jet energy loss
scaled by the initial parton density has a jet momentum dependence
$\langle\Delta E\rangle/\rho\propto p_{T}^{0.13}\ln p_{T}$ that is slightly
stronger than a logarithmic form. It has a system size dependence
$\langle\Delta E\rangle/\rho\propto(\langle N_{\rm
part}\rangle^{1/3})^{0.59}$. This behavior of jet energy loss is different
from that of a single energetic parton. Such a difference can be attributed to
energy loss carried by the radiated gluons and medium response outside the jet
cone.

## Acknowledgements

This work is supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics under Contracts No. DE-AC02-05CH11231 and No.
89233218CNA000001, and by NSF under Grant No. OAC-2004571 within the X-SCAPE
Collaboration. W.K. is also supported by the Laboratory Directed Research and
Development Program at LANL.

## References

  * Gyulassy and Plumer (1990) M. Gyulassy and M. Plumer, Phys. Lett. B 243, 432 (1990).
  * Wang and Gyulassy (1992) X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992).
  * Adcox  _et al._ (2002) K. Adcox _et al._ (PHENIX), Phys. Rev. Lett. 88, 022301 (2002), arXiv:nucl-ex/0109003 .
  * Adler  _et al._ (2002) C. Adler _et al._ (STAR), Phys. Rev. Lett. 89, 202301 (2002), arXiv:nucl-ex/0206011 .
  * Wang (2005) X.-N. Wang, Nucl. Phys. A 750, 98 (2005), arXiv:nucl-th/0405017 .
  * Aad  _et al._ (2010) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 105, 252303 (2010), arXiv:1011.6182 [hep-ex] .
  * Aamodt  _et al._ (2011a) K. Aamodt _et al._ (ALICE), Phys. Lett. B 696, 30 (2011a), arXiv:1012.1004 [nucl-ex] .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Eur. Phys. J. C 72, 1945 (2012), arXiv:1202.2554 [nucl-ex] .
  * Khachatryan  _et al._ (2017) V. Khachatryan _et al._ (CMS), Phys. Rev. C 96, 015202 (2017), arXiv:1609.05383 [nucl-ex] .
  * Acharya  _et al._ (2020) S. Acharya _et al._ (ALICE), Phys. Rev. C 101, 034911 (2020), arXiv:1909.09718 [nucl-ex] .
  * Chen  _et al._ (2010) X.-F. Chen, C. Greiner, E. Wang, X.-N. Wang, and Z. Xu, Phys. Rev. C 81, 064908 (2010), arXiv:1002.1165 [nucl-th] .
  * Burke  _et al._ (2014) K. M. Burke _et al._ (JET), Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003 [nucl-th] .
  * Cao  _et al._ (2021) S. Cao _et al._ (JETSCAPE), Phys. Rev. C 104, 024905 (2021), arXiv:2102.11337 [nucl-th] .
  * Liu  _et al._ (2022) F.-L. Liu, W.-J. Xing, X.-Y. Wu, G.-Y. Qin, S. Cao, and X.-N. Wang, Eur. Phys. J. C 82, 350 (2022), arXiv:2107.11713 [hep-ph] .
  * Apolinário  _et al._ (2022) L. Apolinário, Y.-J. Lee, and M. Winn, Prog. Part. Nucl. Phys. 127, 103990 (2022), arXiv:2203.16352 [hep-ph] .
  * Xie  _et al._ (2023) M. Xie, W. Ke, H. Zhang, and X.-N. Wang, Phys. Rev. C 108, L011901 (2023), arXiv:2206.01340 [hep-ph] .
  * Xie  _et al._ (2022) M. Xie, W. Ke, H. Zhang, and X.-N. Wang, (2022), arXiv:2208.14419 [hep-ph] .
  * Baier  _et al._ (1997) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, Nucl. Phys. B 483, 291 (1997), arXiv:hep-ph/9607355 .
  * Zakharov (1996) B. G. Zakharov, JETP Lett. 63, 952 (1996), arXiv:hep-ph/9607440 .
  * Wiedemann (2000) U. A. Wiedemann, Nucl. Phys. B 588, 303 (2000), arXiv:hep-ph/0005129 .
  * Gyulassy  _et al._ (2000) M. Gyulassy, P. Levai, and I. Vitev, Phys. Rev. Lett. 85, 5535 (2000), arXiv:nucl-th/0005032 .
  * Guo and Wang (2000) X.-F. Guo and X.-N. Wang, Phys. Rev. Lett. 85, 3591 (2000), arXiv:hep-ph/0005044 .
  * Wang and Guo (2001) X.-N. Wang and X.-F. Guo, Nucl. Phys. A 696, 788 (2001), arXiv:hep-ph/0102230 .
  * Majumder (2012) A. Majumder, Phys. Rev. D 85, 014023 (2012), arXiv:0912.2987 [nucl-th] .
  * Gyulassy  _et al._ (2002) M. Gyulassy, I. Vitev, X.-N. Wang, and P. Huovinen, Phys. Lett. B 526, 301 (2002), arXiv:nucl-th/0109063 .
  * Arleo and Falmagne (2022) F. Arleo and G. Falmagne, (2022), arXiv:2212.01324 [hep-ph] .
  * He  _et al._ (2015) Y. He, T. Luo, X.-N. Wang, and Y. Zhu, Phys. Rev. C 91, 054908 (2015), [Erratum: Phys.Rev.C 97, 019902 (2018)], arXiv:1503.03313 [nucl-th] .
  * Wang  _et al._ (2017) X.-N. Wang, S.-Y. Wei, and H.-Z. Zhang, Phys. Rev. C 96, 034903 (2017), arXiv:1611.07211 [hep-ph] .
  * He  _et al._ (2019a) Y. He, S. Cao, W. Chen, T. Luo, L.-G. Pang, and X.-N. Wang, Phys. Rev. C 99, 054911 (2019a), arXiv:1809.02525 [nucl-th] .
  * He  _et al._ (2019b) Y. He, L.-G. Pang, and X.-N. Wang, Phys. Rev. Lett. 122, 252302 (2019b), arXiv:1808.05310 [hep-ph] .
  * Kang  _et al._ (2016) Z.-B. Kang, F. Ringer, and I. Vitev, JHEP 10, 125 (2016), arXiv:1606.06732 [hep-ph] .
  * Kang  _et al._ (2017) Z.-B. Kang, F. Ringer, and I. Vitev, Phys. Lett. B 769, 242 (2017), arXiv:1701.05839 [hep-ph] .
  * Adhya  _et al._ (2022) S. P. Adhya, C. A. Salgado, M. Spousta, and K. Tywoniuk, Eur. Phys. J. C 82, 20 (2022), arXiv:2106.02592 [hep-ph] .
  * Eskola  _et al._ (2017) K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, Eur. Phys. J. C 77, 163 (2017), arXiv:1612.05741 [hep-ph] .
  * Xie  _et al._ (2021) M. Xie, X.-N. Wang, and H.-Z. Zhang, Phys. Rev. C 103, 034911 (2021), arXiv:2003.02441 [hep-ph] .
  * Aad  _et al._ (2015a) G. Aad _et al._ (ATLAS), Phys. Lett. B 748, 392 (2015a), arXiv:1412.4092 [hep-ex] .
  * Spousta and Cole (2016) M. Spousta and B. Cole, Eur. Phys. J. C 76, 50 (2016), arXiv:1504.05169 [hep-ph] .
  * Mehtar-Tani and Tywoniuk (2018) Y. Mehtar-Tani and K. Tywoniuk, Phys. Rev. D 98, 051501 (2018), arXiv:1707.07361 [hep-ph] .
  * Arleo (2017) F. Arleo, Phys. Rev. Lett. 119, 062302 (2017), arXiv:1703.10852 [hep-ph] .
  * Baier  _et al._ (2001) R. Baier, Y. L. Dokshitzer, A. H. Mueller, and D. Schiff, JHEP 09, 033 (2001), arXiv:hep-ph/0106347 .
  * Sjostrand  _et al._ (2008) T. Sjostrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * Pablos and Soto-Ontoso (2023) D. Pablos and A. Soto-Ontoso, Phys. Rev. D 107, 094003 (2023), arXiv:2210.07901 [hep-ph] .
  * Pratt  _et al._ (2015) S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Phys. Rev. Lett. 114, 202301 (2015), arXiv:1501.04042 [nucl-th] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Bernhard  _et al._ (2019) J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nature Phys. 15, 1113 (2019).
  * Everett  _et al._ (2021a) D. Everett _et al._ (JETSCAPE), Phys. Rev. C 103, 054904 (2021a), arXiv:2011.01430 [hep-ph] .
  * Everett  _et al._ (2021b) D. Everett _et al._ (JETSCAPE), Phys. Rev. Lett. 126, 242301 (2021b), arXiv:2010.03928 [hep-ph] .
  * Nijs  _et al._ (2021a) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. C 103, 054909 (2021a), arXiv:2010.15134 [nucl-th] .
  * Nijs  _et al._ (2021b) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. Lett. 126, 202301 (2021b), arXiv:2010.15130 [nucl-th] .
  * Parkkila  _et al._ (2021) J. E. Parkkila, A. Onnerstad, and D. J. Kim, Phys. Rev. C 104, 054904 (2021), arXiv:2106.05019 [hep-ph] .
  * Parkkila  _et al._ (2022) J. E. Parkkila, A. Onnerstad, S. F. Taghavi, C. Mordasini, A. Bilandzic, M. Virta, and D. J. Kim, Phys. Lett. B 835, 137485 (2022), arXiv:2111.08145 [hep-ph] .
  * Xu  _et al._ (2018) Y. Xu, J. E. Bernhard, S. A. Bass, M. Nahrgang, and S. Cao, Phys. Rev. C 97, 014907 (2018), arXiv:1710.00807 [nucl-th] .
  * Ke  _et al._ (2018) W. Ke, Y. Xu, and S. A. Bass, Phys. Rev. C 98, 064901 (2018), arXiv:1806.08848 [nucl-th] .
  * Ke and Wang (2021) W. Ke and X.-N. Wang, JHEP 05, 041 (2021), arXiv:2010.13680 [hep-ph] .
  * Andrieu  _et al._ (2003) C. Andrieu, N. de Freitas, A. Doucet, and M. I. Jordan, Machine Learning 50, 5 (2003).
  * Goodman and Weare (2010) J. Goodman and J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Aad  _et al._ (2015b) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 114, 072302 (2015b), arXiv:1411.2357 [hep-ex] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Lett. B 790, 108 (2019), arXiv:1805.05635 [nucl-ex] .
  * Adam  _et al._ (2020) J. Adam _et al._ (STAR), Phys. Rev. C 102, 054913 (2020), arXiv:2006.00582 [nucl-ex] .
  * Abelev  _et al._ (2013) B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044909 (2013), arXiv:1301.4361 [nucl-ex] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), Phys. Rev. Lett. 116, 222302 (2016), arXiv:1512.06104 [nucl-ex] .
  * Adare  _et al._ (2016) A. Adare _et al._ (PHENIX), Phys. Rev. C 93, 024901 (2016), arXiv:1509.06727 [nucl-ex] .
  * Aamodt  _et al._ (2011b) K. Aamodt _et al._ (ALICE), Phys. Rev. Lett. 106, 032301 (2011b), arXiv:1012.1657 [nucl-ex] .