# Semantic Workflows and Machine Learning for the Assessment of Carbon Storage
by Urban Trees

Juan Carrillo  University of WaterlooWaterlooCanada<EMAIL_ADDRESS>,
Daniel Garijo  Information Sciences InstituteUniversity of Southern
CaliforniaUSA<EMAIL_ADDRESS>,  Mark Crowley  University of
WaterlooWaterlooCanada<EMAIL_ADDRESS>,  Rober Carrillo  United
World CollegesBogotaColombia<EMAIL_ADDRESS>,  Yolanda Gil
Information Sciences InstituteUniversity of Southern CaliforniaUSA<EMAIL_ADDRESS>and  Katherine Borda  WaterlooCanada<EMAIL_ADDRESS>
###### Abstract.

Climate science is critical for understanding both the causes and consequences
of changes in global temperatures and has become imperative for decisive
policy-making. However, climate science studies commonly require addressing
complex interoperability issues between data, software, and experimental
approaches from multiple fields. Scientific workflow systems provide
unparalleled advantages to address these issues, including reproducibility of
experiments, provenance capture, software reusability and knowledge sharing.
In this paper, we introduce a novel workflow with a series of connected
components to perform spatial data preparation, classification of satellite
imagery with machine learning algorithms, and assessment of carbon stored by
urban trees. To the best of our knowledge, this is the first study that
estimates carbon storage for a region in Africa following the guidelines from
the Intergovernmental Panel on Climate Change (IPCC).

Reproducibility, scientific workflows, machine learning, land cover mapping,
carbon assessment, Sentinel-2

Copyright ©2019 for this paper by its authors. Use permitted under Creative
Commons License Attribution 4.0 International (CC BY 4.0).

††conference: Sciknow ’19: Third International Workshop on Capturing
Scientific Knowledge; November 19, 2019; Los Angeles, California, USA

##  1\. Introduction

Climate science requires modeling natural and man-made processes that are
highly complex, exhibit non-linear dynamics and possess disparate spatial and
temporal scales. Handling this complexity requires a holistic approach among
multiple disciplines (Meadow et al., 2015), but scientists from different
fields may also need to use domain-specific data sources, methods, and
computational models. The integration of their knowledge and experiments is a
challenging task (Daron et al., 2015), especially when a study is expected to
provide actionable insights for decision making at regional and local scale.

Scientific workflows have emerged as an integrated solution to manage this
challenge, as they capture the computational steps and data dependencies
required to carry out a computational experiment (Taylor et al., 2014).
Scientific workflows ease data handling (metadata, provenance), component
versioning (parametrization, calibration) and have a clear separation between
workflow design and workflow execution (Barseghian et al., 2010) (Wolstencroft
et al., 2013). One of the major advantages of scientific workflow systems is
their role in improving the reproducibility of scientific studies.
Reproducibility plays a critical role in climate sciences due to their impact
in our society (Schulson, 2018). In fact, due to issues with the documentation
of experiments, some climate science studies have been re-examined lately due
to their impact in global policy-making (Gillingham et al., 2018) and water
resources management (Kundzewicz et al., 2018).

In this paper we describe the process we followed to design and implement
reusable scientific workflows for the climate sciences. In particular, we
focus on carbon storage assessment by using urban trees, a common requirement
for cities to reduce their carbon emissions globally. Our contributions
include the development of a library of components for preparing geospatial
data by doing coordinate transformations, the integration of machine learning
components to classify trees in satellite images and the creation of workflows
for carbon storage assessment in cities. In order to implement these
workflows, we use the Workflow Instance Generation and Selection (WINGS)
system (Gil et al., 2011), which has been successfully used for applications
in domains ranging from Genomics (Gil et al., 2012b) to Geosciences (Gil et
al., 2012a).

The paper starts by giving an overview of previous research on the assessment
of carbon storage by urban trees according to the guidelines published by the
Intergovernmental Panel on Climate Change (IPCC), highlighting the advantages
and limitations of the most recent methods. We then describe the design
considerations of our scientific workflows as well as their experimental
implementation and evaluation. The paper continues with a discussion of
results and suggested future work.

##  2\. Background: Carbon emissions and storage

There is an increasing interest among scientific organizations and national
governments regarding Carbon emissions and their role in climate science
(Figueres et al., 2017). The consequences of higher concentrations of carbon
gases in the atmosphere are now more explicit and international organizations
such as the Intergovernmental Panel on Climate Change (IPCC) are leading
initiatives to monitor national efforts to lower emissions and increase carbon
storage (Edenhofer, 2015). Monitoring carbon emissions is fundamental to
inform government policies in topics such as renewable energy, transportation,
and manufacturing technologies. Similarly, the assessment of carbon storage is
equally important to guide hazard mitigation efforts (Gough, 2016).

Most studies in climate science require domain knowledge to design and run the
experiments (Incropera, 2016). But the increasing concerns about the changing
climate require strategies to streamline the replication of assessment
studies, the reusability of data, methods, and results. The use of scientific
workflows can significantly improve the implementation and reproducibility of
carbon assessment studies, with additional gains in data and model sharing as
well as knowledge capture through semantic representations.

###  2.1. Carbon assessment by urban trees

Urban trees provide a natural and cost-effective alternative to capture and
store carbon in cities. Having trees in densely populated areas also improve
human health and biodiversity and provide benefits for flood prevention and
reduced cooling costs, among other benefits (Livesley et al., 2016). In 2003,
the IPCC published the Good Practice Guidance for Land Use, Land-Use Change
and Forestry (Penman et al., 2003) and in 2006 the IPCC Guidelines for
National Greenhouse Gas Inventories (Eggleston et al., 2006). These guidelines
suggest the use of area covered by trees, shrubs, and herbaceous (perennial)
plants to determine the amount of carbon stored as biomass in settlements.
However, due to the limited availability of detailed data the majority of
published studies focus only on tree cover. While these two documents describe
the stages of an assessment study, including aspects such as data collection
and uncertainty estimation, they suggest governments to deal with minor
implementation details according to their technical capacity and available
resources.

Published carbon assessment studies have used different combinations of data,
methods, and software. These studies can be divided into two major groups
according to the data collection approaches and models they use. Assessments
in the first group use a statistics point sampling technique to estimate tree
density from aerial imagery (McGovern and Pasher, 2016)(Pasher et al.,
2014)(Nowak et al., 2013). This method is easy to implement and only requires
imagery for sample areas, but the outcome is a percentage value that does not
describe the spatial distribution of trees. The second group of methods use
LiDAR, aerial or satellite imagery to provide a comprehensive assessment of
tree coverage, including their spatial distribution (Tigges et al.,
2017)(Schreyer et al., 2014)(Davies et al., 2011)(Raciti et al., 2014).
However, the second approach requires imagery for the complete the area of
study as well as the configuration of more complex methods for detection or
classification of urban trees. Later in this document we introduce our own
method, which belongs to the second group and uses freely available satellite
imagery and a carefully designed workflow to facilitate implementation and
reuse.

One common limitation of carbon assessment studies is the lack of a systematic
approach to share data, models, software, and results. Regardless of the
specific data source, technology, or processing method, most reports only
contain descriptions of the work done, which are not enough to replicate the
experiments or reuse the software tools. It is difficult for other scientists,
especially in developing countries, to translate the information from those
reports into actionable knowledge that serves them to efficiently incorporate
published research into their own studies.

##  3\. Scientific Workflow Design

With the advantages and limitations of published methods in mind, we design a
new workflow to efficiently determine tree cover for urban areas We start by
presenting the advantages of knowledge capture systems to represent models in
geosciences as scientific workflows and then describe how we leverage previous
research on carbon assessment and tree mapping to design our own workflow.

The workflow is designed as multiple interconnected components in WINGS that
operate in three consecutive stages as seen in Figure 1, data preprocessing,
mapping of tree coverage, and assessment of carbon storage. Our workflow is
based on previous work in which high resolution satellite imagery are used to
produce land cover maps over urban areas (Raciti et al., 2014) (Schreyer et
al., 2014). However, we focus on using freely available medium resolution
satellite imagery from the Sentinel-2 sensors (Drusch et al., 2012) to
facilitate replication by other researchers.

Figure 1. Stages of our carbon assessment workflow

Stages of our carbon assessment workflow

Spatial data preprocessing stage involves common operations for experiments
across earth sciences, such as transformation of coordinate systems and
conversion between file formats. We have designed all these data preprocessing
steps as reusable building blocks so they can be included in other workflows.
These preprocessing operations are also known as Extract Transform Load (ETL)
tasks and are implemented using the Geospatial Data Abstraction Library GDAL
(GDAL/OGR contributors, 2019).

We map the tree coverage using satellite image classification and design
multiple components to train machine learning algorithms, classify the image
over an area of interest, produce a visualization ready tree cover map, and
determine the resulting accuracy. The Machine Learning algorithms we implement
and calibrate are Random Forests and Support Vector Machines, both available
in the Orfeo Remote Sensing toolbox (Grizonnet et al., 2017). We train the
algorithms using sample points collected through visual inspection.

Our workflow generates a map that includes other land cover categories such as
water, grass and built areas; which may additionally serve for studies in
hydrology, planning, and other applications (as seen in Figure 2). Moreover,
this map is generated in a standard format for further use in Geographic
Information Systems or other scientific platforms.

Figure 2. Some applications of land cover mapping

Some applications of land cover mapping

The assessment of carbon storage is completed following the IPCC guidelines to
calculate carbon stored based on urban canopy area. In the calculation we
multiply the canopy area by a conversion factor to estimate carbon stored in
the form of biomass. Since no values are published specifically for Africa
(our initial region of interest) we use a default value suggested by the IPCC.

Figure 3. City of Juba in South Sudan, Africa

Juba_in_Africa

##  4\. Area of Study and Data

Our area of study is the city of Juba in South Sudan. It is the current
capital of the country and serves as its main commercial and transportation
hub with an estimated population of nearly 386,000 inhabitants (Agency, 2018).
Juba is located in the southern region of the country and has an extension of
103 km2 according to the urban boundary retrieved in July 2019 from Open
Street Maps (Haklay and Weber, 2008). Figure 3 shows the location of Juba and
South Sudan in the African continent.

The country of South Sudan currently faces multiple issues, including
political instability (Sørbø and Ahmed, 2013), poor health services (Macharia
et al., 2017), and a lack of infrastructure, especially for storage and
distribution of water (Nassif et al., 2016). Few mapping projects have been
conducted in Africa, and especially in South Sudan (Baariu et al.,
2019)(Kagawa and Le Sourda, 2018), which hinder the work of humanitarian and
non government organizations; therefore, any study aiming to understand the
characteristics of the rural or urban territories will likely produce positive
outcomes in the short and long term.

The data we use includes 4000 sample points digitized through visual
inspection from Google Earth high resolution satellite imagery. The sample
points are spatially dispersed across the city and distributed in four
categories, tree canopy, grass, impervious, and water, with 1000 points per
category to have a balanced dataset. Impervious includes all built or bare
areas not covered by vegetation, water, or agriculture. Examples of impervious
areas are constructions, parking lots, roads, airport runways, and rocky
surfaces.

Additionally, we use a multi-spectral satellite image from the Sentinel-2
sensor made available by the European Space Agency (ESA) (Drusch et al.,
2012). While this image originally includes 13 bands we use only the four
bands with a pixel size equivalent to 10 meters on the ground. The first three
bands collect data in the visible light spectrum and the fourth in the near
infrared spectrum, with the latter being particularly useful for detecting
vegetation.

Figure 4 shows some sample points of the four land cover categories in a
neighborhood close to the White Nile river in the city of Juba, overlaid on
high resolution satellite imagery from Google Earth and the same area as seen
on the medium resolution Sentinel-2 image.

Figure 4. Sample points and satellite imagery of Juba

Sample points and satellite imagery of Juba

##  5\. Implementation and Results

The workflow is implemented as a series of software components in WINGS.
Additionally, each component is designed to run a particular task and
configured with relevant parameters to calibrate their functionality.
Components are coded as Python scripts that serve as a high level interface
for us to use libraries such as GDAL for data preprocessing and Orfeo command
line tools for satellite image classification.

These components are connected through intermediate datasets that in turn are
outputs and inputs for the previous and following components. We created 14
components in total, with eight dedicated for data preparation, five for
mapping of tree coverage, and one for carbon assessment. The use of components
as modular pieces of software to accomplish specific data processing tasks
creates opportunities for reusability across a variety of models in Earth
Sciences and Geospatial technologies.

The eight components designed for data preparation allow researchers to handle
datasets in the most common file formats and transform them according to the
particular goals of the study. Some operations correspond to sub-setting the
spatial extent of the data, changing its coordinate reference system, and
editing the attributes for tabular datasets. When running the workflow for a
new area of study, the task of setting parameters such as the coordinate
reference system is facilitated by the WINGS system, which suggest the most
appropriate value according to the geographic region as configured by the
workflow designer. Figure 5 shows a fragment of our workflow where we use
multiple data preparation components to perform a format transformation and
reprojecting a file.

Figure 5. Workflow fragment for data preparation of sample points

Workflow fragment for data preparation of sample points

Mapping of tree coverage includes four components focused on training the
Machine Learning algorithms to classify the Sentinel-2 satellite image using
the sample points digitized through visual inspection. Initially, one
component extracts the pixel values of the satellite image for the 4000 sample
locations. Next, another component uses 80% of these values as training and
20% as validation data to train the Random Forest and Support Vector Machine
algorithms. Then we use the trained algorithms to perform pixel level
classification of the Sentinel-2 image for the entire area of interest and
output the tree cover map. Subsequently, one more component evaluates the
classification accuracy. Figure 6 shows a fragment of our workflow where we
use multiple components to map tree cover.

Figure 6. Workflow fragment for mapping tree cover

Workflow fragment for mapping tree cover

Tables 1 and 2 show the confusion matrix with and without normalization
resulting from the evaluation of the Random Forest classification algorithm.
We see that the tree cover category is the one with the lowest accuracy with
only 108 out of the 200 trees in the test set correctly classified. This is
likely the result of using medium resolution satellite imagery (Sentinel-2)
with a ground pixel size of 10 m, in other words, the canopy area of a tree
should be about 100 m2 to be easily identifiable in at least one pixel,
without considering boundary issues between adjacent pixels. The grass and
impervious land cover categories exhibit a comparable accuracy of 65% and 73%
consequently. While areas corresponding to these two categories show a
slightly better accuracy they are still hard to differentiate, presumably due
to grass patches and house rooftops with a size smaller than the area of a
pixel (100 m2). For the water land cover category the algorithm reaches an
almost perfect accuracy, which is anticipated due to the significant
difference in the way it reflects the light compared to the other three land
cover categories, especially in the near infrared band (B4).

Table 1. Normalized confusion matrix True label | Trees | 0.54 | 0.18 | 0.28 | 0
---|---|---|---|---|---
Grass | 0.32 | 0.65 | 0.14 | 0
Impervious | 0.17 | 0.09 | 0.73 | 0
Water | 0.01 | 0 | 0 | 0.99
| Trees | Grass | Impervious | Water
Predicted label
Table 2. Confusion matrix, without normalization True label | Trees | 108 | 37 | 55 | 0
---|---|---|---|---|---
Grass | 42 | 131 | 27 | 0
Impervious | 35 | 17 | 147 | 1
Water | 1 | 0 | 0 | 199
| Trees | Grass | Impervious | Water
Predicted label

The assessment of carbon storage is performed by one component that receives
the tree cover map shown in Figure 7, calculates the total canopy area and
uses a carbon removal factor to determine the total amount of carbon stored as
biomass in the city of Juba. We use a default carbon removal factor of 2.9
tonnes C (ha crown cover)-1 yr -1, which is the value suggested by the IPCC
for Tier 2a studies (Eggleston et al., 2006). Multiplying the default carbon
removal factor by the total tree cover area (10,519 ha) from the previous
stage, we calculate that trees in the city of Juba remove 30,506 tonnes C yr
-1. This amount is equivalent to the carbon dioxide emitted by 6632 passenger
vehicles per year (EPA, 2019). This value along with the classification
accuracy from the previous stage are valuable information for countries to
prepare their carbon assessment reports according to the IPCC guidelines. The
code used to create the geospatial transformations
111https://github.com/jmcarrillog/geospatial-etl and the carbon assessment
workflows 222https://github.com/jmcarrillog/machine-learning-workflow-for-
carbon-assessment is available online.

Figure 7. Resulting land cover map for the city of Juba

Resulting land cover map for the city of Juba

##  6\. Conclusions and Future Work

In this paper we introduced our work to create a library of workflow
components to perform spatial data transformations, land cover mapping and
assessment of carbon storage. By leveraging scientific workflows, we aim to
ease the reusability of these components in other workflows and the
reproducibility and transparency of carbon assessment studies. Our future work
will focus on two main areas: first, we aim to test our workflow using data
from other locations around the globe, which requires additional training data
points. Second, we will focus on calibration of the parameters for the
classifiers to improve our classification accuracy during the land cover
mapping stage.

##  7\. Acknowledgments

Juan Carrillo was financially supported for this project by the Mitacs
Globalink Research Award and the University of Waterloo Machine Learning Lab.
Juan Carrillo gives special thanks to Mark Crowley, Daniel Garijo, and Yolanda
Gil for their mentoring and valuable contributions during this research
project.

## References

  * (1)
  * Agency (2018) Central Intelligence Agency. 2018\.  The World Factbook: South Sudan.  https://www.cia.gov/library/publications/the-world-factbook/geos/od.html
  * Baariu et al. (2019) Sabina N Baariu, Galcano C Mulaku, David N Siriba, et al. 2019\.  State of Cartographic Services among the East African Community Member States.  _Journal of Geographic Information System_ 11, 01 (2019), 56\.
  * Barseghian et al. (2010) Derik Barseghian, Ilkay Altintas, Matthew B. Jones, Daniel Crawl, Nathan Potter, James Gallagher, Peter Cornillon, Mark Schildhauer, Elizabeth T. Borer, Eric W. Seabloom, and Parviez R. Hosseini. 2010\.  Workflows and extensions to the Kepler scientific workflow system to support environmental sensor data access and analysis.  _Ecological Informatics_ 5, 1 (2010), 42 – 50.  https://doi.org/10.1016/j.ecoinf.2009.08.008 Special Issue: Advances in environmental information management.
  * Daron et al. (2015) Joseph David Daron, Kate Sutherland, Christopher Jack, and Bruce C. Hewitson. 2015.  The role of regional climate projections in managing complex socio-ecological systems.  _Regional Environmental Change_ 15, 1 (01 Jan 2015), 1–12.  https://doi.org/10.1007/s10113-014-0631-y
  * Davies et al. (2011) Zoe G Davies, Jill L Edmondson, Andreas Heinemeyer, Jonathan R Leake, and Kevin J Gaston. 2011\.  Mapping an urban ecosystem service: quantifying above-ground carbon storage at a city-wide scale.  _Journal of applied ecology_ 48, 5 (2011), 1125–1134.
  * Drusch et al. (2012) Matthias Drusch, Umberto Del Bello, Sébastien Carlier, Olivier Colin, Veronica Fernandez, Ferran Gascon, Bianca Hoersch, Claudia Isola, Paolo Laberinti, Philippe Martimort, et al. 2012\.  Sentinel-2: ESA’s optical high-resolution mission for GMES operational services.  _Remote sensing of Environment_ 120 (2012), 25–36.
  * Edenhofer (2015) Ottmar Edenhofer. 2015\.  _Climate change 2014: mitigation of climate change_. Vol. 3.  Cambridge University Press.
  * Eggleston et al. (2006) Simon Eggleston, Leandro Buendia, Kyoko Miwa, Todd Ngara, and Kiyoto Tanabe. 2006.  _2006 IPCC guidelines for national greenhouse gas inventories_. Vol. 5.  Institute for Global Environmental Strategies Hayama, Japan.
  * EPA (2019) United States Environmental Protection Agency EPA. 2019.  Greenhouse Gas Emissions from a Typical Passenger Vehicle.  https://www.epa.gov/greenvehicles/greenhouse-gas-emissions-typical-passenger-vehicle
  * Figueres et al. (2017) Christiana Figueres, Hans Joachim Schellnhuber, Gail Whiteman, Johan Rockström, Anthony Hobley, and Stefan Rahmstorf. 2017.  Three years to safeguard our climate.  _Nature News_ 546, 7660 (2017), 593.
  * GDAL/OGR contributors (2019) GDAL/OGR contributors. 2019\.  _GDAL/OGR Geospatial Data Abstraction software Library_.  Open Source Geospatial Foundation.  https://gdal.org
  * Gil et al. (2012a) Yolanda Gil, Ewa Deelman, Ibrahim Demir, Chris Duffy, Marlon Pierce, Suresh Marru, and Gerry Wiener. 2012a.  Designing a Roadmap for Workflow Cyberinfrastructure in the Geosciences: From Big Data to the Long Tail.  American Geophysical Union Fall Meeting, San Francisco, CA.
  * Gil et al. (2012b) Yolanda Gil, Varun Ratnakar, Ewa Deelman, and Christopher Mason. 2012b.  Using Semantic Workflows for Genome-Scale Analysis.  International Conference on Intelligent Systems for Molecular Biology (ISMB), Long Beach, CA.
  * Gil et al. (2011) Yolanda Gil, Varun Ratnakar, Jihie Kim, Pedro Antonio Gonzalez-Calero, Paul Groth, Joshua Moody, and Ewa Deelman. 2011\.  Wings: Intelligent Workflow-Based Design of Computational Experiments.  _IEEE Intelligent Systems_ 26, 1 (2011).  http://www.isi.edu/~gil/papers/gil-etal-ieee-is-11.pdf
  * Gillingham et al. (2018) Kenneth Gillingham, William Nordhaus, David Anthoff, Geoffrey Blanford, Valentina Bosetti, Peter Christensen, Haewon McJeon, and John Reilly. 2018.  Modeling uncertainty in integrated assessment of climate change: A multimodel comparison.  _Journal of the Association of Environmental and Resource Economists_ 5, 4 (2018), 791–826.
  * Gough (2016) Clair Gough. 2016\.  _Carbon capture and its storage: an integrated assessment_.  Routledge.
  * Grizonnet et al. (2017) Manuel Grizonnet, Julien Michel, Victor Poughon, Jordi Inglada, Mickaël Savinaud, and Rémi Cresson. 2017.  Orfeo ToolBox: Open source processing of remote sensing images.  _Open Geospatial Data, Software and Standards_ 2, 1 (2017), 15\.
  * Haklay and Weber (2008) Mordechai Haklay and Patrick Weber. 2008.  Openstreetmap: User-generated street maps.  _IEEE Pervasive Computing_ 7, 4 (2008), 12–18.
  * Incropera (2016) Frank P Incropera. 2016\.  _Climate change: a wicked problem: complexity and uncertainty at the intersection of science, economics, politics, and human behavior_.  Cambridge University Press.
  * Kagawa and Le Sourda (2018) Ayako Kagawa and Guillaume Le Sourda. 2018.  Mapping the world: cartographic and geographic visualization by the United Nations Geospatial Information Section (formerly Cartographic Section). In _Proceedings of the ICA_ , Vol. 1. 58.
  * Kundzewicz et al. (2018) ZW Kundzewicz, V Krysanova, RE Benestad, Ø Hov, M Piniewski, and IM Otto. 2018\.  Uncertainty in climate change impacts on water resources.  _Environmental Science & Policy_ 79 (2018), 1–8.
  * Livesley et al. (2016) SJ Livesley, EG McPherson, and C Calfapietra. 2016\.  The urban forest and ecosystem services: impacts on urban water, heat, and pollution cycles at the tree, street, and city scale.  _Journal of environmental quality_ 45, 1 (2016), 119–124.
  * Macharia et al. (2017) Peter M Macharia, Paul O Ouma, Ezekiel G Gogo, Robert W Snow, and Abdisalan M Noor. 2017.  Spatial accessibility to basic public health services in South Sudan.  _Geospatial health_ 12, 1 (2017), 510.
  * McGovern and Pasher (2016) Mark McGovern and Jon Pasher. 2016.  Canadian urban tree canopy cover and carbon sequestration status and change 1990–2012.  _Urban forestry & urban greening_ 20 (2016), 227–232.
  * Meadow et al. (2015) Alison M Meadow, Daniel B Ferguson, Zack Guido, Alexandra Horangic, Gigi Owen, and Tamara Wall. 2015\.  Moving toward the deliberate coproduction of climate science knowledge.  _Weather, Climate, and Society_ 7, 2 (2015), 179–191.
  * Nassif et al. (2016) Ayman Nassif, Ann Stewart, Millie Mutepfe, and Petros Christou. 2016.  Infrastructure needs in sub Saharan Africa with particular reference to South Sudan. In _Third Australasia and South-East Asia Structural Engineering and Construction Conference: ASEA-SEC-3_. ISEC Press.
  * Nowak et al. (2013) David J Nowak, Eric J Greenfield, Robert E Hoehn, and Elizabeth Lapoint. 2013.  Carbon storage and sequestration by trees in urban and community areas of the United States.  _Environmental pollution_ 178 (2013), 229–236.
  * Pasher et al. (2014) Jon Pasher, Mark McGovern, Michael Khoury, and Jason Duffe. 2014.  Assessing carbon storage and sequestration by Canada’s urban forests using high resolution earth observation data.  _Urban forestry & urban greening_ 13, 3 (2014), 484–494.
  * Penman et al. (2003) Jim Penman, Michael Gytarsky, Taka Hiraishi, Thelma Krug, Dina Kruger, Riitta Pipatti, Leandro Buendia, Kyoko Miwa, Todd Ngara, Kiyoto Tanabe, et al. 2003\.  Good practice guidance for land use, land-use change and forestry.  _Good practice guidance for land use, land-use change and forestry._ (2003).
  * Raciti et al. (2014) Steve M Raciti, Lucy R Hutyra, and Jared D Newell. 2014\.  Mapping carbon storage in urban trees with multi-source remote sensing data: Relationships between biomass, land use, and demographics in Boston neighborhoods.  _Science of the Total Environment_ 500 (2014), 72–83.
  * Schreyer et al. (2014) Johannes Schreyer, Jan Tigges, Tobia Lakes, and Galina Churkina. 2014.  Using airborne LiDAR and QuickBird data for modelling urban tree carbon storage and its distribution—A case study of Berlin.  _Remote Sensing_ 6, 11 (2014), 10636–10655.
  * Schulson (2018) Michael Schulson. 2018\.  Science’s Reproducibility Crisis Is Being Used as Political Ammunition.  _Wired_ (2018).  https://www.wired.com/story/sciences-reproducibility-crisis-is-being-used-as-political-ammunition/
  * Sørbø and Ahmed (2013) Gunnar M Sørbø and Abdel Ghaffar M Ahmed. 2013.  _Sudan divided: continuing conflict in a contested state_.  Springer.
  * Taylor et al. (2014) Ian J. Taylor, Ewa Deelman, Dennis B. Gannon, and Matthew Shields. 2014.  _Workflows for e-Science: Scientific Workflows for Grids_.  Springer Publishing Company, Incorporated.
  * Tigges et al. (2017) Jan Tigges, Galina Churkina, and Tobia Lakes. 2017\.  Modeling above-ground carbon storage: a remote sensing approach to derive individual tree species information in urban settings.  _Urban ecosystems_ 20, 1 (2017), 97–111.
  * Wolstencroft et al. (2013) Katherine Wolstencroft, Robert Haines, Donal Fellows, Alan Williams, David Withers, Stuart Owen, Stian Soiland-Reyes, Ian Dunlop, Aleksandra Nenadic, Paul Fisher, et al. 2013\.  The Taverna workflow suite: designing and executing workflows of Web Services on the desktop, web or in the cloud.  _Nucleic acids research_ 41, W1 (2013), W557–W561.


# Leading order track functions in a hot and dense QGP

João Barata<EMAIL_ADDRESS>Physics Department, Brookhaven National
Laboratory, Upton, NY 11973, USA  Robert Szafron<EMAIL_ADDRESS>Physics
Department, Brookhaven National Laboratory, Upton, NY 11973, USA

###### Abstract

We study the modifications to the fragmentation pattern of partons into
charged particles in the presence of a hot and dense Quark Gluon Plasma. To
this end, we analyze the perturbative renormalization group equations of the
track functions, which describe the energy fraction carried by charged
hadrons. Focusing on pure Yang-Mills theory, we compute the lowest order
moments of the medium-modified track functions, which are found to be
sensitive to the reduced phase space for emissions in the medium and to energy
loss. We use the extracted moments to calculate the Energy Energy Correlator
(EEC) on tracks in the collinear limit. The EEC on medium-evolved tracks does
not differ qualitatively from the EEC on vacuum tracks despite being sensitive
to the color decoherence transition and suppressing the distribution due to
quenching, as seen in other jet observables.

Introduction: Jets are an ideal tool to study the time evolution and spatial
structure of the Quark Gluon Plasma (QGP) produced in heavy ion collisions.
When using jets in this manner, two problems need to be addressed: how the QCD
matter modifies the jets’ structure Mehtar-Tani  _et al._ (2013); Qin and Wang
(2015) and which observables best reveal the medium imprints in the measured
distributions Cunqueiro and Sickles (2022); Apolinário  _et al._ (2022).

Regarding the latter point, we address an open question about Energy Energy
Correlator (EEC) measurements inside jets, i.e. in the collinear limit,
traveling through the QGP Andres  _et al._ (2023a, b); Barata  _et al._
(2023a); Barata and Mehtar-Tani (2023); Barata  _et al._ (2023b); Yang  _et
al._ (2023). The EEC is the simplest correlation function constructed from
quantum field theory light-ray operators measuring the asymptotic energy flow
passing through idealized calorimeter cells, see e.g. Chen (2023); Chen  _et
al._ (2020); Hofman and Maldacena (2008). It can be written in terms of the
inclusive two-particle cross-section as Basham  _et al._ (1978)

$\displaystyle\frac{d\Sigma^{(n)}}{d\chi}$ $\displaystyle=\sum_{\\{i,j\\}}\int
d\sigma\,\frac{E_{i}^{n}E_{j}^{n}}{p_{t}^{2n}}\delta(\vec{n}_{i}\cdot\vec{n}_{j}-\cos\chi)\,,$
(1)

where $\vec{n}_{i}$ denotes the spatial direction of the measured energy flow,
$E_{i}$ is the energy of particle $i$, $p_{t}$ is the total energy of the
system, and the sum is taken over all particles in the final state. The
integer $n\geq 1$ is introduced as an effective way to further suppress
contributions to the EEC from soft emissions.

Values of $n>1$ are especially relevant in the heavy ion context where
uncorrelated soft sources contaminate jets. In such cases, the EEC is
infrared-and-collinear (IRC) unsafe, and Eq. (1) cannot be evaluated
perturbatively. Nonetheless, one can separate the perturbative and non-
perturbative contributions to the EEC Chen  _et al._ (2020); Li  _et al._
(2022). The matching coefficients between the parton level EEC and the full
distribution are given by the moments of the so-called track functions (TFs)
Chang  _et al._ (2013a).111This issue can also be avoided by defining the EEC
on objects which are naturally IRC safe Barata  _et al._ (2023b).

Figure 1: Diagrammatic representation of the LO evolution kernel (see Eq.
(Leading order track functions in a hot and dense QGP)) and the different
medium modifications we take into account: energy loss (top left), induced
radiation (bottom left), and reduced phase space (right).

Several theoretical and phenomenological aspects of the track functions have
been explored for pure vacuum evolution, see e.g. Chang  _et al._ (2013a);
Chen  _et al._ (2020); Li  _et al._ (2022); Chang  _et al._ (2013b); Lee and
Moult (2023). They have not been studied in the heavy ion context, and the in-
medium EEC calculations have so far neglected the matching to the hadronic
level observable Andres  _et al._ (2023a); Barata  _et al._ (2023b). The goal
of this work is to provide the first phenomenologically driven study of the
medium effects on the track functions and explore how these affect the jet EEC
on tracks.

Leading order vacuum track functions: Tracks functions were first introduced a
decade ago Chang  _et al._ (2013a, b), allowing to compute jet observables
based on hadrons with specific quantum numbers in a theoretically well-defined
way; in particular TFs opened the possibility to compare theory calculations
to experimental measurements on charged tracks. In the simplest case, a track
function $T_{f}(x)$ describes the fragmentation of a parton of flavor $f$ into
a subset of hadrons222Here, we have in mind the fragmentation into charged
hadrons, but the calculation is valid for other conserved quantum numbers.
carrying an energy fraction $x$. Similarly to collinear fragmentation
functions (FFs) Collins and Soper (1982), TFs are intrinsically non-
perturbative. However, unlike FFs, they obey a highly non-linear
renormalization group (RG) evolution Jaarsma  _et al._ (2022); Li  _et al._
(2022); Chen  _et al._ (2022), more closely related to multi-hadron
fragmentation Sukhatme and Lassila (1980).333A detailed description of the RG
group for track functions has recently become available beyond leading order
(LO) Li  _et al._ (2022). At leading order (LO) in the strong coupling
$\alpha_{s}$, $T_{i}(x)$ satisfies the RG evolution equation Chang  _et al._
(2013a); Jaarsma  _et al._ (2022)

$\displaystyle\mu\frac{dT_{i}(x,\mu)}{d\mu}=\sum_{j,k}\int_{0}^{1}dz\,\frac{\alpha_{s}}{2\pi}P_{i\to
jk}(z)$
$\displaystyle\times\int_{x_{1},x_{2}}T_{j}(x_{1},\mu)T_{k}(x_{2},\mu)\delta(x-zx_{1}-(1-z)x_{2})\,,$
(2)

where $P_{i\to jk}$ are the LO regularized QCD splitting functions Altarelli
and Parisi (1977), and the evolution is expressed in terms of the
renormalization scale $\mu$, see Fig. 1. The DGLAP evolution Dokshitzer
(1977); Altarelli and Parisi (1977); Gribov and Lipatov (1972) can be
recovered by restricting Eq. (Leading order track functions in a hot and dense
QGP) to a one-body evolution, i.e., when it is fully inclusive by replacing
$T_{k}(x)\to\delta(x)$ Chen  _et al._ (2023), for a subset of particles.444 We
use the notation $\int_{x}$ to denote the integration over the full phase
space in the energy fraction $x$.

The RG equation simplifies for the track functions’ moments555The dependence
on the evolution scale $\mu$ is made implicit.

$\displaystyle T_{i}^{[N]}\equiv\int_{0}^{1}dz\,z^{N}T_{i}(z)\,.$ (3)

In Yang-Mills (YM) theory, it follows from Eq. (Leading order track functions
in a hot and dense QGP) that the moments satisfy

$\displaystyle\mu\frac{dT^{[N]}}{d\mu}=\frac{\alpha_{s}}{2\pi}\int_{0}^{1}dzP(z)$
$\displaystyle\times\int_{x_{1},x_{2}}T(x_{1},\mu)T(x_{2},\mu)(zx_{1}+(1-z)x_{2})^{N}\,,$
(4)

where $P(z)\equiv P_{g\to
gg}(z)=2N_{c}(z/(1-z)_{+}+(1-z)/z+z(1-z))+11N_{c}/6\,\delta(1-z)$ and
$T(z)=T_{g}(z)$. Using conservation of energy to impose the sum rule
$T^{[0]}=1$, and the symmetry of the splitting function under the integral
after shifting the poles to $z=1$, i.e., $\int_{z}P(z)(1-z)=\int_{z}P(z)z$,
the leading moments’ evolution reads

$\displaystyle\mu\frac{dT^{[1]}}{d\mu}=-\frac{\alpha_{s}}{\pi}\gamma(2)T^{[1]}\,,$
$\displaystyle\mu\frac{dT^{[2]}}{d\mu}=-\frac{\alpha_{s}}{\pi}\left(\gamma(3)T^{[2]}-(\gamma(2)-\gamma(3))T^{[1]}T^{[1]}\right)\,.$
(5)

Here, we have introduced the anomalous dimensions

$\displaystyle\gamma(j)$ $\displaystyle\equiv-\int_{0}^{1}dz\,z^{j-1}P(z)\,,$
(6)

such that in the vacuum $\gamma_{\rm
vac}(j)=-11N_{c}/6-2N_{c}(1/(j^{2}-j)+1/((j+1)(j+2))-\psi(j+1)-\gamma_{E}]$,
with $\psi(j)\equiv\Gamma^{\prime}(j)/\Gamma(j)$ the digamma function and
$\gamma_{E}\approx 0.577$ the Euler-Mascheroni constant. Equations (Leading
order track functions in a hot and dense QGP) can be further simplified in YM
theory, where $\gamma(2)=0$. This sum rule follows from energy conservation;
as a result, it provides a non-perturbative constraint that must also be
satisfied for in-medium evolution, as long as energy loss effects are not
present. The related RG invariance under the shift symmetry $T(x)\to T(x+b)$
Li  _et al._ (2022); Jaarsma  _et al._ (2022) constrains higher moments’
evolution, but there are no further non-trivial constraints in YM theory.

Medium modifications to the track functions: We assume that the partonic
cascade occurs in a homogeneous and isotropic gluonic QCD medium, with a
sizeable longitudinal extension $L$.666For recent developments in the
description of fragmentation in more realistic matter models see e.g. Sadofyev
_et al._ (2021); Kuzmin  _et al._ (2023); Barata  _et al._ (2023c).  We take
the interactions with the underlying matter to be dominated by multiple soft
gluon exchanges, see Fig. 1, which are captured in the BDMPS-Z/ASW
phenomenological model Zakharov (1996); Wiedemann and Gyulassy (1999); Baier
_et al._ (1995); Salgado and Wiedemann (2003). Large momentum exchanges with
the medium are neglected; the discussion can be easily extended to such cases.

We first consider the hardest/earliest branchings in the cascade, which are
primarily vacuum-like Caucal  _et al._ (2018). The phase space for these
emissions, at leading logarithmic accuracy, can be captured by introducing the
semi-classical formation time $t_{f}=2/(z(1-z)\theta^{2}p_{t})$; the radiation
pattern is modified depending on the ordering of $t_{f}$ with respect to the
characteristic time scale generated by the medium, $t_{c}$. In particular, the
vacuum-like region corresponds to $t_{f}<t_{c}$ Caucal  _et al._ (2018).
Alternatively, we can use angles to define the vacuum-like sector as
$\theta>\theta_{c}$, where $\theta_{c}$ is the characteristic medium angular
scale. For the matter model considered, the characteristic angular scale is
$\theta_{c}\equiv 2/\sqrt{\hat{q}_{\rm eff}L^{3}}$ Casalderrey-Solana  _et
al._ (2013); Mehtar-Tani  _et al._ (2012), with $\hat{q}_{\rm
eff}=\hat{q}(1-z+z^{2})\equiv\hat{q}f(z)$ the effective jet quenching
transport coefficient in YM theory, and $\hat{q}$ is the scalar adjoint
transport coefficient; the corresponding characteristic time scale is
$t_{c}=\sqrt{2z(1-z)p_{t}/\hat{q}_{\rm eff}}$ Mehtar-Tani and Tywoniuk
(2018a); Caucal  _et al._ (2018).

We consider two different phase spaces for vacuum emission, differing mainly
in the inclusion of out-of-the-medium fragmentation. Since the cascade still
obeys vacuum evolution, the functional form of Eq. (Leading order track
functions in a hot and dense QGP) is not affected, and the phase space
modification can be formally absorbed into the anomalous dimension Barata  _et
al._ (2023b)

$\displaystyle\gamma(j,\mu)=-\int_{0}^{1}dz\,P(z)\Theta_{\rm
PS}(\mu,z)z^{j-1}\,,$ (7)

where $\Theta_{\rm PS}$ denotes the allowed phase space for fragmentation. The
first model only allows for branchings to occur if they are sufficiently
short-lived and wide-angled to be resolved Mehtar-Tani and Tywoniuk (2018a):
$\Theta^{(1)}_{\rm PS}=\Theta(\mu-p_{t}\theta_{c})\Theta(t_{c}-t_{f})$, where
the scale $\mu=p_{t}\theta$. The second phase space we study allows
fragmentation in the region $t_{f}>L$ Caucal  _et al._ (2018):
$\Theta^{(2)}_{\rm PS}=1-\Theta(t_{f}-t_{c})\Theta(L-t_{f})$. Both phase
spaces are formally only well described in the $z\ll 1$ limit, but we naively
extend them to full kinematics here. A study of the phase space structure at
finite energy fractions is still missing from the literature.

To better understand the medium modifications, we first compute the anomalous
dimensions in the soft limit. Using the soft (BFKL) form of the vacuum
anomalous dimension, $\gamma^{\rm BFKL}_{\rm vac}(j)=\frac{2N_{c}}{1-j}$, we
find that:

$\displaystyle\gamma^{(i)}(j,\mu)$ $\displaystyle=\gamma^{\rm BFKL}_{\rm
vac}(j)\left(\lambda_{i+}^{j-1}-\lambda_{i-}^{j-1}+\delta_{i,2}\right)\,,$ (8)

where $i=1,2$ denotes the phase space model,
$\lambda_{1-}=\lambda_{2+}=\sqrt[3]{2\hat{q}p_{t}/\mu^{4}}$, $\lambda_{1+}=1$
, $\lambda_{2-}=2p_{t}/(\mu^{2}L)$, and the phase space constraints depending
only on $\mu$ are left implicit. Inside the medium, the BFKL pole vanishes
since soft gluons cannot go on-shell. Also, Eq. (8) ensures the strongest
medium modifications occur at smaller integer values $j>1$. In Fig. 2, we show
the anomalous dimensions for both models at finite $z$, while imposing the sum
rule $\gamma(2)=0$. On top, we show the evolution with the scale $\mu$. As
expected, if $\mu\approx p_{t}$ (i.e. $\theta\approx 1$), we recover the
vacuum result as there is no loss of phase space, i.e., highly energetic jets
are not sensitive to the medium. Taking the evolution scale down, the medium
reduces the phase space available to radiate, and thus, the anomalous
dimension also decreases. The significant difference between the two models
resides in the $\mu\to 0$ region (at fixed $p_{t}$): while $\gamma^{(1)}\to
0$, $\gamma^{(2)}\to\gamma_{\rm vac}$ since at small angles out-of-the-medium
emissions are included, see Fig. 3 in Barata  _et al._ (2023b). At fixed
$\mu$, the denser the medium, the larger the modification, as expected. These
are more prominent for $\gamma^{(1)}$, where no out-of-the-medium evolution is
incorporated, and for larger $j$. The latter originate from intermediate
values of $z$, explaining the difference to the soft limit.

Figure 2: Evolution of the anomalous dimensions as a function of the evolution
scale $\mu$ and the jet quenching parameter $\hat{q}$.

As illustrated in Fig. 1, the evolution in the medium also results in the
production of stimulated bremsstrahlung radiation, in parallel to the
virtuality cascade. This mechanism has two significant phenomenological
consequences: the enhancement of the splitting function at large angles and
radiative energy loss to the matter. We now describe how these effects can be
incorporated into the evolution of the track functions.

The multiple soft momentum gluons exchanged with the medium result in a
characteristic interference pattern, leading to an enhanced branching rate,
first computed by Landau, Pomeranchuk, and Migdal in QED Landau and
Pomeranchuk (1953); Migdal (1956); Klein  _et al._ (1994). This regime,
captured in the BDMPS-ASW approximation, is valid for intermediate energies of
the radiated partons, $\omega_{\rm BH}<z(1-z)p_{t}<\omega_{c}$, where
$\omega_{\rm BH}$ denotes the Bethe-Heitler frequency, below which the
radiative spectrum is dominated by incoherent contributions Bethe and Heitler
(1934), and $\omega_{c}=\frac{1}{2}\hat{q}L^{2}$ denotes the characteristic
frequency above which the gluon resolves the medium coherently. A closed-form
for the radiation rate is not known in QCD in this region of phase space, and
we consider an asymptotic limit of exactly collinear branchings followed by
in-medium diffusion:777See e.g. Barata  _et al._ (2021a) for further
discussion on this approximation.

$\displaystyle\frac{dI}{dzd\mu\,z^{2}(1-z)^{2}\mu}=\int_{0}^{L}dt\,\frac{{\cal
P}(\mu z(1-z),L-t)}{2\pi}\frac{dI}{dzdt}\,,$ (9)

where $dI/dzdt$ denotes the medium induced radiation rate and ${\cal P}$ is
the so-called momentum broadening distribution. The transverse momentum is
only acquired due to the evolution between the time the outgoing states go on-
shell and the end of the medium; the in-medium rate determines the energy
distribution. In the BDMPS-Z/ASW approximation the broadening distribution is
purely diffusive Blaizot  _et al._ (2014); Barata  _et al._ (2021b), ${\cal
P}(k,t)=4\pi/(\hat{q}_{\rm eff}t)\exp(-k^{2}/(\hat{q}_{\rm eff}t))$; the
radiative rate reads Blaizot  _et al._ (2013a, 2014)

$\displaystyle\frac{dI}{dzdt}$
$\displaystyle=\frac{\alpha_{s}N_{c}}{2\pi}\frac{f^{\frac{5}{2}}(z)}{(z(1-z))^{\frac{3}{2}}}\sqrt{\frac{\hat{q}}{p_{t}}}\,.$
(10)

Eq. (10) can be regularized by a plus distribution while shifting the
singularities to $z=1$, and using the integral identity Blaizot  _et al._
(2014)

$\displaystyle\int_{\varepsilon}^{1-\varepsilon}dz\,\frac{dI}{dzdt}=\int_{0}^{1-\varepsilon}dz\,2z\,\frac{dI}{dzdt}\,.$
(11)

Combining all these elements, Eq. (9) becomes

$\displaystyle\frac{dI}{dzd\mu}=\frac{\alpha_{s}N_{c}}{\pi}\sqrt{\frac{1}{\hat{q}\,p_{t}}}\mu$
$\displaystyle\times\frac{(2z)\,z^{2}(1-z)^{2}f^{\frac{3}{2}}(z)}{(z(1-z))^{\frac{3}{2}}}\Gamma_{0}\left(\frac{z^{2}(1-z)^{2}\mu^{2}}{f(z)\hat{q}L}\right)\,,$
(12)

where $\Gamma_{0}$ denotes the incomplete Gamma function.

We include the medium-induced kernel in the splitting functions entering the
RG evolution in Eq. (Leading order track functions in a hot and dense QGP).
Such a choice has been applied in different jet quenching phenomenological
studies, see e.g Ke and Vitev (2023); Deng and Wang (2010); Majumder (2013);
Armesto  _et al._ (2009); Chien  _et al._ (2016). Note that the in-medium
effects impact both the anomalous dimensions and matrix elements of the
operators, i.e., the boundary conditions for the evolution. We start the
evolution from the hard scale, where the in-medium effects can be treated
perturbatively. Thus, in an effective field theory spirit, the hard scale is
assumed to be parametrically much larger than medium-induced scales, and the
dominant effect comes from the evolution, i.e. from in-medium modification of
the anomalous dimensions. A complete effective field theory formulation of in-
medium effects is still missing; although such an effort would be a worthwhile
endeavor, it goes past the scope of this study, and we leave it for future
work. Thus, the anomalous dimensions in Eq. (6) become

$\displaystyle\gamma(j,\mu)=-\int_{0}^{1}\Bigg{\\{}P(z)+N_{c}\,\mu^{2}\sqrt{\frac{1}{\hat{q}\,p_{t}}}$
$\displaystyle\times
2z\sqrt{z(1-z)}\,f^{\frac{3}{2}}(z)\Gamma_{0}\left(\frac{z^{2}(1-z)^{2}\mu^{2}}{f(z)\hat{q}L}\right)\Bigg{\\}}z^{j-1}\,,$
(13)

up to a $\delta(z-1)$ term, whose coefficient can be obtained by imposing the
sum rule $\gamma(2)=0$. In Fig. 3 we evaluate Eq. (Leading order track
functions in a hot and dense QGP). To ensure that the collinear cascade
approximation is valid, we require $p_{t}\ll\hat{q}L^{2}$ Blaizot  _et al._
(2013a). Therefore, we take a long medium with $L=10$ fm and a small jet
energy with $p_{t}=15$ GeV. The results show a mild dependence on $j$, while
the evolution in $\mu$ is fast, leading to a substantial enhancement that
competes with the phase space constraint effect in Fig. 2. The magnitude of
the increase is related to the asymptotic approximation used; although other
calculations of the in-medium rate at finite $z$ exist Isaksen and Tywoniuk
(2020, 2023); Domínguez  _et al._ (2020); Apolinário  _et al._ (2015); Blaizot
_et al._ (2013b), it is unclear how to incorporate them in the RG evolution at
present.888This is possible for thin or dilute matter, see Sievert and Vitev
(2018). Due to the high sensitivity to the modeling, we will not discuss these
contributions further; nonetheless, we expect this modification to yield a
small effect due to the absence of a collinear pole in the medium.

Figure 3: Anomalous dimension evolution including medium induced radiation.

So far, we have considered the contribution of bremsstrahlung radiation, which
hadronizes into charged hadrons. However, a substantial part of the softer
induced radiation feeds down to the medium scale, and it should be viewed as
part of the QGP while quenching the jet energy. We describe these lost
emissions in terms of a distribution function $D(\varepsilon)$, with
$\varepsilon\ll p_{t}$ representing the amount of quenched energy, see Fig. 1.
This distribution can be computed perturbatively under certain limits Baier
_et al._ (1995); Mehtar-Tani and Tywoniuk (2018b, 2017), but it is expected to
have a significant non-perturbative contribution. Therefore, we remain
agnostic to its functional form and only impose overall energy conservation,
i.e., $\int_{0}^{\infty}d\varepsilon D(\varepsilon)=1$.

To implement the effect of energy loss in the TFs’ RG, we first consider the
simpler case of DGLAP evolution; at $\mathcal{O}(\alpha_{s}^{2})$ the real
contribution to the gluon distribution $G$ reads Altarelli and Parisi (1977)

$\displaystyle G(x,\mu)\Big{|}_{\rm
real}^{\mathcal{O}(\alpha_{s}^{2})}=\frac{\alpha_{s}^{2}}{(2\pi)^{2}}\int_{z_{1},z_{2},x_{1}}\,\hat{P}(z_{1})\hat{P}(z_{2})$
$\displaystyle\times\int_{\mu_{0}}^{\mu}d\ln\mu^{\prime}\int_{\mu_{0}}^{\mu^{\prime}}d\ln\mu^{\prime\prime}G(x_{1},\mu^{\prime\prime})\delta(x-x_{1}z_{1}(1-z_{2}))\,,$
(14)

where $\hat{P}(z)=\hat{P}(1-z)$ is the unregularized splitting function. We
choose $z_{2}$ to represent the fraction of energy lost to the matter. The RHS
of Eq. (Leading order track functions in a hot and dense QGP) can be
manipulated the following way (dropping the $\mu$ dependence):

$\displaystyle\frac{\alpha_{s}^{2}}{(2\pi)^{2}}\int_{x}^{1}dz_{1}\int_{0}^{1-\frac{x}{z_{1}}}dz_{2}\,\frac{\hat{P}(z_{1})\hat{P}(z_{2})}{z_{1}(1-z_{2})}G\left(\frac{x}{z_{1}(1-z_{2})}\right)$
$\displaystyle\approx\frac{\alpha_{s}}{2\pi}\int_{x}^{1}\frac{dz_{1}}{z_{1}}\int_{0}^{(z_{1}-x)p_{t}}d\varepsilon\,\hat{P}\left(z_{1}+\frac{\varepsilon}{p_{t}}\right)D\left(\varepsilon\right)G\left(\frac{x}{z_{1}}\right)\,,$
(15)

where we replaced the second splitting function by the energy loss
distribution $D(\varepsilon)$, modeling energy transferred to the matter, see
Fig. 1. Assuming $1>x\leavevmode\nobreak\ \sim
z_{1}\gg\frac{\varepsilon}{p_{t}}$, the upper cut-off of the $\varepsilon$
integration should be extended to infinity. Treating the energy loss as an
additional convolution within the standard DGLAP cascade and including the
virtual corrections, we find

$\displaystyle\mu\frac{dG(x,\mu)}{d\mu}$
$\displaystyle\approx\frac{\alpha_{s}}{2\pi}\int_{x_{1},z}\int_{0}^{\infty}d\varepsilon\,D(\varepsilon)P\left(z+\frac{\varepsilon}{p_{t}}\right)$
$\displaystyle\times G(x_{1})\delta(x-x_{1}z)\,,$ (16)

which agrees with the Sudakov factor found in Mehtar-Tani and Tywoniuk (2017).
A similar procedure for the TFs, with the initial truncation in Eq. (Leading
order track functions in a hot and dense QGP) at
$\mathcal{O}(\alpha_{s}^{3})$, with
$\delta(x-x_{1}z(1-z_{2}))\to\delta(x-x_{1}z_{1}(1-z_{2})-x_{2}(1-z_{1})(1-z_{3}))$,
where $z_{2},z_{3}$ are the energy loss fractions from uncorrelated emissions,
leads to

$\displaystyle\mu\frac{dT}{d\mu}=\frac{\alpha_{s}}{2\pi}\int_{\varepsilon_{1},\varepsilon_{2}}D(\varepsilon_{1})D(\varepsilon_{2})\int_{0}^{1}dzP\left(z+\frac{\varepsilon_{1}}{p_{t}}+\frac{\varepsilon_{2}}{p_{t}}\right)$
$\displaystyle\int_{x_{1},x_{2}}T(x_{1},\mu)T(x_{2},\mu)\,\delta(x-zx_{1}-(1-z)x_{2})\,.$
(17)

The upper limit of the $\varepsilon_{1,2}$ integrals is set to infinity. The
evolution for the moments at LO is unchanged with respect to Eq. (Leading
order track functions in a hot and dense QGP) if the anomalous dimensions
include the double convolution with the energy loss probabilities.

Without specifying the form for $D(\varepsilon)$, but assuming that the energy
loss distribution has a small dispersion, we simplify the evolution equation
replacing $\varepsilon$ by its mean value
$\langle\varepsilon\rangle\sim\hat{q}L^{2}$. We note that this approximation
might be insufficient for phenomenological applications, see e.g. Baier  _et
al._ (2001). With this caveat, the anomalous dimension reduces to

$\displaystyle\gamma(j)$
$\displaystyle=-\int_{\frac{2\langle\varepsilon\rangle}{p_{t}}}^{1}dz\,P(z)\left(z-\frac{2\langle\varepsilon\rangle}{p_{t}}\right)^{j-1}\,,$
(18)

up to a $\delta(1-z)$ term. We can no longer fix its coefficient via the sum
rule $\gamma(2)=0$ since we assumed $1>x\gg\frac{\varepsilon}{p_{t}}$.
Nonetheless, we impose the sum rule in Fig. 4 to directly compare energy-loss
driven modifications to Figs. 2 and 3. Compared to the naive estimation from
the BDMPS-Z/ASW approximation, $\langle\varepsilon\rangle_{\rm
BDMPS-Z/ASW}\sim\mathcal{O}(20)$ GeV, we observe significant effects already
for $\langle\varepsilon\rangle$ one order of magnitude lower (red markers).
This is likely due to the mean value approximation, which typically
overestimates the quenching effects Baier  _et al._ (2001). Consequently, we
restrict ourselves to small values for $\langle\varepsilon\rangle/p_{t}$, such
that non-physical results are avoided.

Figure 4: Anomalous dimension evolution including average energy loss effects.
We use the same parameters as in Fig. 2.

Application to EECs: Having considered the modifications to the TFs’ RG
evolution, we study their impact on the jet EEC. The LO EEC on tracks is
defined as Chen  _et al._ (2020); Li  _et al._ (2022)

$\displaystyle\frac{d\Sigma^{(n)}_{\rm trk}}{d\chi}$
$\displaystyle\equiv\int_{z,x_{1},x_{2}}x_{1}^{n}T(x_{1})x_{2}^{n}T(x_{2})z^{n}(1-z)^{n}\frac{d\sigma}{\sigma
dzd\chi}$ $\displaystyle=\int_{z}T^{[n]}(\chi p_{t})T^{[n]}(\chi
p_{t})z^{n}(1-z)^{n}\frac{d\sigma}{\sigma dzd\chi}\,,$ (19)

where we only include terms for $\chi>0$. In YM theory, the lowest order non-
trivial track function appears at $n=2$; see Eq. (Leading order track
functions in a hot and dense QGP), which we consider in what follows. We work
at fixed coupling, such that we are only sensitive to the evolution of the
track functions, and for the initial condition, we take the simplistic form
$T(x,\mu_{0})=252x^{2}(1-x)^{6}$ at $\mu_{0}=p_{t}$ Chang  _et al._ (2013b).
Equation (Leading order track functions in a hot and dense QGP) is computed
using the in-medium cross-section obtained in the limit of highly energetic
final states, see Domínguez  _et al._ (2020); Isaksen and Tywoniuk (2020,
2023) for details. Writing $d\sigma=d\sigma_{\rm vac}\left(1+F_{\rm
med}\right)$, we have

$\displaystyle F_{\rm med}$
$\displaystyle=\frac{2}{t_{f}}\Bigg{\\{}\int_{0}^{L}dt\,\int_{t}^{L}\frac{dt^{\prime}}{t_{f}}\,\left[\cos\left(\frac{t^{\prime}-t}{t_{f}}\right)\mathcal{C}_{3}(t^{\prime},t)\mathcal{C}_{4}(L,t^{\prime})\right]$
$\displaystyle-\sin\left(\frac{L-t}{t_{f}}\right)\mathcal{C}_{3}(L,t)\Bigg{\\}}\,.$
(20)

The correlators $\mathcal{C}_{3,4}$ in the BDMPS-Z/ASW approximation take the
form $\mathcal{C}_{3}(t^{\prime},t)\approx e^{-\frac{1}{12}\hat{q}\chi^{2}\
(t^{\prime}-t)^{3}}$ and $\mathcal{C}_{4}(L,t^{\prime})\approx
e^{-\frac{1}{4}\hat{q}\chi^{2}(L-t^{\prime})(t^{\prime}-t)^{2}}$, where we
neglected the $z$ dependence in the phases.

Figure 5: Result for the EEC with $n=2$, including phase space effects. Top:
$p_{t}=50$ GeV; Bottom: $p_{t}=200$ GeV.

The results for $\Sigma^{(2)}_{\rm trk}$ are shown in Figs. 5 and 6. For the
case of reduced phase space, including the medium-modified track functions
leads to small changes compared to vacuum TFs. The major differences arise at
smaller values of $p_{t}$, in the region where the medium leads to an
enhancement of the EEC. Importantly, including the in-medium tracks does not
alter the shape of the distribution, which allows us to qualitatively
visualize the color decoherence transition. The distinction between the phase
space models is clear for small $\chi$, where out-of-the-medium fragmentation
becomes crucial. Mean energy loss shifts the distribution below the vacuum
case, as is typical of energy loss effects in jet observables. Again, the
overall shape of the distribution is only barely modified, even for the
largest $\langle\varepsilon\rangle$.

Figure 6: Result for the EEC with $n=2$, including mean energy loss. Top:
$\langle\varepsilon\rangle=0.4$ GeV; Bottom: $\langle\varepsilon\rangle=1$
GeV.

Conclusions: We conducted the first study of the track functions RG evolution
in a dense QGP environment. Comprehending such effects is vital in heavy ions,
where soft, uncorrelated sources contaminate the jet EEC. We showed that, at a
qualitative level, it is adequate to use the vacuum TFs as the QGP medium does
not drastically impact the shape of the EEC distribution through the modified
TFs.

Our analysis is phenomenologically driven and faces several shortcomings, such
as the need to describe the phase space for vacuum-like emissions at leading
logarithmic accuracy and the lack of suitable treatment of the in-medium
matrix elements. From a phenomenological point of view, it would be worthwhile
to extract the in-medium TFs from available data or MC samples. We hope to
address these aspects in the future.

Acknowledgments: JB and RS are supported by the United States Department of
Energy under Grant Contract DESC0012704. We are grateful to Swagato Mukherjee,
Guilherme Milhano, Yacine Mehtar-Tani, Ian Moult, Wenqing Fan, Andrey
Sadofyev, Xin-Nian Wang, and Mateusz Ploskon for discussions. We particularly
thank Paul Caucal, Pier Monni, and Alba Soto-Ontoso for their collaborative
efforts on closely related work and helpful discussions.

## References

  * Mehtar-Tani  _et al._ (2013) Y. Mehtar-Tani, J. G. Milhano, and K. Tywoniuk, Int. J. Mod. Phys. A 28, 1340013 (2013), arXiv:1302.2579 [hep-ph] .
  * Qin and Wang (2015) G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E 24, 1530014 (2015), arXiv:1511.00790 [hep-ph] .
  * Cunqueiro and Sickles (2022) L. Cunqueiro and A. M. Sickles, Prog. Part. Nucl. Phys. 124, 103940 (2022), arXiv:2110.14490 [nucl-ex] .
  * 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] .
  * Andres  _et al._ (2023a) C. Andres, F. Dominguez, J. Holguin, C. Marquet, and I. Moult, JHEP 09, 088 (2023a), arXiv:2303.03413 [hep-ph] .
  * Andres  _et al._ (2023b) C. Andres, F. Dominguez, R. Kunnawalkam Elayavalli, J. Holguin, C. Marquet, and I. Moult, Phys. Rev. Lett. 130, 262301 (2023b), arXiv:2209.11236 [hep-ph] .
  * Barata  _et al._ (2023a) J. Barata, J. G. Milhano, and A. V. Sadofyev, (2023a), arXiv:2308.01294 [hep-ph] .
  * Barata and Mehtar-Tani (2023) J. Barata and Y. Mehtar-Tani, in  _11th International Conference on Hard and Electromagnetic Probes of High-Energy Nuclear Collisions: Hard Probes 2023_ (2023) arXiv:2307.08943 [hep-ph] .
  * Barata  _et al._ (2023b) J. Barata, P. Caucal, A. Soto-Ontoso, and R. Szafron, (2023b), arXiv:2312.12527 [hep-ph] .
  * Yang  _et al._ (2023) Z. Yang, Y. He, I. Moult, and X.-N. Wang, (2023), arXiv:2310.01500 [hep-ph] .
  * Chen (2023) H. Chen, (2023), arXiv:2311.00350 [hep-ph] .
  * Chen  _et al._ (2020) H. Chen, I. Moult, X. Zhang, and H. X. Zhu, Phys. Rev. D 102, 054012 (2020), arXiv:2004.11381 [hep-ph] .
  * Hofman and Maldacena (2008) D. M. Hofman and J. Maldacena, JHEP 05, 012 (2008), arXiv:0803.1467 [hep-th] .
  * Basham  _et al._ (1978) C. L. Basham, L. S. Brown, S. D. Ellis, and S. T. Love, Phys. Rev. D 17, 2298 (1978).
  * Li  _et al._ (2022) Y. Li, I. Moult, S. S. van Velzen, W. J. Waalewijn, and H. X. Zhu, Phys. Rev. Lett. 128, 182001 (2022), arXiv:2108.01674 [hep-ph] .
  * Chang  _et al._ (2013a) H.-M. Chang, M. Procura, J. Thaler, and W. J. Waalewijn, Phys. Rev. Lett. 111, 102002 (2013a), arXiv:1303.6637 [hep-ph] .
  * Chang  _et al._ (2013b) H.-M. Chang, M. Procura, J. Thaler, and W. J. Waalewijn, Phys. Rev. D 88, 034030 (2013b), arXiv:1306.6630 [hep-ph] .
  * Lee and Moult (2023) K. Lee and I. Moult, (2023), arXiv:2308.00746 [hep-ph] .
  * Collins and Soper (1982) J. C. Collins and D. E. Soper, Nucl. Phys. B 194, 445 (1982).
  * Jaarsma  _et al._ (2022) M. Jaarsma, Y. Li, I. Moult, W. Waalewijn, and H. X. Zhu, JHEP 06, 139 (2022), arXiv:2201.05166 [hep-ph] .
  * Chen  _et al._ (2022) H. Chen, M. Jaarsma, Y. Li, I. Moult, W. J. Waalewijn, and H. X. Zhu, (2022), arXiv:2210.10061 [hep-ph] .
  * Sukhatme and Lassila (1980) U. P. Sukhatme and K. E. Lassila, Phys. Rev. D 22, 1184 (1980).
  * Altarelli and Parisi (1977) G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
  * Dokshitzer (1977) Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977).
  * Gribov and Lipatov (1972) V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972).
  * Chen  _et al._ (2023) H. Chen, M. Jaarsma, Y. Li, I. Moult, W. J. Waalewijn, and H. X. Zhu, JHEP 07, 185 (2023), arXiv:2210.10058 [hep-ph] .
  * Sadofyev  _et al._ (2021) A. V. Sadofyev, M. D. Sievert, and I. Vitev, Phys. Rev. D 104, 094044 (2021), arXiv:2104.09513 [hep-ph] .
  * Kuzmin  _et al._ (2023) M. V. Kuzmin, X. Mayo López, J. Reiten, and A. V. Sadofyev, (2023), arXiv:2309.00683 [hep-ph] .
  * Barata  _et al._ (2023c) J. Barata, X. Mayo López, A. V. Sadofyev, and C. A. Salgado, Phys. Rev. D 108, 034018 (2023c), arXiv:2304.03712 [hep-ph] .
  * Zakharov (1996) B. G. Zakharov, JETP Lett. 63, 952 (1996), arXiv:hep-ph/9607440 .
  * Wiedemann and Gyulassy (1999) U. A. Wiedemann and M. Gyulassy, Nucl. Phys. B 560, 345 (1999), arXiv:hep-ph/9906257 .
  * Baier  _et al._ (1995) R. Baier, Y. L. Dokshitzer, S. Peigne, and D. Schiff, Phys. Lett. B 345, 277 (1995), arXiv:hep-ph/9411409 .
  * Salgado and Wiedemann (2003) C. A. Salgado and U. A. Wiedemann, Phys. Rev. D 68, 014008 (2003), arXiv:hep-ph/0302184 .
  * Caucal  _et al._ (2018) P. Caucal, E. Iancu, A. H. Mueller, and G. Soyez, Phys. Rev. Lett. 120, 232001 (2018), arXiv:1801.09703 [hep-ph] .
  * Casalderrey-Solana  _et al._ (2013) J. Casalderrey-Solana, Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, Phys. Lett. B 725, 357 (2013), arXiv:1210.7765 [hep-ph] .
  * Mehtar-Tani  _et al._ (2012) Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, JHEP 10, 197 (2012), arXiv:1205.5739 [hep-ph] .
  * Mehtar-Tani and Tywoniuk (2018a) Y. Mehtar-Tani and K. Tywoniuk, Phys. Rev. D 98, 051501 (2018a), arXiv:1707.07361 [hep-ph] .
  * Landau and Pomeranchuk (1953) L. D. Landau and I. Pomeranchuk, Dokl. Akad. Nauk Ser. Fiz. 92, 535 (1953).
  * Migdal (1956) A. B. Migdal, Phys. Rev. 103, 1811 (1956).
  * Klein  _et al._ (1994) S. R. Klein _et al._ , AIP Conf. Proc. 302, 172 (1994).
  * Bethe and Heitler (1934) H. Bethe and W. Heitler, Proc. Roy. Soc. Lond. A 146, 83 (1934).
  * Barata  _et al._ (2021a) J. Barata, Y. Mehtar-Tani, A. Soto-Ontoso, and K. Tywoniuk, JHEP 09, 153 (2021a), arXiv:2106.07402 [hep-ph] .
  * Blaizot  _et al._ (2014) J.-P. Blaizot, F. Dominguez, E. Iancu, and Y. Mehtar-Tani, JHEP 06, 075 (2014), arXiv:1311.5823 [hep-ph] .
  * Barata  _et al._ (2021b) J. Barata, Y. Mehtar-Tani, A. Soto-Ontoso, and K. Tywoniuk, Phys. Rev. D 104, 054047 (2021b), arXiv:2009.13667 [hep-ph] .
  * Blaizot  _et al._ (2013a) J.-P. Blaizot, E. Iancu, and Y. Mehtar-Tani, Phys. Rev. Lett. 111, 052001 (2013a), arXiv:1301.6102 [hep-ph] .
  * Ke and Vitev (2023) W. Ke and I. Vitev, (2023), arXiv:2301.11940 [hep-ph] .
  * Deng and Wang (2010) W.-t. Deng and X.-N. Wang, Phys. Rev. C 81, 024902 (2010), arXiv:0910.3403 [hep-ph] .
  * Majumder (2013) A. Majumder, Phys. Rev. C 88, 014909 (2013), arXiv:1301.5323 [nucl-th] .
  * Armesto  _et al._ (2009) N. Armesto, L. Cunqueiro, and C. A. Salgado, Eur. Phys. J. C 63, 679 (2009), arXiv:0907.1014 [hep-ph] .
  * Chien  _et al._ (2016) Y.-T. Chien, A. Emerman, Z.-B. Kang, G. Ovanesyan, and I. Vitev, Phys. Rev. D 93, 074030 (2016), arXiv:1509.02936 [hep-ph] .
  * Isaksen and Tywoniuk (2020) J. H. Isaksen and K. Tywoniuk, JHEP 21, 125 (2020), arXiv:2107.02542 [hep-ph] .
  * Isaksen and Tywoniuk (2023) J. H. Isaksen and K. Tywoniuk, JHEP 09, 049 (2023), arXiv:2303.12119 [hep-ph] .
  * Domínguez  _et al._ (2020) F. Domínguez, J. G. Milhano, C. A. Salgado, K. Tywoniuk, and V. Vila, Eur. Phys. J. C 80, 11 (2020), arXiv:1907.03653 [hep-ph] .
  * Apolinário  _et al._ (2015) L. Apolinário, N. Armesto, J. G. Milhano, and C. A. Salgado, JHEP 02, 119 (2015), arXiv:1407.0599 [hep-ph] .
  * Blaizot  _et al._ (2013b) J.-P. Blaizot, F. Dominguez, E. Iancu, and Y. Mehtar-Tani, JHEP 01, 143 (2013b), arXiv:1209.4585 [hep-ph] .
  * Sievert and Vitev (2018) M. D. Sievert and I. Vitev, Phys. Rev. D 98, 094010 (2018), arXiv:1807.03799 [hep-ph] .
  * Mehtar-Tani and Tywoniuk (2018b) Y. Mehtar-Tani and K. Tywoniuk, Nucl. Phys. A 979, 165 (2018b), arXiv:1706.06047 [hep-ph] .
  * Mehtar-Tani and Tywoniuk (2017) Y. Mehtar-Tani and K. Tywoniuk, JHEP 04, 125 (2017), arXiv:1610.08930 [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 .