# Probing the interaction of semi-hard quarks and gluons with the underlying
event in light- and heavy-flavor triggered proton-proton collisions

Róbert Vértesi  Gyula Bencédi<EMAIL_ADDRESS>Anett Misák  Antonio
Ortiz  Wigner Research Centre for Physics, Hungary  Eötvös Loránd University
Faculty of Sciences  Instituto de Ciencias Nucleares, Universidad Nacional
Autónoma de México

###### Abstract

We study underlying-event observables in inelastic proton-proton (pp)
collisions at a centre-of-mass energy of $\sqrt{s}=13$ TeV with identified
light and heavy-flavor triggers using the PYTHIA 8 event generator. The study
is performed as a function of the transverse momentum of the leading particle
($p_{\mathrm{T}}^{\rm trigger}$). While at high $p_{\mathrm{T}}^{\rm trigger}$
($>10$ GeV/$c$) the underlying-event activity is independent of the leading
particle species, at intermediate $p_{\mathrm{T}}^{\rm trigger}$
($2<p_{\mathrm{T}}^{\rm trigger}<8$ GeV/$c$) it is larger in pion-triggered
events than in events triggered with B mesons. Moreover, the underlying event
in pion-triggered events, the majority of which are initiated by gluons, shows
a stronger effect of color reconnection than events triggered with B-hadrons,
that are mostly initiated by quark jets. The effect is observed at both
hadronic and partonic level. Given that color reconnection affects the
interaction among final partons before the hadronization, and that in the
string model quarks (gluons) are connected to one (two) string piece(s), we
conclude that the observed effect can be attributed to differences in the
interactions of gluon and quark jets with the underlying event.

###### keywords:

Heavy flavor, underlying event, proton-proton collision, PYTHIA 8 Monte Carlo

††journal: Journal of LaTeX Templates

##  1 Introduction

Long-range multi-particle correlations in the final state have been considered
as a signature of hydrodynamic evolution of the strongly coupled quark-gluon
plasma (sQGP) [1, 2, 3, 4]. LHC results, however, unexpectedly revealed that
high-multiplicity proton-proton events exhibit similar collective behavior,
with substantial azimuthal anisotropy ($v_{n}$) [5, 6]. Part of these effects
can be explained with vacuum-QCD processes at the soft-hard boundary, such as
multiple-parton interactions (MPI) [7], although the creation of sQGP in a
small volume cannot be completely ruled out [8, 9, 10, 11]. It is worth
mentioning that the experimental data support the presence of MPI in high-
energy hadronic interactions, see e.g. Refs. [12, 13]. In models that handle
color flow in softer MPI processes in a simplified manner, an additional step
is applied [14]. This step, called the color reconnection (CR), has a major
role in forming the multiplicity as well as the transverse momentum
($p_{\mathrm{T}}$) distribution of the final state particles, and may be a key
ingredient in the formation of collective behavior [15, 16]. Hence vacuum QCD
effects establish connection between the leading hard process and the
underlying event (UE) from particle production in soft and secondary hard
processes as well as from beam remnants [17]. The possibility of this
connection resulting in the observed collectivity was explored in
phenomenology calculations [18, 19, 20] and later corroborated in the
experiment [21, 22]. These works, however, mostly concentrated on light-flavor
and strangeness observables. Jet fragmentation is affected by the flavor of
the initial parton in several ways. First, light-flavor jets are initiated
either by a gluon or a quark, while heavy-flavor jets are mostly initiated by
quarks, therefore the average color charge of the initiating parton is
different. Then, the heavy-quark fragmentation functions differ from those of
light quarks and gluons [23]: heavy quark fragmentation is harder on the
average than light-quark fragmentation, meaning that a larger part of the
momentum is taken away by the heavy quarks [24]. Finally, heavy-flavor jet
shapes at low $p_{\mathrm{T}}$ are influenced by the dead cone effect, or the
suppression of forward gluon emissions [25].

In this work, we use PYTHIA 8 simulations to explore the connection between
the underlying event and the jet-initiating parton, using identified heavy and
light-flavor triggers. We argue that the connection of the underlying event to
the leading process depends on whether the initiating parton is a quark or a
gluon, and flavor-dependent studies are an excellent experimental tool to
understand this dependency. The article is organised as follows: Section 2
provides information about the analysis approach, as well as the simulations
using the PYTHIA 8 Monte Carlo generator. Section 3 presents the results and
discussion, and finally Section 4 summarizes the results.

##  2 Analysis method

In our simulations we used PYTHIA 8.240 [26] with the Monash 2013 tune [27]
and soft QCD settings. PYTHIA is an event generator that models a basic hard
scattering process with leading-order pQCD calculations, amended with initial-
state radiation (ISR), final-state radiation (FSR) to account for higher-order
processes. Multiple-parton interactions are also handled on the partonic
level, and CR is also applied at a later step. The hadronic final state is
produced using Lund string fragmentation, and then secondary decays and
rescattering between hadrons are computed. The Monash 2013 tune is mainly
focused on describing the minimum-bias and UE distributions accurately. We
simulated 100 million proton-proton events at $\sqrt{s}$ = 13 TeV with the
following settings.

To characterize the events by underlying-event activity, in each event we
identified the highest transverse momentum particle in the central
pseudorapidity window $|\eta|<0.8$ corresponding to the acceptance of the main
tracking detector of the ALICE experiment [28]. Only triggers above
$p_{\mathrm{T}}^{\rm trigger}=0.5$ GeV/$c$ were considered. We categorized
events based on the trigger particles: as for heavy-flavor, we collected
events where the triggers were either D-mesons ($\mathrm{D}^{0}$,
$\mathrm{\bar{D}}^{0}$ and $\mathrm{D}^{\pm}$) or B-mesons ($\mathrm{B}^{0}$,
$\mathrm{\bar{B}}^{0}$ and $\mathrm{B}^{\pm}$). For comparison, we identified
charged pions ($\pi^{\pm}$) to represent light flavor. Instead of a
reconstruction of heavy-flavor hadrons from the final state, we disabled their
decays. This way we also avoid accounting for the feed-down contribution from
B to D. Excited states, on the other hand, were allowed to decay, to simplify
the analogous experimental environment. Within each event we correlated all
associated particles to the trigger that had a minimal transverse momentum of
$p_{\mathrm{T}}^{\rm assoc.}$ = 0.5 GeV/$c$, at central pseudorapidity
($|\eta|<0.8$). We applied a spatial division based on the azimuth angle
($\Delta\varphi$) between the trigger and the associated particles. Following
the definition by the CDF Collaboration, the direction of the trigger particle
is used to define regions in the azimuthal plane that have different
sensitivity to the UE [29]. The azimuthal angular difference between the
trigger particles and the associated particles is used to define the following
azimuthal regions considered in this study: $|\Delta\varphi|<\frac{\pi}{3}$
for the toward region and $\frac{\pi}{3}<|\Delta\varphi|<\frac{2\pi}{3}$ for
the transverse region, dominated by the UE. PYTHIA 8 is known to reproduce the
charged-hadron multiplicities in the jet as well as in the underlying-event
regions [21]. Besides our reference simulations we also generated samples with
non-physical settings, to understand the individual effect of model
components: without FSR and ISR, disabled and enhanced CR.

##  3 Results and Discussion

We characterized the particle production both with the average density of
associated particles $\langle\mathrm{d}^{2}N_{\rm
ch}/\mathrm{d}\eta\mathrm{d}\varphi\rangle\equiv n_{\rm ch}^{\rm assoc.}$,
where $N_{\rm ch}$ is the number of particles in the selected region, and the
average transverse momentum sum density $\langle\mathrm{d}^{2}\sum{p}_{\rm
T}/\mathrm{d}\eta\mathrm{d}\varphi\rangle\equiv\sum p_{\mathrm{T}}^{\rm
assoc.}$ carried by these particles. These distributions are shown in Fig. 1
for both the toward and the transverse regions, and they are compared to
recent experimental data (open circles) for charged particles measured by the
ALICE Collaboration [21].

Figure 1:  Average charged-particle number (panels (a) and (b)) and
summed-$p_{\mathrm{T}}$ (panels (c) and (d)) densities as a function of
$p_{\mathrm{T}}^{\rm trigger}$ for charged pions, D- and B-meson triggers in
the toward (left column) and transverse (right column) regions. The PYTHIA 8
predictions are compared to experimental data (open circles) for charged
particles measured by the ALICE Collaboration [21].

In the toward region, a monotonic rise of both $n_{\rm ch}^{\rm assoc.}$ and
$\sum p_{\mathrm{T}}^{\rm assoc.}$ is present with $p_{\mathrm{T}}^{\rm
trigger}$. In case of $\sum p_{\mathrm{T}}^{\rm assoc.}$ this rise is
approximately linear, revealing that the trigger particle carries an
approximately constant fraction of the full jet momentum in the observed
momentum range. Since higher-momentum jets fragment into particles with higher
average $p_{\mathrm{T}}$, the steepness of the $n_{\rm ch}^{\rm assoc.}$ curve
decreases toward higher $p_{\mathrm{T}}^{\rm trigger}$. The ALICE data are
well described by simulations in which charged pions are the trigger
particles. However, if one considers heavy-flavour triggered events, the
activity is significantly reduced. This remarkable flavor hierarchy represents
the flavor (mass) dependent fragmentation. Similarly to what has been observed
for light flavor [21], in the transverse region there is a near-saturation
towards high $p_{\mathrm{T}}^{\rm trigger}$. We see that the saturation value
does not depend significantly on the flavor of the trigger. This indicates
that the leading process is largely independent of the UE. However, in the
intermediate transverse-momentum range $2<p_{\mathrm{T}}^{\rm trigger}<8$
GeV/$c$ the transverse event activity is different for each flavor: the
underlying event activity in pion-triggered events is higher than that for
heavy-flavor triggers. The effect is similar, albeit much larger, than what
has been found for different light-flavor hadrons [19], which shows that the
effect cannot be explained by differences in the the mass of the (final-state)
trigger hadron. In the same study, however, color reconnection (CR) was found
to be the reason for the hadron-quality-dependent behavior.

Figure 2:  Average number of multiple-parton interactions as a function of
$p_{\mathrm{T}}^{\rm trigger}$ for charged pions, D- and B-meson triggers
without and with color reconnection. Panels show different values of color
reconnection range ($RR$): $RR=0$, $RR=1.8$, $RR=10$ in the left, middle and
right panels, respectively.

The event activity in the transverse region is related to the number of multi-
parton interactions, $N_{\rm MPI}$, in an event [18]. In order to gain a
deeper insight, we investigate $\left<N_{\rm MPI}\right>$ in simulations with
different strengths of color reconnection, controlled in PYTHIA with the
Reconnection Range parameter ($RR$). Figure 2 shows the average number of
multiple-parton interactions per trigger, $\left<N_{\rm MPI}\right>$, as a
function of $p_{\mathrm{T}}^{\rm trigger}$, for the different identified
triggers. We see that if there is no color reconnection ($RR=0$), then there
is no significant difference between $N_{\rm MPI}$ for any of the trigger
hadron flavors. If the color reconnection matches the physical tune
($RR=1.8$), a clear difference can be observed in $\left<N_{\rm MPI}\right>$
depending on the trigger particle species in the momentum range
$2<p_{\mathrm{T}}^{\rm trigger}<8$ GeV/$c$, similarly to what is present for
$n_{\rm ch}^{\rm assoc.}$ and $\sum p_{\mathrm{T}}^{\rm assoc.}$. The MPI
activity corresponding to light-flavor triggers is enhanced in this region
compared to heavy-flavor triggered events, forming a visible bump structure.
This bump is already much less prominent in the case of D-mesons, and for
B-mesons the MPI activity grows monotonically until saturation. As can be
expected, the effect is even more pronounced in the enhanced color-
reconnection scenario. It is noteworthy that in the toward region the flavor
ordering is present even if there is no color reconnection. Therefore, it can
be assumed that the ordering in the transverse region arises in connection
with the ordering observed in the jet region, via final-state collective
traits produced by CR [15].

In the string fragmentation model of PYTHIA, the main difference between quark
and gluon hadronization is that gluons are connected to two string pieces,
while quarks are only connected to a single string piece [30]. One would
expect that selecting gluon jets biases the sample towards events with a
larger number of strings than selecting enhanced quark jet samples. Given the
larger number of strings in events with gluon jets than with quark jets, color
reconnection is expected to give a larger flow-like effect if enhanced gluon
jet samples are analyzed. To understand whether what we observed using pion or
B meson trigger particles is connected with this, we estimate the fraction of
gluon (quark) jets in triggered leading-pion (-B meson) events.

In order to estimate the gluon (quark) jet contribution in event samples
containing a pion (B meson) trigger particle, we implemented an algorithm that
identifies events based on whether the parton that eventually fragments into
the trigger particle is a quark or a gluon. The initiating parton is not
necessarily well-defined in a complex interaction. In PYTHIA there is no a
priori way to determine the initiating parton type with full certainty. As an
approximation, we applied the following method: we look for a
high-$p_{\mathrm{T}}$ parton in a cone region around the trigger particle,
with a radius $R=\sqrt{\Delta\varphi^{2}+\Delta\eta^{2}}<0.5$. Since the
trigger particle is most likely a jet constituent and it is expected to be
close to the jet axis, this region is sensitive to jet fragmentation. The
value of $R$ ensures that the dominant contribution in this region comes from
the jet, and the jet will be mostly contained in this cone [31]. To prevent
the association of soft partons to the trigger particle, the presence of a
high-$p_{\mathrm{T}}$ parton with a $p_{\mathrm{T}}$ of at least $30\%$ of the
$p_{\mathrm{T}}$ of the trigger particle was required. Moreover, the selected
parton should be marked “ancestor” by PYTHIA. If no such parton is found, then
we look for the closest parton in the event record that has
$p_{\mathrm{T}}>0.3\cdot p_{\mathrm{T}}^{\rm trigger}$ and is within the
$R<0.5$ cone. If no such parton exists, then we finally look for the closest
parton in the event record that is marked “ancestor” by PYTHIA. As a next step
we assess the effectivity of our sorting algorithm. Figure 3 shows the
fractions of charged pion, D-meson and B-meson triggers that are identified as
gluon-initiated to the case without the selection of the initiating parton.

Figure 3:  Ratio of $p_{\mathrm{T}}^{\rm trigger}$ distributions with and
without the selection on gluon-initiated triggers for charged pions, D-, and
B-mesons. The left and the right panel shows the cases with and without ISR
and FSR.

The left panel shows the physical case, while in the right panel the case is
shown where initial-state radiation (ISR) and final-state radiation (FSR) are
turned off. Since in this case no radiated gluons are present, we expect that
our selection algorithm works reliably. In the case where ISR and FSR are
allowed, the algorithm identifies up to 30% of D-meson triggers and 10–15% of
B-meson triggers as gluon-initiated. On the other hand, without radiation
contribution the fraction of heavy-flavor triggers identified as gluon-
initiated is below 5–10% at low-$p_{\mathrm{T}}^{\rm trigger}$ and vanishes
towards higher $p_{\mathrm{T}}^{\rm trigger}$. Regarding pion triggers, the
algorithm identifies more than 60% of the events as gluon-initiated in the
$0.5<p_{\mathrm{T}}^{\rm trigger}<5$ GeV/$c$ interval. As a conclusion, B-jets
can be used as proxies for quark jets within an error margin below 5–10%, and
pions can be used as proxies for gluon jets albeit with a higher
misidentification rate, between 10-25% depending on the observed
$p_{\mathrm{T}}^{\rm trigger}$ range. It is to be noted that higher-order hard
processes, such as gluon splitting and flavor excitation, comprise a
significant contribution to heavy flavor production [32]. Therefore the parton
that initiates the triggering jet cannot be directly identified to the hardest
pQCD process. Nevertheless, we also cross-checked our algorithm with PYTHIA
where gluon and quark processes were explicitly selected first as a PYTHIA
hard QCD process, and also as the first scattering process in soft QCD. We
found similar patterns in both cases to our initiating parton selection
algorithm. However, a larger fraction (up to 35%) of heavy-flavor triggers
were recorded to come from events with a gluon leading process.

In order to see how well the flavor separation might work in an experimental
environment, we compare the behavior of B-mesons to pions, as proxies for
quark-initiated and gluon-initiated jets, for different CR strengths. In Fig.
4 we show $\left<N_{\rm MPI}\right>$ as a function of $p_{\mathrm{T}}^{\rm
trigger}$ for B-mesons and charged pions. For pion-triggered events, the
increase of the color reconnection range causes a depletion of the
$\left<N_{\rm MPI}\right>$ at low $p_{\mathrm{T}}^{\rm trigger}$
($p_{\mathrm{T}}^{\rm trigger}\approx 3$ GeV/$c$) and an enhancement at high
$p_{\mathrm{T}}^{\rm trigger}$ ($p_{\mathrm{T}}^{\rm trigger}\approx 5$
GeV/$c$). The bump is not observed for B-meson triggered events.

Figure 4:  Average number of multiple-parton interactions for B-meson (left
panel) and charged pion (right panel) triggers for different values of color
reconnection ranges ($RR$): $RR=0$, $RR=1.8$, and $RR=10$.

For the purpose of providing an experimentally accessible quantity, we show
$\sum p_{\mathrm{T}}^{\rm assoc.}$ for different CR settings in a similar
manner in Fig. 5. Our simulation results are also compared to recent
experimental data (full circles) for charged particles measured by the ALICE
Collaboration [21]. Note that while the pattern seen in $\left<N_{\rm
MPI}\right>$ (in Fig. 2) is somewhat better represented in the transverse
particle density than in the sum transverse momentum (shown in Fig. 1), the
presence of color reconnection is known to affect event multiplicity
distributions [16] and therefore the number of particles cannot be used when
comparing different CR settings, while the momentum is conserved.

Figure 5:  Average summed-$p_{\mathrm{T}}$ densities for B-meson (left panel)
and charged pion (right panel) triggers for different values of color
reconnection ranges ($RR$): $RR=0$, $RR=1.8$, and $RR=10$. The PYTHIA 8
predictions are compared to experimental data (open circles) for charged
particles measured by the ALICE Collaboration [21].

While in the case of summed-$p_{\mathrm{T}}$ a weak but steady rise is present
in the range where $\left<N_{\rm MPI}\right>$ saturates, the observed trends
are again very similar. There is no significant difference above the point of
onset for the $\pi^{\pm}$ case, while in the semi-hard $p_{\mathrm{T}}^{\rm
trigger}$ range the summed-$p_{\mathrm{T}}$ distributions of light $\pi^{\pm}$
mesons show a characteristic dependence on CR. These observations corroborate
the assumption that the leading hard process is connected to the UE via color
reconnection and that the flavor quality of the leading process has a strong
influence on the development of the UE. The observed pattern suggests that the
flavor dependence is driven in a large part by the color composition of the
initiating parton. Experimental confirmation of the results presented here can
be useful to improve our knowledge on the particle production mechanism in pp
interactions, and their connection with the heavy-ion-like collective effects
discovered at the LHC.

##  4 Conclusions

We investigated the particle production in the underlying event corresponding
to identified light and heavy-flavor triggers, using PYTHIA 8 simulations. The
particle yield in the transverse region, as well as the number of multiple-
parton interactions that is strongly correlated to it, saturates towards
higher $p_{\mathrm{T}}^{\rm trigger}$, independently of particle species.
However, as a consequence of color reconnection, flavor hierarchy shows up in
the region $p_{\mathrm{T}}^{\rm trigger}<8$ GeV/$c$ where soft processes play
an important role. It can be assumed that this arises from the jet region
being linked to the transverse region via collectivity generated by color
reconnection. In the case of light flavor, $\left<N_{\rm MPI}\right>$
overshoots the saturation value, while this is not seen in heavy flavor. This
“bump” structure observed at $p_{\mathrm{T}}^{\rm trigger}\approx 4$ GeV/$c$
can be associated with gluon-initiated jets. We find that while pion triggers
are mostly gluon-initiated, B-mesons serve as a good proxy for quark jets.
This provides an opportunity to evaluate the connection of gluon- and quark-
initiated jets to the underlying event by comparing the particle production in
the transverse region by light-hadron and B-meson triggered events in the
experiment, and thus gives new means for model development.

## Acknowledgement

This work was supported by the Hungarian National Research, Development and
Innovation Office (NKFIH) under the contract numbers OTKA FK131979 and
K135515, and the NKFIH grant 2019-2.1.11-TÉT-2019-00078. A. Ortiz acknowledges
the support from CONACyT under the Grants CB No. A1-S-22917 and CF No. 2042.
The authors also acknowledge the computational resources provided by the
Wigner GPU Laboratory and research infrastructure provided by the Eötvös
Loránd Research Network (ELKH).

## References

## References

  * [1] K. Adcox, et al., Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration, Nucl. Phys. A 757 (2005) 184–283.  arXiv:nucl-ex/0410003, doi:10.1016/j.nuclphysa.2005.03.086.
  * [2] J. Adams, et al., Experimental and theoretical challenges in the search for the quark gluon plasma: The STAR Collaboration’s critical assessment of the evidence from RHIC collisions, Nucl. Phys. A 757 (2005) 102–183.  arXiv:nucl-ex/0501009, doi:10.1016/j.nuclphysa.2005.03.085.
  * [3] B. B. Back, et al., The PHOBOS perspective on discoveries at RHIC, Nucl. Phys. A 757 (2005) 28–101.  arXiv:nucl-ex/0410022, doi:10.1016/j.nuclphysa.2005.03.084.
  * [4] I. Arsene, et al., Quark gluon plasma and color glass condensate at RHIC? The Perspective from the BRAHMS experiment, Nucl. Phys. A 757 (2005) 1–27.  arXiv:nucl-ex/0410020, doi:10.1016/j.nuclphysa.2005.02.130.
  * [5] V. Khachatryan, et al., Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC, JHEP 09 (2010) 091.  arXiv:1009.4122, doi:10.1007/JHEP09(2010)091.
  * [6] S. Acharya, et al., Investigations of Anisotropic Flow Using Multiparticle Azimuthal Correlations in pp, p-Pb, Xe-Xe, and Pb-Pb Collisions at the LHC, Phys. Rev. Lett. 123 (14) (2019) 142301.  arXiv:1903.01790, doi:10.1103/PhysRevLett.123.142301.
  * [7] C. Bierlich, Microscopic collectivity: The ridge and strangeness enhancement from string–string interactions, Nucl. Phys. A 982 (2019) 499–502.  arXiv:1807.05271, doi:10.1016/j.nuclphysa.2018.07.015.
  * [8] C. Bierlich, G. Gustafson, L. Lönnblad, Collectivity without plasma in hadronic collisions, Phys. Lett. B 779 (2018) 58–63.  arXiv:1710.09725, doi:10.1016/j.physletb.2018.01.069.
  * [9] B. Blok, C. D. Jäkel, M. Strikman, U. A. Wiedemann, Collectivity from interference, JHEP 12 (2017) 074.  arXiv:1708.08241, doi:10.1007/JHEP12(2017)074.
  * [10] L. Yan, J.-Y. Ollitrault, Universal fluctuation-driven eccentricities in proton-proton, proton-nucleus and nucleus-nucleus collisions, Phys. Rev. Lett. 112 (2014) 082301.  arXiv:1312.6555, doi:10.1103/PhysRevLett.112.082301.
  * [11] T. A. Trainor, Questioning quark-gluon plasma formation in small collision systemsarXiv:1905.11585.
  * [12] A. Ortiz, E. Zepeda, Extraction of the multiplicity dependence of Multiparton Interactions from LHC pp data using Machine Learning techniquesarXiv:2101.10274, doi:10.1088/1361-6471/abef1e.
  * [13] A. Ortiz, A. Paz, J. D. Romo, S. Tripathy, E. A. Zepeda, I. Bautista, Multiparton interactions in $pp$ collisions from machine learning-based regression, Phys. Rev. D 102 (7) (2020) 076014.  arXiv:2004.03800, doi:10.1103/PhysRevD.102.076014.
  * [14] T. Sjöstrand, The Development of MPI Modeling in Pythia, Adv. Ser. Direct. High Energy Phys. 29 (2018) 191–225.  arXiv:1706.02166, doi:10.1142/9789813227767_0010.
  * [15] A. Ortiz, G. Bencedi, H. Bello, Revealing the source of the radial flow patterns in proton–proton collisions using hard probes, J. Phys. G 44 (6) (2017) 065001.  arXiv:1608.04784, doi:10.1088/1361-6471/aa6594.
  * [16] Z. Varga, R. Vértesi, G. Gábor Barnaföldi, Modification of jet structure in high-multiplicity pp collisions due to multiple-parton interactions and observing a multiplicity-independent characteristic jet size, Adv. High Energy Phys. 2019 (2019) 6731362.  arXiv:1805.03101, doi:10.1155/2019/6731362.
  * [17] A. Ortiz Velasquez, P. Christiansen, E. Cuautle Flores, I. Maldonado Cervantes, G. Paić, Color Reconnection and Flowlike Patterns in $pp$ Collisions, Phys. Rev. Lett. 111 (4) (2013) 042001.  arXiv:1303.6326, doi:10.1103/PhysRevLett.111.042001.
  * [18] T. Martin, P. Skands, S. Farrington, Probing Collective Effects in Hadronisation with the Extremes of the Underlying Event, Eur. Phys. J. C 76 (5) (2016) 299.  arXiv:1603.05298, doi:10.1140/epjc/s10052-016-4135-4.
  * [19] A. Ortiz, L. Valencia Palomo, Probing color reconnection with underlying event observables at the LHC energies, Phys. Rev. D 99 (3) (2019) 034027.  arXiv:1809.01744, doi:10.1103/PhysRevD.99.034027.
  * [20] A. Ortiz, L. Valencia Palomo, Universality of the underlying event in pp collisions, Phys. Rev. D 96 (11) (2017) 114019.  arXiv:1710.04741, doi:10.1103/PhysRevD.96.114019.
  * [21] S. Acharya, et al., Underlying Event properties in pp collisions at $\sqrt{s}$ = 13 TeV, JHEP 04 (2020) 192.  arXiv:1910.14400, doi:10.1007/JHEP04(2020)192.
  * [22] S. Tripathy, Topological studies of light-flavor hadron production in high multiplicity pp collisions with ALICE at the LHC, in: 40th International Conference on High Energy Physics, 2020.  arXiv:2011.00741.
  * [23] L. Dai, C. Kim, A. K. Leibovich, Heavy Quark Jet Fragmentation, JHEP 09 (2018) 109.  arXiv:1805.06014, doi:10.1007/JHEP09(2018)109.
  * [24] R. Vértesi, Jet measurements with ALICE: substructure, dead cone, charm jets, PoS LHCP2020 (2021) 143.  arXiv:2009.11228, doi:10.22323/1.382.0143.
  * [25] Y. L. Dokshitzer, V. A. Khoze, S. I. Troian, On specific QCD properties of heavy quark fragmentation (’dead cone’), J. Phys. G 17 (1991) 1602–1604.  doi:10.1088/0954-3899/17/10/023.
  * [26] T. Sjöstrand, S. Mrenna, P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178 (2008) 852–867.  arXiv:0710.3820, doi:10.1016/j.cpc.2008.01.036.
  * [27] P. Skands, S. Carrazza, J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C 74 (8) (2014) 3024.  arXiv:1404.5630, doi:10.1140/epjc/s10052-014-3024-y.
  * [28] C. W. Fabjan, et al., ALICE: Physics performance report, volume II, J. Phys. G 32 (2006) 1295–2040.  doi:10.1088/0954-3899/32/10/001.
  * [29] C. M. Buttar, et al., The Underlying Event, in: HERA and the LHC: A Workshop on the Implications of HERA for LHC Physics: CERN - DESY Workshop 2004/2005 (Midterm Meeting, CERN, 11-13 October 2004; Final Meeting, DESY, 17-21 January 2005), CERN, Geneva, 2005.
  * [30] T. Sjöstrand, S. Mrenna, P. Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026.  arXiv:hep-ph/0603175, doi:10.1088/1126-6708/2006/05/026.
  * [31] S. Chatrchyan, et al., Shape, Transverse Size, and Charged Hadron Multiplicity of Jets in pp Collisions at 7 TeV, JHEP 06 (2012) 160.  arXiv:1204.3170, doi:10.1007/JHEP06(2012)160.
  * [32] E. Norrbin, T. Sjöstrand, Production and hadronization of heavy quarks, Eur. Phys. J. C 17 (2000) 137–161.  arXiv:hep-ph/0005110, doi:10.1007/s100520000460.


# CR-ENTREES - Cosmic-Ray ENergy TRansport in timE-Evolving astrophysical
Settings

Anita Reimer  anita.reimer(at)uibk.ac.at Lukas Merten  Margot Boughelilba
Paolo Da Vela  Serguei Vorobiov  Jon Paul Lundquist

###### Abstract

In order to understand observable signatures from putative cosmic-ray (CR)
sources in-source acceleration of particles, their energy and time-dependent
transport including interactions in an evolving environment and their escape
from source have to be considered, in addition to source-to-Earth propagation.

We present the code CR-ENTREES (Cosmic-Ray ENergy TRansport in timE-Evolving
astrophysical Settings) that evolves the coupled time- and energy-dependent
kinetic equations for cosmic-ray nucleons, pions, muons, electrons, positrons,
photons and neutrinos in a one-zone setup of (possibly) non-constant size,
with user-defined particle and photon injection laws. All relevant
interactions, particle/photon escape and adiabatic losses are considered in a
radiation-dominated, magnetized astrophysical environment that is itself
evolving in time. Particle and photon interactions are pre-calculated using
event generators assuring an accurate interactions and secondary particle
production description. We use the matrix multiplication method for fast
radiation and particle energy transport which allows also an efficient
treatment of transport non-linearities due to the produced particles/photons
being fed back into the simulation chain.

Examples for the temporal evolution of the non-thermal emission from AGN jet-
like systems with focus on proton-initiated pair cascades inside an expanding
versus straight jet emission region, are further presented.

††journal: PoS: 38th ICRC 2023

[a]organization=Universitaet Innsbruck, Institut fuer Astro- und
Teilchenphysik,city=Innsbruck, postcode=6020, country=Austria
[b]organization=Theoretical Physics IV, Plasma Astroparticle Physics, Faculty
for Physics and Astronomy, Ruhr University Bochum,city=Bochum, postcode=44780,
country=Germany [c]organization=Ruhr Astroparticle and Plasma Physics Center
(RAPP Center),city=Bochum, postcode=44780, country=Germany
[d]organization=Center for Astrophysics and Cosmology (CAC), University of
Nova Gorica,city=Nova Gorica, postcode=SI-5000, country=Slovenia

##  1 Introduction

With more than 60% of all the detected sources in the $\gamma$-ray sky
belonging to the class of jetted active galactic nuclei (AGN), and their
proposed contribution to the PeV neutrino and ultra-high energy cosmic-ray
(UHECR) sky, a deep exploration of the multi-messenger nature of this source
class is central. Where and how are the high-energy messengers produced? How
are the charged particles accelerated to such extreme energies and how can
they escape the magnetized jet environment? What is the overall composition of
the jet? A tool that supports the investigation of these questions is
presented here: CR-ENTREES, a code for fully time-dependent Cosmic-Ray ENergy
TRansport in timE-Evolving astrophysical Settings.

CR-ENTREES is used as the base code for the heavy nuclei propagation code of
[9] (see this proceedings).

##  2 Propagation Physics

CR-ENTREES solves the following nonlinear system of coupled (integro-
differential) Fokker-Planck transport equations:

$\partial_{t}F_{N}+\dot{F}^{\rm esc}_{N}+\partial_{E}[(\dot{E}_{\rm
loss}F_{N})]+\dot{F}^{\rm dec}_{N}=Q^{\rm inj,pr}_{N}$
$\partial_{t}F_{\mu,\pi,K}+\dot{F}^{\rm
esc}_{\mu,\pi,K}+\partial_{E}[(\dot{E}_{\rm loss}F_{\mu,\pi,K})]+\dot{F}^{\rm
dec}_{\mu,\pi,K}=\dot{F}^{p\gamma;h}_{\mu,\pi,K}$
$\partial_{t}F_{e}+\dot{F}^{\rm esc}_{e}+\partial_{E}[(\dot{E}_{\rm
loss}F_{e})]=Q^{\rm
inj,pr}_{e}+\dot{F}_{e}^{\gamma\gamma}+\dot{F}_{e}^{p\gamma}$
$\partial_{t}F_{\gamma}+\dot{F}^{\rm
esc}_{\gamma}+\dot{F}_{\gamma}^{\gamma\gamma}=\dot{F}_{\gamma}^{\rm
em}+\dot{F}_{\gamma}^{p\gamma;h}$

Here, $F_{\rm X}=F_{\rm X}(E,t)$ is the energy (E) and time (t) dependent
density of the cosmic-ray (CR) nucleons $X=N$, electrons and positrons $X=e$,
muons, pions, kaons $X=\mu,\pi,K$, respectively, and photons $X=\gamma$.
$Q^{\rm inj,pr}_{N,e}=Q^{\rm inj,pr}_{N,e}(E,t)$ describes the source function
of the primary nucleons (N) and pairs (e), and $\dot{E}_{\rm
loss}=\dot{E}_{\rm loss}(F_{\gamma}(\epsilon,t),B(t);E,t)$ the continuous loss
processes that potentially depend on the (possibly evolving) magnetic field
strength $B(t)$ and a energy ($\epsilon$) and time dependent target radiation
field. Injection of secondary particles and photons produced in hadronic (h)
and electromagnetic (em) particle-photon, photon-photon and particle-field
interactions are considered with the density rates
$\dot{F}_{e}^{p\gamma}=\dot{F}_{e}^{p\gamma}(F_{\gamma}(\epsilon,t);E,t)$,
$\dot{F}^{p\gamma;h}_{\mu,\pi,K}$,
$\dot{F}_{e}^{\gamma\gamma}=\dot{F}_{e}^{\gamma\gamma}(F_{\gamma}(\epsilon,t);E,t)$,
$\dot{F}_{\gamma}^{\gamma\gamma}=\dot{F}_{\gamma}^{\gamma\gamma}(F_{\gamma}(\epsilon,t);\epsilon,t)$,
$\dot{F}_{\gamma}^{p\gamma}=\dot{F}_{\gamma}^{p\gamma}(F_{\gamma}(\epsilon,t);\epsilon,t)$.
The rigorous treatment of a time-evolving emission volume, magnetic field and
target radiation field within the present setup causes the non-linearity of
this transport equation system. The user of CR-ENTREES can choose to consider
the continuous loss processes inverse Compton scattering, Bethe-Heitler pair
production and synchrotron radiation of all charged particles, and the
catastrophic loss processes photomeson production, particle decay and escape.
For a conical jet setup adiabatic losses are taken into account as well. We
note that in the absence of nucleons the system becomes identical to a
Synchrotron-Self Compton model.

##  3 The Model

While CR-ENTREES (thanks to its modular implementation) is flexible to be
adapted to a large range of geometrical setups where CR transport occurs, its
first purpose is to describe a moving homogeneous one-zone emission region
within a straight or conical outflow/jet. Here, CR transport is treated in the
co-moving frame. Hence, all input parameters (except the outflow speed - see
below) are considered in this frame.

###  3.1 Input: Characterizing geometry and environment

The moving (with given constant speed $\beta_{\rm J}c$) spherical emission
region of fixed or time-evolving (in case of a conical outflow) radius $R$
contains a homogeneous magnetic field of strength $B$ and (assumed) isotropic
target radiation field for particle-photon interactions. Magnetic field
evolution is currently implemented with a $B(t)\propto t^{-1}$ scaling (with
$t$ the co-moving propagation time) in case of a conical jet setup. The target
photon field density distribution is discretized on a fixed 161 log-equal
spacing energy grid in the range $10^{-10...6}$eV. A (diluted) blackbody with
given temperature $T$ or power-law spectra (normalized using a given energy
density) with up to 2 break energies (and corresponding power-law indices) can
be chosen to fill this target field. Alternatively, this field can be filled
for each energy bin by the user. The total target radiation field is then
determined from the internal (jet) radiation field (calculated in each time
step) and the user-defined target field.

###  3.2 Input: Characterizing particle injection

Power-law spectra (of possibly exponential cutoff) of any particle type (or
photon) within a user-defined energy range and index are injected into the
simulation chain. These injection spectra (and all particle spectra during
propagation) are discretized on a fixed 300 log-equal spacing energy grid in
the range $10^{-3...12}$GeV. Currently, up to 2 particle populations can be
injected. For their normalization the number ratio of these populations and
the particle-to-field energy density ratio must be provided by the user.

All input parameter values are provided by the user in a dedicated steering
file.

###  3.3 Energy loss processes and secondary particle production

We use Monte Carlo event generators (photomeson production and decay
processes: [1]; Bethe-Heitler pair production: modified version of [2];
inverse Compton scattering and photon-photon pair production: [3]) to pre-
calculate the yields and interaction rates (which are stored in HDF5-files) of
each interaction type, and discretized on a fixed 300 log-equal spacing energy
grid in the range $10^{-3...12}$GeV for a range of energies for the target
particle/photon. The corresponding yield and interaction rate for the (assumed
isotropic) target radiation field in each time step are then determined by
convolving over this target field. Synchrotron radiation yields and
corresponding loss rates (using pitch-angle averaged terms) of all produced
charged particles (i.e., charged pions, muons, kaons and electrons/positrons)
are calculated following [4], [5], and [6] (for the self-absorption process).
Note that particles suffering from continuous losses in combination with
catastrophic losses on a fixed energy grid require a dedicated numerical
treatment to let these particles move down this energy grid from high to low
energies.

The adiabatic loss time scale for particles of Lorentz factor $\gamma$ in a
conical jet of opening angle $0.26/\Gamma_{\rm J}$ [7], with $\Gamma_{\rm
J}=(1-\beta_{\rm J}^{2})^{-0.5}$ the bulk Lorentz factor of the moving
emission region, is calculated to $t_{\rm
ad}=\frac{\gamma^{2}}{\gamma^{2}-1}t\left(1-\frac{R_{0}}{R(t)}\right)^{-1}$
with $R_{0}$ the size of the emission region at propagation time $t=0$. Note
the decrease of the adiabatic loss rate with increasing propagation time.

Finally, neutral particles (including photons) escape on a time scale $t_{\rm
esc,n}=\frac{3}{4}R(t)/c$, while for the charged particles’ escape time scale
$t_{\rm esc,c}=\eta t_{\rm esc,n}$ with $\eta\geq 1$, is used.

All particle energy, adiabatic and escape losses are tracked as well
throughout the entire simulation chain to allow for the verification of energy
conservation in each simulation time step.

###  3.4 Propagation method

We propagate $\gamma$-ray, proton, neutron, electronic pair, muon, pion, kaon,
muon and electron neutrino populations, discretized on a fixed 300 log-equal
spacing energy grid in the range $10^{-3...12}$GeV, and lower energy photons
on a fixed 300 log-equal spacing energy grid in the range $10^{-18...-3}$GeV,
using the matrix multiplication method of [3], [8], [2] in the framework of CR
transport. Here, transfer matrices are created from the aforementioned yields
and interaction probabilities which describe the change after a given time
step $\delta t$ of the density of a given particle type upon all the
interactions pre-set by the user. Such explicit integration scheme, while
extremely fast, requires to use time steps that are smaller than the smallest
time scale of the system (the Courant-Friedrichs-Levy condition). Still, with
the matrix doubling method of [8] applied, we found correct and stable results
also for somewhat larger time steps. Here, energy conservation is verified
after each time step.

###  3.5 Code output

CR-ENTREES’ output encompasses the co-moving density of all propagated
particle types after the pre-set propagation time, and interaction rates of
all chosen processes, on the aforementioned energy grid. Post-processing of
the output density within the emission volume to transform into the observer
frame is carried out outside the CR-ENTREES framework.

##  4 Examples

As an example for the use of CR-ENTREES we show here the development of a pair
cascade initiated by a cosmic-ray proton population inside a magnetized,
relativistically moving emission region along a straight, and a conical jet.
In our example model the initial size of the comoving emission region amounts
to $R_{\rm emi}=3\times 10^{16}$cm, located at a distance of $10^{17}$cm away
from the central engine, and moves with a speed that corresponds to a bulk
Lorentz factor of $\Gamma_{\rm J}=22.4=D$ (where $D$ is the Doppler factor).
We instantaneously inject a relativistic electron population, that follows a
$\propto\gamma_{e}^{-2}$ particle energy distribution between Lorentz factors
$\gamma_{e}=1$ and $\gamma_{e}=100$, into the magnetized (with initial
magnetic field strength of $10$G) emission region. These electrons build up a
synchrotrotron radiation field, which has its power peak in the X-ray domain
(see Figs. 1, 2), resembling an HBL-like synchrotron spectrum to some extend.
Cosmic-ray protons, injected along with the electron population with the same
spectral index but in an extended particle energy range from $\gamma_{p}=1$ to
$\gamma_{p}=10^{9}$, initiate then, via proton-photon interactions in this
evolving synchrotron radiation field, a pair cascade, which we follow with
time steps of $\delta t=10^{4}$s. With an initial particle-to-field energy
density ratio of $10$, and a proton-to-electron injection number density ratio
of $10^{5}$, the normalization of the two injected particle populations are
set. In order to focus on the cascade development, we switch off proton
synchrotron radiation in these simulations. All other processes, including
charged particle escape (on a time scale $t_{\rm esc,c}=3/4R_{\rm emi}/c$) are
kept active.

Figure 1 shows the evolution of the broadband photon SED from a straight jet
moving emission region in the galaxy frame at times $t_{\rm obs}\cdot
D=10^{4}$s, $2\times 10^{4}$s, $3\times 10^{4}$s, $10^{5}$s, $2\times
10^{5}$s, $8\times 10^{5}$s (black lines) and $t_{\rm obs}\cdot D=1.8\times
10^{6}$s, $2.8\times 10^{6}$s, $5.9\times 10^{6}$s, $8.9\times 10^{6}$s (grey
lines) after the injection. The synchrotron radiation field, target for
particle-photon interactions, reaches its maximum power within less than one
dynamical time scale, its decline is much slower. Figure 3 shows the
correspondingly strongly asymmetric X-ray light curve (dashed red line). Also
the subsequently developing pair cascade responds with a fast flux increase,
and a slower decrease past its power maximum (see the $\gamma$-ray light
curves in Fig. 3). The 1 PeV neutrino light curve (violet line in Fig. 3)
responds even slower: Its power maximum is reached within less than a day
(galaxy frame), its decline starts clearly after the decrease of the photon
cascade’s emission.

We then study such proton-initiated pair cascade development, using the same
model parameters as above, in a moving emission region within now a conical
jet of opening angle $\sim 0.66^{\rm o}$. The adiabatic expansion of the
emission region leads in this case to further energy losses. At the same time,
the corresponding increase of the comoving escape time scales (from $7.5\times
10^{5}$s at injection time to $3.7\times 10^{6}$s at $1.5\times 10^{7}$s after
injection at the end of our simulations) as well as the decline of the
magnetic field $\propto t^{-1}$ keeps particles and photons inside the
emission region and decreases radiative losses. As a result, the increase and
decrease of the target and cascade photon emission is significantly prolonged,
as can be appreciated from Figs. 2, 3 of our work: E.g., the duration of the
flux enhancement above one order of magnitude below the respective power
maximum in the 3 energy bands shown seems a factor of 2-3 longer for the
expanding conical jet as compared to the straight jet.

From this study we conclude: The temporal evolution of the size of the
emission region, e.g., as a result of the jet’s shape, has a significant
impact on the resulting multi-messenger light curves.

Figure 1: Photon SED (host galaxy frame) from the emission region moving along
a straight jet with parameters as described in Sect. 4, at (comoving) times
$t=10^{4}$s (solid black line), $2\times 10^{4}$s (dotted black line),
$3\times 10^{4}$se (dashed black line), $10^{5}$s (dashed-dotted black line),
$2\times 10^{5}$s (dashed-triple-dotted black line), $8\times 10^{5}$s (long-
dashed black lines) and $t=1.8\times 10^{6}$s (solid grey line), $2.8\times
10^{6}$s (dotted grey line), $5.9\times 10^{6}$s (dashed grey line),
$8.9\times 10^{6}$s (dashed-dotted grey line) after injection. The red dashed
line indicates the 1 keV-slice, the blue dashed line the 100 MeV-slice and the
pale blue dashed line the 1 TeV-slice through the flux-energy-time data cube.
Figure 2: Photon SED (host galaxy frame) from the emission region moving along
a conical jet with parameters as described in Sect. 4, at (comoving) times
$t=10^{4}$s (solid black line), $2\times 10^{4}$s (dotted black line),
$3\times 10^{4}$s (dashed black line), $10^{5}$s (dashed-dotted black line),
$2\times 10^{5}$s (dashed-triple-dotted black line), $8\times 10^{5}$s (long-
dashed black lines) and $t=1.8\times 10^{6}$s (solid grey line), $2.8\times
10^{6}$s (dotted grey line), $5.9\times 10^{6}$s (dashed grey line),
$8.9\times 10^{6}$s (dashed-dotted grey line), $1.2\times 10^{7}$s (dashed-
triple-dotted grey line), $1.5\times 10^{7}$s (long-dashed grey line) after
injection. The red dashed line indicates the 1 keV-slice, the blue dashed line
the 100 MeV-slice and the pale blue dashed line the 1 TeV-slice through the
flux-energy-time data cube. Figure 3: Light curves (galaxy frame) from the
straight jet (dashed lines) and conical jet (solid lines) emission region
taken at photon energies of 1 keV (red lines), 100 MeV (blue lines, 1 TeV
(pale blue lines) and at neutrino energy 1 PeV (violet lines).

##  5 Conclusion

We present CR-ENTREES (Cosmic-Ray ENergy TRansport in timE-Evolving
astrophysical Settings), a cosmic-ray energy transport code that solves the
nonlinear system of coupled Fokker-Planck transport equations for cosmic-ray
nucleons, mesons, leptons, photons and neutrinos. Its flexibility allows to
study the resulting multi-messenger emission from particle and photon energy
propagation in an environment that is itself changing with time.

In this work CR-ENTREES is used to study the temporal evolution of pair
cascades, initiated by the instantaneous injection of relativistic electrons
and protons, in an expanding jet, as compared to a straight jet. We find that
the cascade development in an expanding emission region, accompanied with a
diluted field environment, leads to a significantly prolonged outburst of
multi-messenger emission, as compared to the emission from a fixed-size region
in a straight jet.

## Acknowledgements

Financial support for this project was received from the Austrian Science Fund
(FWF) under grant agreement number I 4144-N27 and the Slovenian Research
Agency-ARRS (project no. N1-0111). MB has for this project received funding
from the European Union’s Horizon 2020 research and innovation program under
the Marie Sklodowska-Curie grant agreement No 847476. The views and opinions
expressed herein do not necessarily reflect those of the European Commission.
LM acknowledges support from the DFG within the collaborative Research Center
SFB 1191 ”Cosmic Interacting Matters - From Source to Signal”.

## References

  * [1] A. Mücke, R. Engel, J.P. Rachen, R.J. Protheroe, T. Stanev, _Monte Carlo simulations of photohadronic processes in astrophysics_ , _CPC_ 124 (2000) 290
  * [2] R.J. Protheroe, P.A. Johnson _Propagation of ultra high energy protons and gamma rays over cosmological distances and implications for topological defect models_ _Astroparticle Physics_ 5 (1996) 215
  * [3] R.J. Protheroe, _Effect of electron-photon cascading on the observed energy spectra ofextragalactic sources of ultra-high-energy gamma-rays_ _MNRAS_ 221 (1986) 769
  * [4] G.B. Pacholczyk, _Radio Astrophysics_ Freeman, San Francisco, 1970.
  * [5] J.J. Brainerd, V. Petrosian _Analytic and Numerical Calculations of Quantum Synchrotron Spectra from Relativistic Electron Distributions_ _ApJ_ 320 (1987) 703
  * [6] G.B. Rybicki, A.P. Lightman _Radiative Processes in Astrophysics_ John Wiley & Sons, 2017.
  * [7] A.B. Pushkarev, et al, _MOJAVE - XIV. Shapes and opening angles of AGN jets_ _MNRAS_ 468 (2017) 4992
  * [8] R.J. Protheroe, T. Stanev _Electron-photon cascading of very high-energy gamma-rays in the infrared background_ _MNRAS_ 264 (1993) 191
  * [9] L. Merten, P. Da Vela, A. Reimer, M. Boughelilba, J.P. Lundquist, S. Vorobiov, J. Becker Tjus, _Efficient Modeling of Heavy Cosmic Rays Propagation in Evolving Astrophysical Environments_ _PoS_ (2023)