††institutetext: Department of Physics and Technology, University of Bergen,
5007 Bergen, Norway

# A unified picture of medium-induced radiation

Johannes Hamre Isaksen  Adam Takacs  and Konrad Tywoniuk
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We revisit the picture of jets propagating in the quark-gluon plasma. In
addition to vacuum radiation, partons scatter on the medium constituents
resulting in induced emissions. Analytical approaches to including these
interactions have traditionally dealt separately with multiple, soft, or rare,
hard scatterings. A full description has so far only been available using
numerical methods. We achieve full analytical control of the relevant scales
and map out the dominant physical processes in the full phase space. To this
aim, we extend existing expansion schemes for the medium-induced spectrum to
the Bethe–Heitler regime. This covers the whole phase space from early to late
times, and from hard splittings to emissions below the thermal scale. Based on
the separation of scales, a space-time picture naturally emerges: at early
times, induced emissions start to build from rare scatterings with the medium.
At a later stage, induced emissions due to multiple soft scatterings result in
a turbulent cascade that rapidly degrades energy down to, and including, the
Bethe–Heitler regime. We quantify the impact of such an improved picture,
compared to the current state-of-the-art factorization that includes only soft
scatterings, by both analytical and numerical methods for the medium-induced
energy distribution function. Our work serves to improve our understanding of
jet quenching from small to large systems and for future upgrades of Monte
Carlo generators.

##  1 Introduction

Short-lived droplets of hot and dense nuclear matter, called the quark-gluon
plasma (QGP), are produced in relativistic heavy-ion collisions at RHIC and
LHC. Embedded in the same high-energy collisions, hard QCD processes are also
present, resulting in the production of collimated sprays of energetic
particles that are commonly referred to as jets Ellis:318585 ;
Dokshitzer:1991wu . Jets are well-understood, perturbative objects within
perturbative QCD and they are described up to high precision in proton-proton
collisions Larkoski:2017jix ; Marzani:2019hun ; Dasgupta:2020fwr . During
their propagation, however, jet particles can interact with the surrounding
nuclear matter. The modification of jet features, therefore, reflects the
properties of the QGP created in heavy-ion collisions Peigne:2008wu ; Mehtar-
Tani:2013pia ; Ghiglieri:2015zma ; Andrews:2018jcm . Currently, a vigorous
experimental program dedicated to quantifying jet modifications is ongoing at
both RHIC and LHC, focusing on a broad set of observables which includes
measurements of the modification of the jet spectrum, jet substructure
observables, and jet correlations Armesto:2015ioy ; Connors:2017ptx ;
Kogler:2018hem ; Cunqueiro:2021wls  (for a selection of predictions see Refs.
Chien:2015hda ; Mehtar-Tani:2016aco ; KunnawalkamElayavalli:2017hxo ;
Caucal:2019uvr ; Casalderrey-Solana:2019ubu ; Ringer:2019rfk ; Caucal:2020xad
; Caucal:2021cfb ; Mehtar-Tani:2021fud ; Takacs:2021bpv ; Attems:2022ubu ).

High-energy jets are particularly suitable probes of the QGP because their
energy scale $Q_{\rm jet}$ is much larger than the typical momentum scale of
the medium $Q_{\rm med}$. If this is the case, the impact of medium
modifications should therefore not affect the internal structure of the jet,
which would still rely on perturbative QCD Casalderrey-Solana:2012evi ;
Caucal:2018dla ; Vaidya:2020lih . A key ingredient when considering jet
modifications is the radiation induced by scatterings with the deconfined
medium constituents. Such emissions typically appear at scales comparable to
$Q_{\text{med}}$. Emission at scales much higher than $Q_{\rm med}$, on the
other hand, are unaffected by the medium, resulting in a factorized picture
between vacuum and medium processes Mehtar-Tani:2017web ; Caucal:2018dla .
Medium-induced emissions redistribute the original jet parton energy to
multiple, soft particles over large angles, including out of the jet cone.
This leads to a net jet energy loss which, in turn, is manifested as a
suppression of the jet spectrum (for an updated discussion of jet quenching
see Refs. Mehtar-Tani:2021fud ; Takacs:2021bpv , and for applications to
substructure see Refs. Mehtar-Tani:2016aco ; Caucal:2019uvr ; Caucal:2021cfb
). Consequently, medium-induced emissions are a crucial component of jet
energy loss and thus of phenomenological studies of jet observables in heavy-
ion collisions.

The medium-induced emission spectrum was formulated a long time ago
Baier:1996kr ; Zakharov:1997uu ; Arnold:2002ja . Previous solutions were
limited to either (i) expanding in the number of scatterings (referred to as
the opacity expansion) Gyulassy:1999zd ; Wiedemann:2000za , or (ii)
considering multiple soft scatterings (called the harmonic oscillator
approximation) Baier:1996kr ; Zakharov:1997uu ; Salgado:2003gb ;
Armesto:2003jh . Meanwhile, several works focused on the underlying scales
that separate the limiting cases Arnold:2008iy ; Arnold:2009mr ;
Blaizot:2012fh ; Kurkela:2014tla ; Dominguez:2019ges . The full problem has
also been tackled by numerical techniques Zakharov:2004vm ; CaronHuot:2010bp
(or more recently in Refs. Feal:2018sml ; Ke:2018jem ; Andres:2020vxs ;
Andres:2020kfg ; Schlichting:2020lef ). Not long ago, analytical techniques
were developed that provided a unified description of the multiple, soft and
rare, hard scatterings in a dense medium Mehtar-Tani:2019tvy ; Mehtar-
Tani:2019ygg ; Barata:2020sav ; Barata:2020rdn ; Barata:2021wuf , which better
match the full numerical solutions. The main challenge, common to both the
numeric and analytic approaches, resides in dealing with multiple interactions
with the underlying medium.

In this paper, we revisit the different analytic approaches to resumming
multiple interactions for calculating the medium-induced emission spectrum.
These include the opacity expansion (OE) and the improved opacity expansion
(IOE), which includes harmonic oscillator approximation as the leading term.
Moreover, we rigorously derive the resummed opacity expansion (ROE) for the
first time, which extends the description of the spectrum to low energy
emissions in the Bethe–Heitler regime. We provide a novel unified picture of
these resummation schemes by identifying their relevant emergent scales and
demonstrating their respective regions of validity. For example, we show that
the single scattering approximation, contained in the leading order of OE, is
valid even for a big medium, where one would expect more than one scattering
if the emitted energy is high enough. We show for the first time that the full
phase space of medium-induced emissions, spanning from the maximal jet energy
to the thermal scale, is covered by a union of these expansions. Each of the
expansions is also associated with the corresponding physical scattering
processes, and thus we reinterpret the frequently used terms such as GLV
emissions, coherent scatterings, and Bethe–Heitler region in a unified
framework.

As a next step, we identify the regimes where not only multiple interactions
are important, but also multiple emissions. These conditions are met for
sufficiently soft emissions in a large medium. The previously established
hierarchy of emergent scales plays a crucial role in mapping out early, rare,
and relatively hard emissions and a successive cascade of soft splittings. In
this context, hard medium-induced splittings can be thought of as extra
sources, in addition to the parent parton, for the full cascade. This
description is realized analytically in a novel scheme that combines a fixed
order expansion of rare emissions with an all-order resummation of soft
splittings. The resulting energy distribution links the asymptotic early and
late time behaviors for which analytical solutions exist. Finally, we resum
multiple induced emissions numerically, using the previously obtained precise
determination of the in-medium splitting rates, to calculate the energy
distribution function. We highlight the interplay of rare hard scatterings,
coherent soft splittings, and Bethe–Heitler emissions in a finite medium,
providing a state-of-the-art resummation.

Our reorganized picture helps not only with the physical understanding of
induced emissions, but provides a fast and efficient way to calculate the
medium-induced spectrum, which is a key ingredient for estimating jet energy
loss. It also serves to inform Monte Carlo algorithms simulating full jet
evolution inside the medium about how to implement multiple medium-induced
emissions and how to combine them with vacuum-like emissions, e.g. see in Ref.
Caucal:2018dla .

The paper is structured as follows. In Sec. 2, as an introduction, we discuss
the structure of the induced spectrum in the various regimes using heuristic
arguments, and we show how the radiation in the different regimes is related
to single soft, multiple soft and single hard scatterings with the medium, see
Fig. 3. The spectrum is calculated in detail in Sec. 3. We revisit the opacity
expansion and the improved opacity expansion schemes, and put on a firm
footing a novel resummation scheme, dubbed resummed opacity expansion, which
is valid for emissions below the Bethe–Heitler scale. Improving on previous
discussions, we provide formulas for the spectrum at arbitrary order and
calculate it exactly up to second order in all the expansions, allowing us for
the first time to establish regions where they converge. Finally, in Sec. 4 we
consider the problem of multiple emissions. We analyze induced particles
coming from the full phase space and confirm the importance of considering
multiple emissions, especially in the soft sector. In order to facilitate an
improved analytical understanding of the problem, we finally suggest a
resummation scheme of multiple emissions by iterating in rare, hard emissions
and including an arbitrary number of soft splittings. This is compared to the
full numerical results. We conclude with an outlook in Sec. 6. The appendix
contains lots of useful formulas, including the rate of emissions and
finite-$z$ corrections that are important for phenomenology. The code we have
developed to calculate the kernels and solve for the energy distribution is
provided at https://github.com/adam-takacs/kernels.git.

Figure 1: Three regimes of the radiative spectrum in dense media,
$L\gg\lambda$: the Bethe–Heitler regime $\omega<\omega_{\text{\tiny BH}}$
(leftmost), the BDMPS-Z regime $\omega_{\text{\tiny BH}}<\omega<\omega_{c}$
(middle), and the hard GLV regime $\omega_{c}<\omega$ (rightmost figure). The
size of the blob represents the typical formation time of the emission, which
is ${t_{\text{f}}}<\lambda$ in the leftmost panel and ${t_{\text{f}}}>\lambda$
in the center and rightmost panels. Black, wavy lines represent soft
interactions with the medium, where the momentum transfer is of the order of
the medium scale $|{\bm{q}}|\lesssim\mu$, while the red, wavy line represents
a hard scattering event, with $|{\bm{q}}|\gg\mu$.

##  2 Heuristic discussion of the medium-induced spectrum

The spectrum of gluon emissions induced by scatterings on a deconfined medium
plays a central role in the phenomenology of jet quenching. However, a full
understanding of all regimes have so far been lacking analytically and was
previously only achievable through numerical methods. Here, we present a
unified view of all relevant medium scales and their related regimes.

We formulate the spectrum of induced emissions by focusing on the relevant
length scales:

  * •

The size of the medium $L$ (or the length of propagation $t<L$).

  * •

The mean free path of the medium $\lambda\sim\frac{1}{n\sigma_{\rm tot}}$,
which combines the density $n$ and the scattering strength $\sigma_{\rm
tot}\sim\int{\rm d}\sigma$, and it describes the distance between typical
scatterings.

  * •

The formation time of an emission $t_{\rm f}=\frac{2\omega}{\bm{k}^{2}}$,
where $\omega$ is the energy and $\bm{k}$ is the transverse momentum of the
emission.

In addition, the relation between the in-medium screening scale $\mu$ and the
range of available transverse momenta $|{\bm{k}}|$ is also important.

The opacity $\chi\equiv L/\lambda$ characterizes the denseness of the medium.
If the opacity is small ($L\ll\lambda$), the medium is “dilute”, or weakly
interacting, while it is “dense”, or strongly interacting if
$L\gg\lambda$.111From this perspective, we have “fixed” $L$ and vary
$\lambda$. Naturally, we could also have identified these two regimes as a
“large” and “small” media, where we “fix” the mean free path and vary $L$
instead. The dilute medium barely consists of scattering centers, however, in
the dense medium one should account for an arbitrary number of interactions.

A big part of this paper will be about calculating the spectrum of medium-
induced emissions. For reference, the vacuum spectrum reads,

$\omega\frac{{\rm d}I}{{\rm d}\omega}\sim\alpha_{s}\int\frac{{\rm
d}{\bm{k}}^{2}}{{\bm{k}}^{2}}\,,$ (1)

where currently we did not specify the limits of the transverse momentum
integral. This contains the well-known soft ($\omega\to 0$) and collinear
(${\bm{k}}^{2}\to 0$) divergences. In contrast, the collinear divergence in
the medium spectrum is removed by the need to exchange transverse momentum
with the medium.

Motivated by this, we introduce a heuristic model that captures the main
features of the medium-induced spectrum, given by

$\omega\frac{{\rm d}I}{{\rm d}\omega}\sim\alpha_{s}L\int{\rm
d}^{2}\bm{k}\,\sigma(\bm{k})\sim\alpha_{s}^{3}nL\int\frac{{\rm
d}\bm{k}^{2}}{(\bm{k}^{2}+\mu^{2})^{2}}\,,$ (2)

where $n$ is the medium density and $\mu$ a screening mass. The behavior at
high-${\bm{k}}$, $\sigma\sim\bm{k}^{-4}$, reproduces the expected Coulomb
tail. The factor $L$ arises since the emission can take place anywhere along
the medium length. For a more accurate description of medium-induced
emissions, we refer the reader to Sec. 3.

In our effective description, we focus on the hierarchy among the introduced
scales and show the separation of different scattering regions. Firstly, in
the $t_{\rm f}\gg L$ limit, the formation of the emission extends beyond the
medium, where one naturally should expect vacuum physics to dominate.222For
such soft emissions, medium effects can influence the color coherence
properties leading to a modification of the phase space Mehtar-Tani:2010ebp ;
Mehtar-Tani:2012mfa . We will hence not consider this possibility here. The
remaining cases are listed below:

####  Dilute media ($t_{\rm f}\leq L\ll\lambda$):

In case of a low medium opacity, we expect that roughly one scattering occurs.
This process will typically transfer a momentum of order of the Debye mass to
the emitted gluon, or $\langle{\bm{k}}^{2}\rangle\sim\mu^{2}$, leading to
${t_{\text{f}}}=2\omega/\mu^{2}$. The formation of the gluon has to take place
inside the medium, giving rise to the characteristic energy scale in the
dilute regime, namely

$\bar{\omega}_{c}=\frac{1}{2}\mu^{2}L\,.$ (3)

This separates two regimes of emissions that are sourced via different
scattering processes: on the one side soft gluons with
$\omega<\bar{\omega}_{c}$, generated via a soft scattering with the medium
$\langle{\bm{k}}^{2}\rangle\lesssim\mu^{2}$. Hard gluons with
$\omega>\bar{\omega}_{c}$ can also be generated, however those demand a large
momentum exchange with the medium, $\langle{\bm{k}}^{2}\rangle>\mu^{2}$, which
is comparatively rare. Let us now consider how the spectrum behaves in these
two distinct regimes.

According to our discussion above, the soft production should be dominated by
soft transverse momentum exchanges with the medium. Hence, we expect that the
spectrum of emitted gluons goes as

$\left.\omega\frac{{\rm d}I}{{\rm
d}\omega}\right|_{\omega<\bar{\omega}_{c}}\sim\alpha_{s}L\,\int_{0}^{\infty}{\rm
d}\bm{k}^{2}\frac{\alpha_{s}^{2}n}{(\bm{k}^{2}+\mu^{2})^{2}}\sim\alpha_{s}\frac{L}{\lambda}\,,$
(4)

where the integral is dominated by ${\bm{k}}^{2}\lesssim\mu^{2}$. This
integral gives the proportionality with the inverse mean free path i.e.
$n/\mu^{2}\sim 1/\lambda$, resulting in an overall factor of medium opacity
$L/\lambda$. This parametric estimate misses an important logarithmic factor
$\sim\ln\bar{\omega}_{c}/\omega$, see a further discussion in Sec. 3.2, which
signals that the simplifications pertaining to the “soft” regime break down at
$\omega\approx\bar{\omega}_{c}$.

For hard emissions, $\omega>\bar{\omega}_{c}$, we instead get that

$\left.\omega\frac{{\rm d}I}{{\rm
d}\omega}\right|_{\omega>\bar{\omega}_{c}}\sim\alpha_{s}^{3}nL\,\int_{\omega/L}^{\infty}\frac{{\rm
d}{\bm{k}}^{2}}{{\bm{k}}^{4}}\sim\alpha_{s}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\,,$
(5)

where we used ${t_{\text{f}}}=2\omega/\bm{k}^{2}<L$, and neglected the
screening mass $\mu^{2}$ in this parametric regime, since
$\langle\bm{k}^{2}\rangle\gg\mu^{2}$. Compared to the soft regime from Eq. (4)
it is suppressed by an additional power of $\bar{\omega}_{c}/\omega\ll 1$. The
complete spectrum in the dilute regime is sketched in Fig. 2 (left).

####  Dense media with long formation time ($\lambda\ll t_{\rm f}\ll L$):

In a dense medium we should expect that typically many scatterings occur
during the emission process, which is illustrated in the middle of Fig. 1.
This allows to approximate the total transferred transverse momentum by
$\langle{\bm{k}}^{2}\rangle\sim\hat{q}t$, which resembles a random walk for
$t$ time in two dimensions, with $\hat{q}$ playing the role of a diffusion
constant. This constant determines the typical transverse momentum accumulated
per unit length, or $\hat{q}\sim n\sim\mu^{2}/\lambda$. In this case the
formation time becomes

${t_{\text{f}}}=\sqrt{\frac{2\omega}{\hat{q}}}\,.$ (6)

This is often called the coherence length, since during the formation time,
interference effects between multiple scattering with the medium are active
and the gluon feels only one effective scattering center. The accumulated
transverse momentum during the splitting process is in this case
$\langle{\bm{k}}^{2}\rangle=\sqrt{2\omega\hat{q}}$, which is the celebrated
Landau-Pomeranchuk-Migdal (LPM) effect.

Again, comparing the formation time to the medium length, leads to the
characteristic energy scale in the dense regime, namely

$\omega_{c}=\frac{1}{2}\hat{q}L^{2}\,,$ (7)

and thus $\omega<\omega_{c}$ for multiple soft scatterings. The maximal
possible momentum accumulated via multiple soft scatterings is denoted
$\langle{\bm{k}}^{2}\rangle\sim Q_{s}^{2}=\hat{q}L$. The other limiting scale
of the multiple scattering regime arises when considering the minimal
formation time in this hierarchy, i.e. ${t_{\text{f}}}>\lambda$, giving rise
the scale $\omega>\omega_{\text{\tiny BH}}$, (see later in Eq. (10)). In this
case, the accumulated transverse momentum squared reduces to a single soft
scattering $\langle{\bm{k}}^{2}\rangle\sim\hat{q}\lambda\sim\mu^{2}$.

In the multiple soft scattering regime, characterized by $\omega_{\text{\tiny
BH}}\ll\omega\ll\omega_{c}$, the mean free path has to be replaced by the
formation time in Eq. (4), leading to

$\left.\omega\frac{{\rm d}I}{{\rm d}\omega}\right|_{\omega_{\text{\tiny
BH}}\ll\omega\ll\omega_{c}}\sim\alpha_{s}\frac{L}{{t_{\text{f}}}}\sim\alpha_{s}\sqrt{\frac{\hat{q}L^{2}}{\omega}}\,.$
(8)

This is also often referred to as the BDMPS-Z spectrum in the soft limit.

For hard gluon emissions, $\omega\gg\omega_{c}$, we also have to demand that
$\langle{\bm{k}}^{2}\rangle\gg Q_{s}$. In other words, only a hard scattering
can provide sufficient transverse momentum to fulfill all the conditions. The
relevant contribution is therefore captured by Eq. (5) and, remarkably, the
spectrum in this limit is identical to the hard tail in the dilute regime,
namely

$\left.\omega\frac{{\rm d}I}{{\rm
d}\omega}\right|_{\omega\gg\omega_{c}}\sim\alpha_{s}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\,.$
(9)

This demonstrates that, even in a dense medium, hard emissions mostly are
driven by single, rare hard scattering events. An illustration of this can be
seen on the right in Fig. 1.

Figure 2: A sketch of the spectrum of medium-induced gluons in a dilute medium
$L\ll\lambda$ (left) and in a dense medium $L\gg\lambda$ (right).

####  Dense media with short formation time ($t_{\rm f}\ll\lambda\ll L$):

The picture we just described should hold as long as there indeed is time for
multiple scatterings during the emission process, namely that
${t_{\text{f}}}>\lambda$. However, when the formation time is short the parton
will only have time to scatter once before it splits. This is illustrated on
the left in Fig. 1. The transverse scale is typically soft
$\langle\bm{k}^{2}\rangle\sim\mu^{2}$, and thus
${t_{\text{f}}}=\frac{2\omega}{\mu^{2}}$. This regime is characterized by
${t_{\text{f}}}\ll\lambda$, or equivalently as a condition on the energy
$\omega\ll\omega_{\text{\tiny BH}}$, where we have defined the scale

$\omega_{\text{\tiny BH}}=\frac{1}{2}\mu^{2}\lambda\,.$ (10)

Note that $\bar{\omega}_{c}(\lambda)\equiv\omega_{\text{\tiny BH}}$. In this
case the spectrum becomes

$\left.\omega\frac{{\rm d}I}{{\rm d}\omega}\right|_{\omega<\omega_{\text{\tiny
BH}}}\sim\alpha_{s}\frac{L}{\lambda}\,,$ (11)

which is similar to the result in Eq. (4) and it is some times referred as
Bethe–Heitler region because the QED analogue. Our heuristic analysis fails to
capture the additional logarithmic term which in this case comes with
$\ln\omega_{\text{\tiny BH}}/\omega$, see Sec. 3.3. A sketch of the spectrum
in a dense medium can be found in Fig. 2 (right).

#### Summary:

Bringing together the heuristic arguments of this section, we show a sketch of
the emission spectrum for dilute and dense media in Fig. 2. The full emission
phase space is divided by three lines corresponding to the emergent scales:
$\bar{\omega}_{c}(t)=\frac{1}{2}\mu^{2}t$ in the dilute regime ($L<\lambda$),
and $\omega_{\text{\tiny BH}}=\frac{1}{2}\mu^{2}\lambda$ and
$\omega_{c}(t)=\frac{1}{2}\hat{q}t^{2}$ in the dense regime ($L>\lambda$) and
they are shown in Fig. 3. The $\omega_{c}(t)$ line is not completely straight
because $\hat{q}$ in general is $\omega$ dependent (see in Sec. 3.4). When the
medium size is of the order of the mean free path they all collapse to the
same value, i.e.
$\bar{\omega}_{c}(\lambda)=\omega_{c}(\lambda)=\omega_{\text{\tiny BH}}$.
Typically, we adopt a notation where the scales written without the
$t$-argument denote their respective values at $L$, e.g.
$\omega_{c}\equiv\omega_{c}(L)$.

These scales delineate three distinct regimes of scattering processes and thus
induced emissions that were discussed in the preceding paragraphs. The areas
between these scales are governed by few soft, multiple soft and rare hard
interactions with the medium, as depicted with colors in Fig. 3 and discussed
above. We also show two length scales: the mean free path $\lambda$ and the
critical medium length $L_{c}$. The mean free path marks the time where
multiple scatterings appear. The critical medium length indicates where rare,
hard scatterings will no longer have an effect, that is where
$\omega_{c}(t)=E$, leading to $L_{c}=\sqrt{2E/\hat{q}}$.333When considering a
finite splitting fraction $z$ the exact definition turns out to be
$L_{c}=\sqrt{E/(2\hat{q})}$.

Figure 3: The phase space for medium-induced emissions designated to the
leading scattering processes. The main length scales of the problem are
$\lambda$ and the critical length $L_{c}$, corresponding to the characteristic
energies $\omega_{\text{\tiny BH}}$ and $E$, see text for further details.

As the parton moves through the medium, at each instant $t$ the emission rate
$\omega{\rm d}I/({\rm d}\omega{\rm d}t)$ is different since the phase space
for available emissions changes. The emission spectrum $\omega{\rm d}I/{\rm
d}\omega$, evaluated at a given $t$ (typically $t=L$), includes the
accumulated range of processes that occurred up to that time. The main goal of
this work is to consider multiple medium-induced emissions in this system of
dynamically evolving scales presented in Fig. 3. In the next section, Sec. 3,
we will formally derive the results we only have argued for in this section.
Finally, in Sec. 4 we will tackle the issue of multiple emissions in this
scheme.

##  3 Spectrum of medium-induced emissions

In this section, we derive the all order emission spectrum induced by elastic
scatterings on a deconfined medium. Of the three expansions we present here
the opacity expansion and the improved opacity expansion have been discussed
in depth in previous works, see Refs. Gyulassy:1999zd ; Wiedemann:2000za  and
Mehtar-Tani:2019tvy ; Mehtar-Tani:2019ygg ; Barata:2020sav ; Barata:2020rdn ;
Barata:2021wuf  respectively. The resummed opacity expansion has been argued
for before, see Wiedemann:2000za ; Andres:2020kfg ; Schlichting:2021idr , but
is here derived rigorously for the first time. This paper strives to be a
comprehensive reference for all of the expansions, and hence they are all
presented in detail. Furthermore, we extend previous calculations to all
orders and present results to order $N=2$ for the opacity expansion and
$N_{r}=2$ for the resummed opacity expansion, and extract the relevant limits
analytically. This provides valuable insight into the underlying structure of
the expansions in different regimes.

###  3.1 General formalism

Currently, we consider the emission of a gluon with energy $\omega$ from a
parent parton with energy $E$ in the soft limit, i.e. $\omega\ll E$. The soft
limit is used in this section because it is much more clear and readable. For
a description beyond the strictly soft limit we refer to App. C, which
includes novel results.

Our starting point is the definition of the spectrum of medium-induced gluons
Baier:1996kr ; Zakharov:1997uu ; Wiedemann:2000za ; Arnold:2002ja ,

$\omega\frac{{\rm d}I}{{\rm d}\omega}=\frac{2\alpha_{s}C_{R}}{\omega^{2}}{\rm
Re}\int_{0}^{\infty}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm
d}t_{1}\,\partial_{\bm{x}}\cdot\partial_{\bm{y}}\big{[}{\cal
K}({\bm{x}},t_{2};{\bm{y}},t_{1})-{\cal
K}_{0}({\bm{x}},t_{2};{\bm{y}},t_{1})\big{]}_{{\bm{x}}={\bm{y}}=0}\,,$ (12)

where $C_{R}$ is the Casimir color factor of the emitting particle
($C_{R}=C_{F}$ for a quark and $C_{R}=N_{c}$ for a gluon).444This expression
can be derived directly from the fully $z$ dependent spectrum in Eq. (138),
see the discussion in App. C. The three-point correlator ${\cal K}$ solves the
Schrödinger equation

$\left[i\partial_{t}+\frac{\partial_{\bm{x}}^{2}}{2\omega}+iv({\bm{x}},t)\right]{\cal
K}({\bm{x}},t;{\bm{y}},t_{0})=i\delta(t-t_{0})\delta({\bm{x}}-{\bm{y}})\,,$
(13)

where the potential $v({\bm{x}},t)$ describes scatterings in a thermal or
quasi-particle like background,

$v({\bm{x}},t)=\int_{\bm{q}}\sigma({\bm{q}},t)\left(1-{\rm
e}^{i{\bm{q}}\cdot{\bm{x}}}\right)\,.$ (14)

Here $\sigma({\bm{q}},t)=N_{c}{\rm d}^{2}\sigma_{\rm el}/{\rm d}^{2}{\bm{q}}$
is the medium elastic scattering potential.555Throughout, we adopt a shorthand
notation, so that $\int_{p}=\int\frac{{\rm d}^{4}p}{(2\pi)^{4}}$,
$\int_{{\bm{p}}}=\int\frac{{\rm d}^{2}{\bm{p}}}{(2\pi)^{2}}$, and
$\int_{\bm{x}}=\int{\rm d}^{2}{\bm{x}}$. The color factor $N_{c}$ appears
because it is only the emitted gluon that is allowed to pick up transverse
momentum in the medium. This potential can be extracted from an effective
theory that accounts for both large and small momentum exchanges with the
medium Moore:2021jwe . In the main part of this paper, we will use the
Gyulassy-Wang potential Gyulassy:1993hr  that both contains a hard Coulomb-
tail and implements screening in the infrared $\mu^{2}\gg{\bm{q}}^{2}$,

$\frac{{\rm d}^{2}\sigma_{\rm el}}{{\rm
d}^{2}{\bm{q}}}=\frac{g^{4}n(t)}{({\bm{q}}^{2}+\mu^{2})^{2}}\,,$ (15)

where $n$ is the density of medium scattering centers and $\mu$ is a screening
mass of the order of the Debye mass of a thermal medium. We also provide the
spectrum in App. E with the LO hard thermal loop (HTL) potential
Aurenche:2002pd . Also, in Eq. (13) we neglect quark and gluon thermal masses,
which corresponds to taking the high-energy limit ($E\gg
m_{q,g}^{4}/\hat{q}\sim\omega_{\text{\tiny BH}}$), see e.g. Refs.
Arnold:2008iy ; Schlichting:2021idr  for further discussion.

When the medium is not present, $v=0$, one recovers the propagation of a
single parton in vacuum, ${\cal K}({\bm{x}},t_{2};{\bm{y}},t_{1})\equiv{\cal
K}_{0}({\bm{x}}-{\bm{y}},t_{2}-t_{1})$ with ${\cal
K}_{0}({\bm{x}},t)=\frac{\omega}{2\pi i\,t}\exp[i\,\omega{\bm{x}}^{2}/(2t)]$.
To only capture medium effects, the vacuum term is explicitly subtracted in
Eq. (12).

On many occasions it is more practical to work in transverse momentum space,

${\cal K}({\bm{p}},t_{2};{\bm{p}}_{0},t_{1})=\int_{{\bm{x}},{\bm{y}}}\,{\rm
e}^{-i{\bm{p}}\cdot{\bm{x}}+i{\bm{p}}_{0}\cdot{\bm{y}}}{\cal
K}({\bm{x}},t_{2};{\bm{y}},t_{1})\,.$ (16)

The vacuum propagator ${\cal K}_{0}$ then becomes a plane wave, i.e. ${\cal
K}_{0}({\bm{p}},t)=\exp[-i\,{\bm{p}}^{2}t/(2\omega)]$. In this representation,
the solution to the Schrödinger equation (13) can be written as the recursive
equation

$\displaystyle{\cal K}({\bm{p}},t_{2};{\bm{p}}_{0},t_{1})$
$\displaystyle=(2\pi)^{2}\delta({\bm{p}}-{\bm{p}}_{0}){\cal
K}_{0}({\bm{p}};t_{2}-t_{1})$ $\displaystyle-\int_{t_{1}}^{t_{2}}{\rm
d}s\int_{\bm{q}}\,{\cal K}_{0}({\bm{p}};t_{2}-s)v({\bm{q}},s){\cal
K}({\bm{p}}-{\bm{q}},s;{\bm{p}}_{0},t_{1})\,,$ (17)

where now

$v({\bm{q}},s)=(2\pi)^{2}\delta({\bm{q}})\Sigma(s)-\sigma({\bm{q}},s)\,,$ (18)

ensures probability conservation. Here,
$\Sigma(s)\equiv\int_{\bm{q}}\,\sigma({\bm{q}},s)$ can be interpreted as the
inverse of the (local) mean free path $\lambda$ along a trajectory of a
propagating parton, or

$\lambda(s)=\frac{1}{\Sigma(s)}\,.$ (19)

In these expressions, we have assumed that the integral over the elastic
scattering cross section exists. In many cases, e.g. for the HTL potential
Aurenche:2002pd , one needs to introduce an IR regulator. However,
$v({\bm{q}},s)$ in Eq. (18) is not sensitive to this IR regulation and
therefore the expansion in Eq. (3.1) is well-defined. We have provided a
further discussion of the HTL potential in App. E.

The medium-induced spectrum now reads

$\omega\frac{{\rm d}I}{{\rm d}\omega}=\frac{2\alpha_{s}C_{R}}{\omega^{2}}{\rm
Re}\int_{0}^{\infty}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm
d}t_{1}\int_{{\bm{p}},{\bm{q}}}\,{\bm{p}}\cdot{\bm{q}}\,\big{[}{{\cal
K}}({\bm{p}},t_{2};{\bm{q}},t_{1})-(2\pi)^{2}\delta({\bm{p}}-{\bm{q}}){\cal
K}_{0}({\bm{p}},t_{2}-t_{1})\big{]}\,,$ (20)

The vacuum contribution can then be removed by inserting Eq. (3.1) into Eq.
(20), yielding

$\omega\frac{{\rm d}I}{{\rm d}\omega}=\frac{4\alpha_{s}C_{R}}{\omega}{\rm
Re}\,i\int_{0}^{L}{\rm d}t_{2}\int_{0}^{t_{2}}{\rm
d}t_{1}\int_{{\bm{p}},{\bm{p}}_{0},{\bm{q}}}\,\frac{{\bm{p}}\cdot{\bm{p}}_{0}}{{\bm{p}}^{2}}v({\bm{q}},t_{2}){\cal
K}({\bm{p}}-{\bm{q}},t_{2};{\bm{p}}_{0},t_{1})\,,$ (21)

where we regulated the integral over the latter time coordinate using an
adiabatic turn-off at infinity (see also Ref. Andres:2020kfg ). The other time
integrals are limited by the extent of the medium $L$. Noticing, that

$\displaystyle\int_{{\bm{p}}}\,\frac{p^{i}}{{\bm{p}}^{2}}v({\bm{p}}-{\bm{k}},s)=\frac{k^{i}}{{\bm{k}}^{2}}\Sigma({\bm{k}}^{2},s)\,,$
(22)

where
$\Sigma({\bm{k}}^{2},s)=\int_{{\bm{q}}}\,\sigma({\bm{q}},s)\Theta({\bm{q}}^{2}-{\bm{k}}^{2})$,666Also
$\Sigma(0,s)=\Sigma(s)$, consistent with the definition in Eq. (18). we obtain

$\omega\frac{{\rm d}I}{{\rm d}\omega}=\frac{4\alpha_{s}C_{R}}{\omega}{\rm
Re}\,i\int_{0}^{L}{\rm d}t_{2}\int_{0}^{t_{2}}{\rm
d}t_{1}\int_{{\bm{p}},{\bm{p}}_{0}}\,\Sigma({\bm{p}}^{2},t_{2})\frac{{\bm{p}}\cdot{\bm{p}}_{0}}{{\bm{p}}^{2}}{\cal
K}({\bm{p}},t_{2};{\bm{p}}_{0},t_{1})\,.$ (23)

While the above results are valid for any medium potential, in this work we
will focus on the GW scattering potential, defined in Eq. (15). In this case,
we find that

$\Sigma({\bm{k}}^{2},s)=\frac{\hat{q}_{0}(s)}{{\bm{k}}^{2}+\mu^{2}}\,,$ (24)

where ${\hat{q}}_{0}(s)=4\pi\alpha_{s}^{2}N_{c}n(s)$ is a measure of the
scattering density. Currently, we consider a medium of constant density,
$n(s)=n_{0}$.

The spectrum, given by Eq. (12) or Eq. (23), can be evaluated using numerical
techniques Zakharov:2004vm ; CaronHuot:2010bp ; Feal:2018sml ; Ke:2018jem ;
Andres:2020vxs ; Andres:2020kfg ; Schlichting:2020lef  or by employing
analytical approximations. As we will see the approximate approaches rely on
expanding the problem as a series, which will give the true answer at infinite
order. The different series have different radiuses of convergence, and none
of them will alone converge for all $L$ and $\omega$, meaning more than one
have to be employed. In most cases, however, the first order expansion is
sufficient to provide an accurate approximation of the all order result. In
the following we discuss three well-defined approaches that together provide
an accurate description of the true problem for all $L$ and $\omega$, called:
the opacity expansion, the resummed opacity expansion, and the improved
opacity expansion. We will derive these, discuss their limits and their
regions of validity.

###  3.2 Opacity expansion (OE)

The opacity expansion of the spectrum arises when inserting Eq. (3.1) directly
into Eq. (23), and was developed in Refs. Gyulassy:1999zd ; Wiedemann:2000za .
The truncation of this series at a given order $n$ in the medium scattering
potential gives the $N=n$ term, which is by definition proportional to
$(L/\lambda)^{n}$ (see Eq. (A.1)). Physically this means, at $N=n$ one counts
$n$ number of scatterings (both with and without momentum exchange) on the
full elastic potential. The relevant energy scales that arise are
$\bar{\omega}_{c}=\frac{1}{2}\mu^{2}L$, and
$\frac{L}{\lambda}\bar{\omega}_{c}=\frac{1}{2}\hat{q}_{0}L^{2}$ as discussed
in Sec. 2. A general formula for the spectrum at any order is derived in App.
A, and with finite-$z$ corrections in App. C. These results are used in the
following calculations, and we will refer to the appendices for more details.

#### First order ($N=1$):

The spectrum at first order of opacity is well known Gyulassy:1999zd ;
Wiedemann:2000za . Since Eq. (23) already includes at least one scattering, we
obtain the $N=1$ term by replacing the full propagator ${\cal K}$ by the
vacuum one. We then find,

$\displaystyle\omega\frac{{\rm d}I^{N=1}}{{\rm
d}\omega}=8\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\int_{\bm{p}}\,\tilde{\Sigma}({\bm{p}}^{2})\,{\rm
Re}\,i\int_{0}^{1}{\rm d}t_{1}\int_{0}^{t_{1}}{\rm d}t_{0}\,{\rm
e}^{-i{\bm{p}}^{2}(t_{1}-t_{0})}\,,$ (25)

where we have switched to dimensionless integration variables by defining
${\bm{p}}^{2}\to{\bm{p}}^{2}L/(2\omega)$ and $t\to t/L$, and where
$\tilde{\Sigma}({\bm{p}}^{2})=({\bm{p}}^{2}+\bar{\omega}_{c}/\omega)^{-1}$.
This expression can also be obtained from the general $N=n$ result in Eq.
(A.1). After simplifications, the spectrum becomes

$\displaystyle\omega\frac{{\rm d}I^{N=1}}{{\rm
d}\omega}=2{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\int_{0}^{\infty}{\rm
d}p\,\frac{1}{p+\frac{\bar{\omega}_{c}}{\omega}}\frac{p-\sin{p}}{p^{2}}\,,$
(26)

where ${\rm\bar{\alpha}}=\alpha_{s}C_{R}/\pi$ and
$\lambda=\mu^{2}/{\hat{q}}_{0}$. We recognize the dependence on the medium
opacity $L/\lambda$ and the ratio $\bar{\omega}_{c}/\omega$. The remaining
integral can be done analytically, but the resulting expression is not very
illuminating. However, the limiting behavior can readily be extracted,

$\omega\frac{{\rm d}I^{N=1}}{{\rm
d}\omega}\simeq\begin{cases}2{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\ln\frac{\bar{\omega}_{c}}{\omega}-1+\gamma_{E}\right),&\qquad\text{for
}\omega\ll\bar{\omega}_{c}\,,\\\\[0.2pt]
\frac{\pi}{2}{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega},&\qquad\text{for
}\omega\gg\bar{\omega}_{c}\,.\end{cases}$ (27)

This agrees well with the heuristic discussion in Sec. 2. In particular, we
identify a logarithmic behavior $\sim\ln\frac{\bar{\omega}_{c}}{\omega}$ in
the infrared. Notice the different expansion structures in the soft
$\sim\bar{\alpha}\frac{L}{\lambda}$ and in the hard
$\sim\bar{\alpha}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}$ limits,
which we will come back to.

#### Second order ($N=2$):

The calculation for $N=2$ follows in a similar way, leading to

$\displaystyle\omega\frac{{\rm d}I^{N=2}}{{\rm d}\omega}$
$\displaystyle=-8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\frac{\bar{\omega}_{c}}{\omega}\int_{{\bm{p}}_{2},{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{2}^{2})\frac{{\bm{p}}_{2}\cdot{\bm{p}}_{1}}{{\bm{p}}_{2}^{2}}\tilde{v}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times{\rm Re}\,i\int_{0}^{1}{\rm d}t_{2}\int_{0}^{t_{2}}{\rm
d}t_{1}\int_{0}^{t_{1}}{\rm d}t_{0}\,{\rm
e}^{-i{\bm{p}}_{2}^{2}(t_{2}-t_{1})}{\rm
e}^{-i{\bm{p}}_{1}^{2}(t_{1}-t_{0})}\,,$ (28)

with dimensionless integration variables, and where
$\tilde{v}({\bm{p}})=(2\pi)^{2}\delta({\bm{p}})-\frac{\bar{\omega}_{c}}{\omega}\tilde{\sigma}({\bm{p}})$.
In the GW model,
$\tilde{\sigma}({\bm{p}})=4\pi/({\bm{p}}^{2}+\frac{\bar{\omega}_{c}}{\omega})^{2}$.
After inserting $\tilde{v}$, doing the time integrals and simplifying this can
be written as

$\omega\frac{{\rm d}I^{N=2}}{{\rm
d}\omega}=-4\bar{\alpha}\left(\frac{L}{\lambda}\right)^{2}\frac{\bar{\omega}_{c}}{\omega}\left[{\cal
I}_{1}\left(\frac{\bar{\omega}_{c}}{\omega}\right)-\frac{\bar{\omega}_{c}}{\omega}{\cal
I}_{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)\right]\,,$ (29)

where we have defined the integrals

$\displaystyle{\cal I}_{1}\left(\frac{\bar{\omega}_{c}}{\omega}\right)$
$\displaystyle=\int_{0}^{\infty}{\rm
d}p\,\frac{1}{p+\frac{\bar{\omega}_{c}}{\omega}}\frac{1-\cos p-\frac{p}{2}\sin
p}{p^{3}}\,,$ (30) $\displaystyle{\cal
I}_{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)$
$\displaystyle=\int_{0}^{\infty}{\rm d}p_{2}\int_{0}^{\infty}{\rm
d}p_{1}\,\frac{p_{1}}{p_{2}+\frac{\bar{\omega}_{c}}{\omega}}\frac{1}{\left[\left(p_{2}+p_{1}+\frac{\bar{\omega}_{c}}{\omega}\right)^{2}-4p_{2}p_{1}\right]^{3/2}}$
$\displaystyle\times\frac{1}{p_{2}-p_{1}}\left[\frac{1}{p_{1}^{2}}\left(1-\cos{p_{1}}\right)-\frac{1}{p_{2}^{2}}\left(1-\cos{p_{2}}\right)\right]\,.$
(31)

The ${\cal I}_{1}$ integral can be done analytically, but ${\cal I}_{2}$ is
more complicated. It can be shown that it is much smaller than ${\cal I}_{1}$
in the soft limit. In the hard limit, ${\cal I}_{1}$ and ${\cal I}_{2}$ cancel
at the order of $\mathcal{O}(\frac{\bar{\omega}_{c}}{\omega})$, leaving a
positive contribution going as
$\mathcal{O}(\frac{\bar{\omega}_{c}}{\omega})^{2}$. In summary,

$\displaystyle\omega\frac{{\rm d}I^{N=2}}{{\rm
d}\omega}\simeq\begin{cases}-{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2},&\text{for
}\omega\ll\bar{\omega}_{c}\,,\\\\[0.5pt]
\sim{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{2},&\text{for
}\omega\gg\bar{\omega}_{c}\,.\end{cases}$ (32)

We notice that the $N=2$ is proportional to
$\bar{\alpha}(\frac{L}{\lambda})^{2}$ in the soft limit, and goes like
$\bar{\alpha}(\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega})^{2}$ in the
hard limit. This immediately implies that $N=2$ is always subleading to $N=1$
if the medium is dilute $L\ll\lambda$ or if the emission is hard
$\omega\gg\frac{L}{\lambda}\bar{\omega}_{c}$. Given the structure of the
expansion, we expect the previous statement to hold at arbitrary $N=n$ order.
This is in agreement with the earlier, heuristic observation in Ref.
Arnold:2008iy .

Based on the limits for $N=1$ and $N=2$, given by Eqs. (27) and (32), in the
regimes where the expansion holds the all order OE spectrum is expected to
take the form

$\omega\frac{{\rm d}I}{{\rm
d}\omega}=\begin{cases}{\bar{\alpha}\sum_{n=1}^{\infty}\left(\frac{L}{\lambda}\right)^{n}h_{n}\left(\frac{\omega}{\bar{\omega}_{c}}\right)}\,,&\omega\ll\bar{\omega}_{c}\,,\\\
{\bar{\alpha}\sum_{n=1}^{\infty}\left(\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\right)^{n}\tilde{h}_{n}\left(\frac{\bar{\omega}_{c}}{\omega}\right)}\,,&\omega\gg\bar{\omega}_{c}\,,\end{cases}$
(33)

where the OE coefficients $h_{n},\tilde{h}_{n}$ are finite and can be
calculated order by order. Note that we have not strictly proven this for all
orders, although our $N=1$ and $N=2$ results strongly indicate this structure.
In the soft limit $\omega\ll\bar{\omega}_{c}$, the OE expansion converges
rapidly, defining the expected “naive” radius of convergence $L/\lambda<1$.
However, in the hard limit $\omega\gg\bar{\omega}_{c}$ there is convergence
even if the medium is big, provided
$\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}<1$. The full region of
convergence is shown in green in the left panel of Fig. 4. Outside of this
region we expect higher orders to grow uncontrollably and hence the OE is not
valid when truncated at any finite order.

The resulting spectrum from Eqs. (26) and (29) is shown in the right panel of
Fig. 4 for different propagation lengths (labeled with $t$). For short lengths
$t<\lambda$ the OE is valid for all $\omega$. For $t>\lambda$, the OE is only
valid if $\frac{t}{\lambda}\frac{\bar{\omega}_{c}}{\omega}<1$ (see also the
green region on the left panel). We note that the $N=2$ correction becomes
important at $t>\lambda$ and $\omega\approx\frac{t}{\lambda}\bar{\omega}_{c}$
(the latter constraint shown as bullets in the figure). For larger media, the
grey bullets, representing the minimal energy for achieving convergence, moves
to higher values, and the truncated OE series at $\omega$ smaller than this
becomes ill-defined. This can be seen, for instance, in the upper line in Fig.
4 (right) for a medium length of $t=4$ fm. We have also compared to a full
numerical evaluation of the spectrum from Ref. Andres:2020kfg ; Andres:2022 .
The figure shows that this indeed is well approximated by the OE in its region
of validity, as we have argued.

Figure 4: Left: The sketch of the region of validity (and convergence) of the
opacity expansion for different propagation length $t$ and emission energy
$\omega$. Right: The induced emission spectrum for gluons in the opacity
expansion. The gray part of the curves denotes regions, where the expansion is
not valid. Using the parameters presented, $\lambda=0.06$ fm. The full numeric
solution is also presented with dashed lines.

###  3.3 Resummed opacity expansion (ROE)

Next, we turn to dense media, $L\gg\lambda$, where multiple scattering have to
be accounted for. However, as depicted to the left in Fig. 1 and discussed in
Sec. 2, for soft emissions with short formation times, a single scattering
still gives the leading contribution to the spectrum. This defines the so-
called Bethe–Heitler regime named after the QED analogue of this process. A
qualitative argument of this regime was first given in Ref. Wiedemann:2000za
and later developed in Ref. Andres:2020kfg , see also in Ref.
Schlichting:2021idr  which coined the name “resummed opacity expansion” (ROE).

Here, we provide for the first time a consistent framework for dealing with an
expansion of the real scatterings with the medium, whereby real we mean
interactions with a finite transverse momentum exchange. At the same time, an
all-order resummation of the corresponding virtual interactions, with zero
transverse momentum exchange, is performed.

Dividing Eq. (3.1) by the vacuum propagator leaves us with

$\displaystyle\frac{{\cal K}({\bm{p}},t;{\bm{p}}_{0},t_{0})}{{\cal
K}_{0}({\bm{p}};t-t_{0})}$
$\displaystyle=(2\pi)^{2}\delta({\bm{p}}-{\bm{p}}_{0})-\int_{t_{0}}^{t}{\rm
d}s\int_{\bm{q}}\,v({\bm{q}},s)\frac{{\cal
K}({\bm{p}}-{\bm{q}},s;{\bm{p}}_{0},t_{0})}{{\cal
K}_{0}({\bm{p}};s-t_{0})}\,.$ (34)

Next, taking a derivative with respect to the latest time results in

$\frac{\partial}{\partial t}\frac{{\cal
K}({\bm{p}},t;{\bm{p}}_{0},t_{0})}{{\cal
K}_{0}({\bm{p}};t-t_{0})}+\Sigma(t)\frac{{\cal
K}({\bm{p}},t;{\bm{p}}_{0},t_{0})}{{\cal
K}_{0}({\bm{p}};t-t_{0})}=\int_{\bm{q}}\,\sigma({\bm{q}})\frac{{\cal
K}({\bm{p}}-{\bm{q}},t;{\bm{p}}_{0},t_{0})}{{\cal
K}_{0}({\bm{p}};t-t_{0})}\,,$ (35)

where $\Sigma(t)=\int_{\bm{q}}\,\sigma({\bm{q}},t)$, as before. This motivates
defining the elastic Sudakov factor,

$\Delta(t,t_{0})\equiv{\rm e}^{-\int_{t_{0}}^{t}{\rm d}s\,\Sigma(s)}={\rm
e}^{-(t-t_{0})\Sigma}\,,$ (36)

where the last equality holds for media with constant density.777For medium
potentials with unscreened soft divergences, such as the HTL potential, one
has to modify this prescription to include an IR regulator. We refer to App. E
for a further discussion. This represents the probability of no elastic
scattering occurring between times $t_{0}$ and $t$. Integrating out the time,
we arrive at a slightly modified iterative equation

$\displaystyle{\cal K}({\bm{p}},t;{\bm{p}}_{0},t_{0})$
$\displaystyle=(2\pi)^{2}\delta({\bm{p}}-{\bm{p}}_{0})\Delta(t,t_{0}){\cal
K}_{0}({\bm{p}};t-t_{0})$ $\displaystyle+\int_{t_{0}}^{t}{\rm
d}s\frac{\Delta(t,t_{0})}{\Delta(s,t_{0})}\int_{\bm{q}}\,{\cal
K}_{0}({\bm{p}};t-s)\sigma({\bm{q}},s){\cal
K}({\bm{p}}-{\bm{q}},s;{\bm{p}}_{0},t_{0})\,.$ (37)

Compared to the standard opacity expansion in Eq. (3.1), this looks very
similar. However, the expansion is not in the potential $v({\bm{q}},s)$, which
contains both a “real” and a “virtual” part, but in the potential scattering
$\sigma({\bm{q}},s)$ which comes from the term that provides a finite momentum
transfer in the process. The virtual contributions, where no net momentum was
exchanged, are accounted for to all orders in the Sudakov factor. This is why
this expansion referred as a resummed opacity expansion.

The relevant scale that appears at high opacity is $\omega_{\text{\tiny
BH}}=\frac{1}{2}\mu^{2}\lambda$, as we discussed in Sec. 2. Interestingly, at
low opacity $L\ll\lambda$, the scale changes to
$\bar{\omega}_{c}=\frac{1}{2}\mu^{2}L$, which we recognize from the opacity
expansion. In this regime the ROE is actually equivalent to the OE if one
gathers all the terms up to the same order in $(\frac{L}{\lambda})^{n}$.
However, as we will see, the terms are reshuffled in the ROE compared to the
OE.

It is possible to reach a general formula for the resummed opacity expansion
at arbitrary order. This was done in App. A.2, and we refer to that section
for detailed calculations.

#### First order ($N_{r}=1$):

The first order can be obtained from Eq. (A.2) with $n=1$, and reads

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=1}}{{\rm
d}\omega}=-8\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\int_{{\bm{p}}}\,\tilde{\Sigma}({\bm{p}}^{2})\,{\rm
Im}\,T\left({\bm{p}}^{2}-i\frac{L}{\lambda}\right)\,,$ (38)

in re-scaled, dimensionless variables. Equation (38) corresponds to the
formula (Eq. (4.6)) in Ref. Andres:2020kfg , but is here derived more
rigorously. Here, we have defined the function

$\displaystyle T(x)=\int_{0}^{1}{\rm d}t_{1}\int_{0}^{t_{1}}{\rm d}t_{0}\,{\rm
e}^{-ix(t_{1}-t_{0})}=\frac{1-ix-{\rm e}^{-ix}}{x^{2}}\,.$ (39)

The real and imaginary parts of $T({\bm{p}}^{2}-i\chi)$ are given in Eq.
(A.2). After doing the angular integral this becomes

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=1}}{{\rm
d}\omega}=-2{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\int_{0}^{\infty}{\rm
d}p\,\frac{1}{p+\frac{\bar{\omega}_{c}}{\omega}}\,{\rm
Im}\,T\left(p-i\frac{L}{\lambda}\right)\,.$ (40)

At low opacity $L\ll\lambda$, the function $T(p-iL/\lambda)$ becomes

$\left.-{\rm Im}\,T(p)\right|_{L\ll\lambda}=\frac{p-\sin(p)}{p^{2}}\,,$ (41)

making it equivalent to the OE result in Eq. (26). The limiting behavior in
the relevant limits of Eq. (40) can be extracted, leading to

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=1}}{{\rm
d}\omega}\simeq\begin{cases}2{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\ln\left(\frac{\bar{\omega}_{c}}{\omega}\right)-1+\gamma_{E}\right)-{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\left(1-\pi\frac{\omega}{\bar{\omega}_{c}}\right),&\text{for
}\omega\ll\bar{\omega}_{c}\,,\\\\[0.5pt]
\frac{\pi{\rm\bar{\alpha}}}{2}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}-\frac{\pi}{6}{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\frac{\bar{\omega}_{c}}{\omega},&\text{for
}\omega\gg\bar{\omega}_{c}\,.\end{cases}$ (42)

At leading order in $\mathcal{O}(\frac{L}{\lambda})$ this is the same as the
$N=1$ opacity expansion, presented in Eq. (27). However, in contrast to the
OE, subleading “$N=2$”-like terms $\sim(\frac{L}{\lambda})^{2}$ appear, which
only will be relevant when compared to higher-order contributions at
$N_{r}=2$.

In the high opacity limit $L\gg\lambda$, we have to extract the relevant limit
of $T(p-i\chi)$ in a careful way, yielding

$\left.-{\rm
Im}\,T(p-i\chi)\right|_{L\gg\lambda}\simeq\frac{p}{(L/\lambda)^{2}+p^{2}}\,.$
(43)

Changing the integration variable to $q=p\lambda/L$, we observe that
$\omega_{\text{\tiny BH}}$ replaces $\bar{\omega}_{c}$ as the relevant scale,
and Eq. (40) becomes

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=1}}{{\rm d}\omega}$
$\displaystyle\simeq
2{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\omega_{\text{\tiny
BH}}}{\omega}\int_{0}^{\infty}{\rm d}q\,\frac{1}{q+\frac{\omega_{\text{\tiny
BH}}}{\omega}}\frac{q}{1+q^{2}}=2{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\omega_{\text{\tiny
BH}}}{\omega}\frac{\frac{\pi}{2}+\frac{\omega_{\text{\tiny
BH}}}{\omega}\ln\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)}{1+\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)^{2}}\,.$ (44)

Finally, one can extract the soft and hard limits of this expression, which
are given by

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=1}}{{\rm
d}\omega}\simeq\begin{cases}2{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\ln\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)+\frac{\pi}{2}\frac{\omega}{\omega_{\text{\tiny
BH}}}\right),&\text{for }\omega\ll\omega_{\text{\tiny BH}}\,,\\\\[0.5pt]
\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\omega_{\text{\tiny
BH}}}{\omega},&\text{for }\omega\gg\omega_{\text{\tiny BH}}\,.\end{cases}$
(45)

The soft limit agrees with the heuristic discussion in Sec. 2. Strikingly, we
see that the behavior in the soft and hard limit takes exactly the same form
as for $N=1$ except that $\bar{\omega}_{c}$ has been replaced by
$\omega_{\text{\tiny BH}}$ (note that
$\bar{\omega}_{c}(L=\lambda)=\omega_{\text{\tiny BH}}$).

#### Second order ($N_{r}=2$):

The second order is found from Eq. (A.2) with $n=2$, and reads

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm d}\omega}$
$\displaystyle=8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{2}\int_{{\bm{p}}_{2},{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{2}^{2})\frac{{\bm{p}}_{2}\cdot{\bm{p}}_{1}}{{\bm{p}}_{2}^{2}}\tilde{\sigma}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times\frac{1}{{\bm{p}}_{2}^{2}-{\bm{p}}_{1}^{2}}\left({\rm
Re}\,T({\bm{p}}_{1}^{2}-i\chi)-{\rm Re}\,T({\bm{p}}_{2}^{2}-i\chi)\right)\,.$
(46)

After going to polar coordinates and doing the angular integrals, this becomes

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm d}\omega}$
$\displaystyle=4{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{2}\int_{0}^{\infty}{\rm
d}p_{2}\int_{0}^{\infty}{\rm
d}p_{1}\,\frac{1}{p_{2}+\frac{\bar{\omega}_{c}}{\omega}}\frac{p_{1}}{\left[\left(p_{1}+p_{2}+\frac{\bar{\omega}_{c}}{\omega}\right)^{2}-4p_{1}p_{2}\right]^{3/2}}$
$\displaystyle\times\frac{1}{p_{2}-p_{1}}\left({\rm Re}\,T(p_{1}-i\chi)-{\rm
Re}\,T(p_{2}-i\chi)\right)\,.$ (47)

We study this expression separately in the low- and high-opacity limits.

In the low opacity limit $L\ll\lambda$, the spectrum becomes

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm d}\omega}$
$\displaystyle\simeq
4{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{2}\int_{0}^{\infty}{\rm
d}p_{2}\int_{0}^{\infty}{\rm
d}p_{1}\,\frac{1}{p_{2}+\frac{\bar{\omega}_{c}}{\omega}}\frac{p_{1}}{\left[\left(p_{1}+p_{2}+\frac{\bar{\omega}_{c}}{\omega}\right)^{2}-4p_{1}p_{2}\right]^{3/2}}$
$\displaystyle\times\frac{1}{p_{2}-p_{1}}\left(\frac{1-\cos{p_{1}}}{p_{1}^{2}}-\frac{1-\cos{p_{2}}}{p_{2}^{2}}\right)\,,$
(48)

where again the only relevant energy scale is $\bar{\omega}_{c}$, as it is in
the OE. The double momentum integral can be recognized as $\mathcal{I}_{2}$
from $N=2$ of the opacity expansion. The soft and hard limits are

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm
d}\omega}\simeq\begin{cases}\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\frac{\omega}{\bar{\omega}_{c}},&\qquad\text{for
}\omega\ll\bar{\omega}_{c}\,,\\\\[0.5pt]
\frac{\pi}{6}{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{2}\frac{\bar{\omega}_{c}}{\omega},&\qquad\text{for
}\omega\gg\bar{\omega}_{c}\,.\end{cases}$ (49)

Summing up the two first orders of the ROE and OE we see that $N_{r}=1+2$
agrees with $N=1+2$, but only when keeping the subleading
$\sim(L/\lambda)^{2}$ terms at order $N_{r}=1$. As mentioned before, the
opacity expansion is arranged so that the order $N=n$ only contains terms
where the opacity scales as $\sim\chi^{n}$, where $\chi=L/\lambda$. The
resummed opacity expansion also includes all of the same terms, but they are
spread out over different orders of the expansion due to the resummation
contained in the Sudakov factor. The orders $N_{r}<n$ do contain terms going
as $\chi^{n}$. To get the right term at order $\chi^{n}$ in the ROE one
therefore has to keep the subleading corrections going as $\chi^{n}$ at all
previous orders of the expansion. For this reason the opacity expansion is
more convenient to use in the low opacity limit, as it does not mix orders of
opacity.

In the high opacity limit $L\gg\lambda$, we get

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm d}\omega}$
$\displaystyle\simeq
4{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)^{2}\int_{0}^{\infty}{\rm d}p_{2}\int_{0}^{\infty}{\rm
d}p_{1}\,\frac{1}{p_{2}+\frac{\bar{\omega}_{c}}{\omega}}\frac{p_{1}}{\left[\left(p_{1}+p_{2}+\frac{\bar{\omega}_{c}}{\omega}\right)^{2}-4p_{1}p_{2}\right]^{3/2}}$
$\displaystyle\times\frac{p_{2}+p_{1}}{(1+p_{2}^{2})(1+p_{1}^{2})}\,.$ (50)

The soft and hard limits of this expression are given by

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=2}}{{\rm
d}\omega}\simeq\begin{cases}\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\omega}{\omega_{\text{\tiny
BH}}},&\qquad\text{for }\omega\ll\omega_{\text{\tiny BH}}\,,\\\
\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\omega_{\text{\tiny
BH}}}{\omega},&\qquad\text{for }\omega\gg\omega_{\text{\tiny
BH}}\,,\end{cases}$ (51)

where similarly to $N_{r}=1$, the relevant scale is now $\omega_{\text{\tiny
BH}}$. Both $N_{r}=1$ and 2 goes as $\sim\bar{\alpha}\frac{L}{\lambda}$,
however, in the soft limit $N_{r}=1$ dominates, while in the hard limit ${\rm
d}I^{N_{r}=2}\sim{\rm d}I^{N_{r}=1}$. This shows that ROE is quickly
convergent if $\omega\ll\omega_{\text{\tiny BH}}$, while the expansion appears
to break down for harder emissions. We expect this structure to appear to all
orders in $N_{r}=n$. The resulting validity of the expansion is shown in the
left of Fig. 5 in blue. Based on our findings, the expansion scheme for the
ROE at high opacity $L\gg\lambda$ is

$\omega\frac{{\rm d}I}{{\rm
d}\omega}=\bar{\alpha}\frac{L}{\lambda}\sum_{n=0}^{\infty}f_{n}\left(\frac{\omega}{\omega_{\text{\tiny
BH}}}\right)\,,$ (52)

where $f_{n}$ is a finite function that can be obtained order by order for
$\omega\ll\omega_{\text{\tiny BH}}$.

The resulting spectrum is shown in the right of Fig. 5 for different
propagation lengths. For short times ($t<\lambda$), the ROE is valid for all
$\omega$ and it gives the same spectrum as the OE (compare to the right panel
of Fig. 4). For longer propagation the ROE is only valid if
$\omega<\omega_{\text{\tiny BH}}$, which is denoted with bullets in the
figure. Outside of the valid region, the curves turn to gray (see also the
left panel). Based on the figure, $N_{r}=2$ has negligible contribution to the
spectrum until $t\approx\lambda$ or $\omega\approx\omega_{\text{\tiny BH}}$.
Again, the dashed line represents the full numerical evaluation of the
spectrum from Ref. Andres:2020kfg ; Andres:2022  which is well approximated by
the ROE in its region of validity.

Figure 5: Left: The sketch of region of validity (and convergence) of the
resummed opacity expansion for different propagation length $t$ and emission
energy $\omega$. Right: The induced emission spectrum for gluons in the
resummed opacity expansion. The gray part of the curves denotes regions, where
the expansion is not valid. With the parameters presented, $\lambda=0.06$ fm.
The full numeric solution is also presented with dashed lines.

###  3.4 Improved opacity expansion (IOE)

The final expansion scheme we consider is the improved opacity expansion,
introduced in Refs. Mehtar-Tani:2019tvy ; Mehtar-Tani:2019ygg ; Barata:2020sav
; Barata:2020rdn ; Barata:2021wuf . We saw that the ROE at high opacity does
not work for energies much higher than $\omega_{\text{\tiny BH}}$. This makes
sense since at higher energies the formation time ${t_{\text{f}}}$ of the
emission becomes bigger than the mean free path, implying that the parton will
scatter many times on the medium. The main motivation for the improved opacity
expansion is to resum multiple soft scatterings, while account perturbatively
for rare, hard scatterings. This is achieved by introducing a scale $Q^{2}$
that separates soft and hard scatterings in the scattering potential,

$v({\bm{x}},t)\approx v_{\text{\tiny HO}}({\bm{x}},t)+\delta v({\bm{x}},t)\,,$
(53)

where $v_{\text{\tiny HO}}({\bm{x}},t)=\frac{{\hat{q}}(t)}{4}{\bm{x}}^{2}$ and
$\delta
v({\bm{x}},t)=\frac{{\hat{q}}_{0}(t)}{4}{\bm{x}}^{2}\ln{\frac{1}{Q^{2}{\bm{x}}^{2}}}$.
Equation (53) is the $\mu|\bm{x}|\ll 1$ expansion of Eq. (14) with the GW
potential. The first term is referred to as the harmonic oscillator
approximation (HO), where the jet quenching parameter is

${\hat{q}}(t)={\hat{q}}_{0}(t)\ln{\frac{Q^{2}}{\mu_{*}^{2}}}\,,$ (54)

where $\mu_{\ast}^{2}=\frac{\mu^{2}}{4}{\rm e}^{-1+2\gamma_{E}}$ for the GW
potential. The logarithm in $\hat{q}$ comes from the fact that the typical
exchanged momentum
$\langle\bm{k}^{2}\rangle=L\,\int_{\bm{q}}\bm{q}^{2}\sigma(\bm{q})$ is
divergent and thus it has to be regulated resulting in the leading logarithmic
form in Eq. (54) (see in Ref. Arnold:2009mr ). As long as
$Q^{2}/\mu_{\ast}^{2}\gg 1/(Q^{2}\bm{x}^{2})$, the HO term dominates over
$\delta v$, and the latter can be treated as a perturbation. This provides a
big advantage, since the multiple scattering in the HO approximation can be
resummed analytically.

The separation scale $Q^{2}$ has to be fixed in a meaningful way to not to
interfere with the expansion Barata:2020sav . A natural choice that achieves
this is evaluating $Q^{2}$ at the typical transverse momentum of the emission
$\bm{k}^{2}\sim\hat{q}{t_{\text{f}}}$, that yields

$Q_{r}^{2}(\omega)=\sqrt{\hat{q}\omega}\,,$ (55)

which constitutes an implicit equation for $Q^{2}$ through $\hat{q}$ itself,
see Eq. (54).

The IOE corresponds to expanding the full medium solution for ${\cal
K}({\bm{x}};{\bm{y}})$ around the harmonic oscillator, in contrast to the
conventional opacity expansion where one expands around the vacuum solution.
It can be cast in the iterative equation

$\displaystyle{\cal K}({\bm{x}},t_{2};{\bm{y}},t_{1})$ $\displaystyle={\cal
K}_{\text{\tiny HO}}({\bm{x}},t_{2};{\bm{y}},t_{1})$
$\displaystyle-\int_{t_{1}}^{t_{2}}{\rm d}s\int_{\bm{z}}\,{\cal
K}_{\text{\tiny HO}}({\bm{x}},t_{2};{\bm{z}},s)\delta v({\bm{z}},s){\cal
K}({\bm{z}},s;{\bm{y}},t_{1})\,.$ (56)

Here, ${\cal K}_{\text{\tiny HO}}({\bm{x}};{\bm{y}})$ is itself the solution
to an iterative equation, namely

$\displaystyle{\cal K}_{\text{\tiny HO}}({\bm{x}},t_{2};{\bm{y}},t_{1})$
$\displaystyle={\cal K}_{0}({\bm{x}}-{\bm{y}},t_{2}-t_{1})$
$\displaystyle-\int_{t_{1}}^{t_{2}}{\rm d}s\,\int_{\bm{z}}\,{\cal
K}_{0}({\bm{x}}-{\bm{z}},t_{2}-s)v_{\text{\tiny HO}}({\bm{z}},s){\cal
K}_{\text{\tiny HO}}({\bm{z}},s;{\bm{y}},t_{1})\,.$ (57)

The formal solution for ${\cal K}_{\text{\tiny
HO}}({\bm{x}},t_{2};{\bm{y}},t_{1})$ can also be cast as a path integral,
namely

${\cal K}_{\text{\tiny
HO}}({\bm{x}},t_{2};{\bm{y}},t_{1})=\int_{{\bm{r}}(t_{1})={\bm{y}}}^{{\bm{r}}(t_{2})={\bm{x}}}{\cal
D}{\bm{r}}\,{\rm e}^{i\int_{t_{1}}^{t_{2}}{\rm
d}s\,\big{[}\frac{\omega}{2}\dot{\bm{r}}^{2}+iv_{\text{\tiny
HO}}({\bm{r}},s)\big{]}}\,,$ (58)

which has a well-known analytical solution in a static medium,

$\mathcal{K}_{\text{\tiny
HO}}\left(\bm{x},t_{2};\bm{y},t_{1}\right)=\frac{\omega\Omega}{2\pi
i\,\sin(\Omega\Delta t)}{\rm e}^{\frac{i\omega\Omega}{2\sin(\Omega\Delta
t)}\left[\cos(\Omega\Delta
t)\,({\bm{x}}^{2}+{\bm{y}}^{2})-2\bm{x}\cdot\bm{y}\right]}\,,$ (59)

where $\Delta t\equiv t_{2}-t_{1}$ and
$\Omega=\frac{1-i}{2}\sqrt{\hat{q}/\omega}$ is the characteristic oscillator
frequency.

Inserting this expansion into the equation for the medium-induced spectrum Eq.
(12) separates it into two parts,

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny HO}}}{{\rm d}\omega}$
$\displaystyle=\frac{2\alpha_{s}C_{R}}{\omega^{2}}{\rm
Re}\int_{0}^{\infty}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm
d}t_{1}\,\partial_{\bm{x}}\cdot\partial_{\bm{y}}\big{[}{\cal K}_{\text{\tiny
HO}}({\bm{x}},t_{2};{\bm{y}},t_{1})-{\cal
K}_{0}({\bm{x}},t_{2};{\bm{y}},t_{1})\big{]}_{{\bm{x}}={\bm{y}}=0}\,,$ (60)
$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny IOE}}}{{\rm d}\omega}$
$\displaystyle=-\frac{2\alpha_{s}C_{R}}{\omega^{2}}{\rm
Re}\int_{0}^{\infty}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm d}s\int_{0}^{s}{\rm
d}t_{1}\,\int_{\bm{z}}\,$
$\displaystyle\times\partial_{\bm{x}}\cdot\partial_{\bm{y}}[{\cal
K}_{\text{\tiny HO}}({\bm{x}},t_{2};{\bm{z}},s)\delta v({\bm{z}},s){\cal
K}({\bm{z}},s;{\bm{y}},t_{1})]_{{\bm{x}}={\bm{y}}=0}\,.$ (61)

The first term gives rise to the well-known HO spectrum, while the second
constitutes an expansion in hard splittings around the harmonic oscillator.
The IOE spectrum can be simplified further, giving

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny IOE}}}{{\rm
d}\omega}=\frac{2{\rm\bar{\alpha}}}{\omega}{\rm Re}\,i\int_{0}^{L}{\rm
d}t_{2}\int^{t_{2}}_{0}{\rm d}t_{1}\int_{\bm{x}}\,{\rm
e}^{-i\frac{\omega\Omega}{2}\tan{(\Omega(L-t_{2})){\bm{x}}^{2}}}\delta
v({\bm{x}})\frac{{\bm{x}}}{{\bm{x}}^{2}}\cdot\partial_{\bm{y}}{\cal
K}({\bm{x}},t_{2};{\bm{y}},t_{1})|_{{\bm{y}}=0}\,,$ (62)

where ${\cal K}({\bm{x}};{\bm{y}})$ should be iterated using Eq. (3.4) in
order to generate higher orders of the expansion. A general formula for the
improved opacity expansion at arbitrary order can also be derived, which was
done in Sec. A.3, see also in Ref. Barata:2020sav .

#### Harmonic oscillator (HO):

The harmonic oscillator approximation resums all coherent soft scatterings
during the formation of the emission. The relevant scale is

$\omega_{c}\equiv\frac{1}{2}\hat{q}L^{2}\,,$ (63)

which was discussed in Sec. 2. The HO approximation is expected to be valid
for $\omega_{\text{\tiny BH}}\ll\omega\ll\omega_{c}$.

We derive here the familiar harmonic oscillator spectrum in a new way. The
vacuum contribution can easily be subtracted by inserting Eq. (3.4) into Eq.
(12), which gives

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny HO}}}{{\rm d}\omega}=$
$\displaystyle-\frac{2\alpha_{s}C_{R}}{\omega^{2}}{\rm
Re}\int_{0}^{\infty}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm d}s\int_{0}^{s}{\rm
d}t_{1}\int_{\bm{z}}$
$\displaystyle\times\partial_{\bm{x}}\cdot\partial_{\bm{y}}\big{[}{\cal
K}_{0}({\bm{x}},t_{2};{\bm{z}},s)v_{\text{\tiny HO}}({\bm{z}},s){\cal
K}_{\text{\tiny
HO}}({\bm{z}},s;{\bm{y}},t_{1})\big{]}_{{\bm{x}}={\bm{y}}=0}\,.$ (64)

Using the fact that $\int_{s}^{\infty}{\rm d}t_{2}\,\partial_{\bm{x}}{\cal
K}_{0}({\bm{x}},t_{2};{\bm{z}},s)|_{{\bm{x}}=0}=-i\frac{\omega}{\pi}\frac{{\bm{z}}}{{\bm{z}}^{2}}$,
the spectrum becomes

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny HO}}}{{\rm d}\omega}=$
$\displaystyle\frac{{\rm\bar{\alpha}}{\hat{q}}}{2\omega}{\rm
Re}\,i\int_{0}^{L}{\rm d}t_{2}\int^{t_{2}}_{0}{\rm
d}t_{1}\int_{\bm{z}}\,{\bm{z}}\cdot\partial_{\bm{y}}{\cal K}_{\text{\tiny
HO}}({\bm{z}},t_{2};{\bm{y}},t_{1})|_{{\bm{y}}=0}\,.$ (65)

This can be further simplified by using that

$\displaystyle\int_{\bm{z}}\,{\bm{z}}\cdot\bm{\partial}_{\bm{y}}{\cal
K}_{\text{\tiny HO}}({\bm{z}},t_{2};{\bm{y}},t_{1})|_{{\bm{y}}=0}$
$\displaystyle=-\frac{(\omega\Omega)^{2}}{2\pi\sin^{2}(\Omega(t_{2}-t_{1}))}\int_{\bm{z}}\,{\bm{z}}^{2}{\rm
e}^{i\frac{\omega\Omega}{2}{\bm{z}}^{2}\cot(\Omega(t_{2}-t_{1}))}\,,$
$\displaystyle=\frac{2}{\cos^{2}(\Omega(t_{2}-t_{1}))}\,.$ (66)

The time integration can be now be dealt with straightforwardly, yielding

$\displaystyle\int_{0}^{L}{\rm d}t_{2}\int_{0}^{t_{2}}{\rm
d}t_{1}\,\frac{1}{\cos^{2}(\Omega(t_{2}-t_{1}))}=-\frac{\ln\cos\Omega
L}{\Omega^{2}}\,,$ (67)

and thus the spectrum becomes

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny HO}}}{{\rm d}\omega}$
$\displaystyle=2\bar{\alpha}\,\ln\left|\cos\Omega L\right|\,,$ (68)

which is the familiar BDMPS-Z spectrum Baier:1996kr ; Zakharov:1996fv . The
limits of this are

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny HO}}}{{\rm
d}\omega}\simeq\begin{cases}{\rm\bar{\alpha}}\sqrt{\frac{2\omega_{c}}{\omega}},&\text{for
}\omega\ll\omega_{c}\,,\\\\[0.7pt]
\frac{{\rm\bar{\alpha}}}{6}\left(\frac{\omega_{c}}{\omega}\right)^{2},&\text{for
}\omega\gg\omega_{c}\,.\end{cases}$ (69)

The soft limit agrees with the discussion in Sec. 2, while the hard limit is
subleading compared to the OE $N=1$ in Eq. (27). Defining $\hat{q}$ with a
logarithm extends the region of validity, which was also found in Ref.
Arnold:2009mr , leading to the curved $\omega_{c}(t)$ line in Fig. 3.

#### Next-to-harmonic oscillator (NHO):

Using the definition in Eq. (62) and the results of Sec. A.3, the first order
of the improved opacity expansion can be written as

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny NHO}}}{{\rm d}\omega}$
$\displaystyle=\frac{2{\rm\bar{\alpha}}}{\pi}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}{\rm
Re}\int_{0}^{1}{\rm
d}s\int_{\bm{u}}\,\frac{1}{2}\ln{\left(\frac{\omega}{\bar{\omega}_{c}}\frac{\mu^{2}}{2Q^{2}}\frac{1}{{\bm{u}}^{2}}\right)}{\rm
e}^{\frac{i}{2}f(s){\bm{u}}^{2}}\,,$ (70)

where we have defined the function
$f(s)=\sigma\sqrt{\omega_{c}/\omega}[\cot(\sigma
s\sqrt{\omega_{c}/\omega})-\tan(\sigma(1-s)\sqrt{\omega_{c}/\omega})]$, and
$\sigma=\frac{1-i}{\sqrt{2}}$. After doing the ${\bm{u}}$ integral this
becomes

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny NHO}}}{{\rm
d}\omega}=2{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}{\rm
Re}\,i\int_{0}^{1}{\rm
d}s\,\frac{1}{f(s)}\left[1-\gamma_{E}+\ln{\left(-i\frac{\omega}{\bar{\omega}_{c}}\frac{\mu_{*}^{2}}{Q^{2}}f(s)\right)}\right]\,.$
(71)

The limits of this expression can readily be extracted. In the soft limit,
$\omega\ll\omega_{c}$, we have $f(s)\to 2i\sigma\sqrt{\omega_{c}/\omega}$
while in the hard limit, $\omega\gg\omega_{c}$, it becomes $f(s)\to 1/s$.
These simplifications make it possible to do the last time integration. Hence,
the extracted limiting behavior is,

$\displaystyle\omega\frac{{\rm d}I^{\text{\tiny NHO}}}{{\rm
d}\omega}\simeq\begin{cases}{\rm\bar{\alpha}}\sqrt{\frac{2\omega_{c}}{\omega}}\,\frac{1}{2\ln
Q^{2}/\mu_{\ast}^{2}}\left(\frac{\pi}{4}+\gamma_{E}+\ln{\left(\frac{\sqrt{{\hat{q}}\omega}}{\sqrt{2}Q^{2}}\right)}\right),&\qquad\text{for
}\omega\ll\omega_{c}\,,\\\
\frac{\pi{\rm\bar{\alpha}}}{2}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega},&\qquad\text{for
}\omega\gg\omega_{c}\,.\end{cases}$ (72)

In the soft limit $\omega\ll\omega_{c}$, NnHO terms will take the form of the
HO by using $Q_{r}$, and thus

$\omega\frac{{\rm d}I}{{\rm
d}\omega}={\rm\bar{\alpha}}\sqrt{\frac{\hat{q}L^{2}}{\omega}}\left(1+\frac{1}{2}\frac{a_{0}}{\ln
Q_{r}^{2}/\mu^{2}_{\ast}}+\mathcal{O}\left(\frac{1}{\ln
Q_{r}^{2}/\mu_{\ast}^{2}}\right)^{2}\right)\,,$ (73)

where we added the HO term, and used Eq. (55). The choice of $Q=Q_{r}(\omega)$
is effective, when the medium is big enough $L\gg\lambda$. It is clear that
the expansion parameter of the IOE is $\ln^{-1}(Q_{r}^{2}/\mu_{\ast}^{2})\ll
1$ in the soft limit. Therefore, NnHO terms can be absorbed into an effective
jet transport parameter,

$\hat{q}_{\text{eff}}(Q^{2})=\hat{q}_{0}\ln\left(\frac{Q_{r}^{2}}{\mu^{2}_{\ast}}\right)\left[1+\frac{a_{0}}{\ln
Q_{r}^{2}/\mu_{\ast}^{2}}+\frac{a_{1}}{\ln^{2}Q_{r}^{2}/\mu_{\ast}^{2}}+\dots\right]\,.$
(74)

The coefficients $a_{0}=1.016$ and $a_{1}=0.316$ of the expansion and higher-
order terms up to N2HO were found in Ref. Barata:2020sav .

In the hard limit of Eq. (72), one can see that the IOE reproduces the hard
limit of $N=1$ in the OE from Eq. (27). Furthermore, it is bigger than the HO
contribution in Eq. (68) and thus NHO dominates for $\omega\gg\omega_{c}$.

Figure 6: Left: The sketch of region of validity (and convergence) of the
improved opacity expansion for different propagation length $t$ and emission
energy $\omega$. Right: The induced emission spectrum for gluons in the
improved opacity expansion. The gray part of the curves denotes regions where
the expansion is not valid. Using the parameters presented, $\lambda=0.06$ fm.
The full numeric solution is also presented with dashed lines.

As mentioned above, the HO is meaningful if $L>\lambda$ and
$\omega>\omega_{\text{\tiny BH}}$. Furthermore, the IOE is expected to
converge if $v_{\text{\tiny HO}}>\delta v$ or equivalently
$\omega>\omega_{\text{\tiny BH}}$. However, as the hard limits of the first
order of the IOE and the OE are equal, and the OE is valid down to
$\bar{\omega}_{c}$, it is reasonable to assume that also the IOE is valid down
to $\bar{\omega}_{c}$. Therefore, the region of validity will be extended to
all $L$ and $\omega>\text{min}(\omega_{\text{\tiny BH}},\bar{\omega}_{c})$, as
shown in the left of Fig. 6 in red.

The spectra obtained with the IOE from Eq. (68) and Eq. (72) are shown in the
right panel of Fig. 6 for different propagation lengths. At early times
$t<\lambda$, the HO approximation is highly suppressed, due to the absence of
multiple scattering. However, the contribution from the NHO makes the total
agree with the hard limit of $N=1$ OE and $N_{r}=1$ ROE (cf. Figs. 4 and 5).
The deviation close to $\omega\approx\bar{\omega}_{c}$ (gray bullet) arises
since $Q_{r}$ was chosen to reproduce the HO spectrum which is strictly valid
for $L\gg\lambda$. For later times $t>\lambda$, both the HO and NHO will give
sizable contributions, where the HO dominates if $\omega<\omega_{c}$ (red
bullets) and NHO dominates if $\omega>\omega_{c}$. The HO approximation breaks
down if $\omega<\omega_{\text{\tiny BH}}$ (gray bullets). The dashed line is
the full numerical solution from Refs. Andres:2020kfg ; Andres:2022  and the
IOE well captures it in its region of validity.

###  3.5 Summary of the regimes and the induced emission spectrum

In this section, we have presented three distinct perturbative expansions (OE,
ROE, and IOE) that provide different ways of calculating the induced emission
spectrum $\omega\frac{{\rm d}I}{{\rm d}\omega}$ in their respective regimes of
convergence. The OE and IOE are expansion schemes that were first developed in
previous works, while the ROE is rigorously derived in this section for the
first time. Their regions of validity are sketched previously in Figs. 4–6,
and at least one of the expansions is valid at every point in the phase space
$(\omega,t)$. Here $\omega$ is the emitted energy and $t$ is the propagated
length ($L$ is the maximal length of the medium and $E$ is the energy of the
emitting particle). As a consequence, our description of the spectrum is
complete in the full phase space, as one can always use one of the expansions
to reach an approximation of the true spectrum, and one can reach better
accuracy by including higher orders. Note that the expansions are overlapping:
for $L<\lambda$ both OE and ROE are valid, and for $\omega>\omega_{c}$ both
IOE the OE can be used.

Figure 7: Summary of the induced emission spectrum for gluons, combined from
the three expansion schemes at different propagation length (different
panels). The black lines are our final forms from Eq. (75) that uses $N_{r}=1$
and HO+NHO. The shaded areas denote the leading scattering process and the
corresponding dotted lines are the limiting $N=1$, HO and $N_{r}=1$
contributions. The dashed lines are the numerical solution of Eq. (12) from
Ref. Andres:2020kfg .

To describe the spectrum in the whole phase space, we use (to first order)

$\frac{{\rm d}I^{\rm Full}}{{\rm d}\omega}=\begin{cases}\frac{{\rm d}I^{\rm
ROE}}{{\rm d}\omega}\,,&\omega<\text{min}(\omega_{\text{\tiny
BH}},\bar{\omega}_{c}(t))\,,\\\ \frac{{\rm d}I^{\rm IOE}}{{\rm
d}\omega}\,,&{\rm otherwise}\,.\end{cases}$ (75)

Based on Figs. 4–6 (and the all order expansion formulas), the first-order
terms already capture the most important effects.

To summarize, the limiting behaviour of the spectrum in different regions of
the phase space is

$\left.\omega\frac{{\rm d}I}{{\rm
d}\omega}\right|_{L\ll\lambda}=\begin{cases}2{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\ln\frac{\bar{\omega}_{c}}{\omega}-1+\gamma_{E}\right)\,,&\quad\text{for
}\omega\ll\bar{\omega}_{c}\,,\\\
\frac{\pi}{2}{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\,,&\quad\text{for
}\bar{\omega}_{c}\ll\omega\,,\end{cases}$ (76)

for $L\ll\lambda$, and

$\left.\omega\frac{{\rm d}I}{{\rm
d}\omega}\right|_{L\gg\lambda}=\begin{cases}2{\rm\bar{\alpha}}\frac{L}{\lambda}\ln\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)\,,&\quad\text{for }\omega\ll\omega_{\text{\tiny
BH}}\,,\\\
{\rm\bar{\alpha}}\sqrt{\frac{2\omega_{c}}{\omega}}\,,&\quad\text{for
}\omega_{\text{\tiny BH}}\ll\omega\ll\omega_{c}\,,\\\
\frac{\pi}{2}{\rm\bar{\alpha}}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\,,&\quad\text{for
}\omega_{c}\ll\omega\,,\end{cases}$ (77)

for $L\gg\lambda$. This agrees with the formulas from the heuristic discussion
in Sec. 2. In Fig. 7, we evaluated Eq. (75) (black curve) up to $N_{r}=1$ and
HO+NHO for different times. The dotted curves are the limits of $N_{r}=1$, HO
and $N=1$, shown in blue, red and green respectively. The regions are shaded
with the same colors as in Fig. 3, visualizing the regions of the distinct
scattering processes. At the transition point $\text{min}(\omega_{\text{\tiny
BH}},\bar{\omega}_{c}(t))$, the spectrum is not completely smooth, and the
difference is expected to vanish as one goes to higher orders in the
perturbative expansion. We defined a switching function that makes the
transition smoother, which is described in App. D.

In Fig. 7, we have also plotted the full numerical spectra from Ref.
Andres:2020kfg ; Andres:2022  with dashed lines. The excellent agreement with
our curves corroborates the validity of our formula in Eq. (75), where we used
$N_{r}=1$ and HO+NHO in the plot. In the ratio panel one can see, how the
different regions are adding up to give an overall very accurate description
of the full spectrum. Moreover, as we saw in Figs. 4–6, by including higher
orders e.g. $N_{r}=2$, the description becomes more accurate, smoothing the
transition around $\omega_{\text{\tiny BH}}$.

##  4 Resumming multiple emissions in the medium

Section 3 presents a theoretical framework consisting of different
perturbative expansions (namely the opacity expansion (OE), the resummed OE
(ROE), and the improved OE (IOE)) to describe the full phase space
($\omega,t$) of the medium-induced gluon emission spectrum. The emitted energy
$\omega$ is limited by the energy of the emitter $\omega\ll E$, and the
propagation time (length) $t$ is in turn limited by the medium length $t<L$.
Our effective formalism accounts for arbitrary many scatterings that can be
arbitrarily hard or soft. Finally, Eq. (75) describes the emission spectrum up
to arbitrary precision and recover the full solution of Eq. (12) that has only
been achieved numerically before CaronHuot:2010bp ; Feal:2018sml ;
Andres:2020vxs .

It is now time to explore what consequences the full induced spectrum instills
on a parton propagating through the medium. In this section we will present
analytical solutions of the evolution equation of the medium-induced cascade.
We present the inclusive gluon energy distribution for different medium
lengths, and we will focus for gluons for transparency. The numerical
solutions are presented in Sec. 5. Vacuum emissions belong outside of the
scope of the current work.

###  4.1 The necessity of multiple emissions

Multiple emissions have to be taken into account whenever the multiplicity of
gluons is large. We define the multiplicity of gluons above the energy
$\omega$ in terms of the spectrum ${\rm d}I/{\rm d}\omega$, as

$N(\omega)=\int_{\omega}^{\infty}{\rm d}\omega^{\prime}\frac{{\rm d}I}{{\rm
d}\omega^{\prime}}\,.$ (78)

The upper limit of the integral is taken to infinity because we will currently
assume that the energy of the emitter $E$ is much larger than the largest
available medium energy scale. Also, for our current purposes, it suffices to
consider the leading behavior of the spectrum in the various scattering
regimes presented in Fig. 3.

At low opacity $L\ll\lambda$ and starting from $\omega<\omega_{\text{\tiny
BH}}$, we find

$N(\omega)\simeq\bar{\alpha}\frac{L}{\lambda}\ln^{2}\frac{\bar{\omega}_{c}}{\omega}+\frac{\pi\bar{\alpha}}{2}\frac{L}{\lambda}\,.$
(79)

where $\bar{\omega}_{c}=\frac{1}{2}\mu^{2}L$ and we have only kept the leading
terms. The maximal multiplicity in the hard regime, see the second term, is
always small for perturbative splittings with $\bar{\alpha}\ll 1$.
Furthermore, the multiplicity in the soft (Bethe–Heitler) regime, given by the
first term, becomes large only at very small energies, i.e. $\omega<{\rm
e}^{\sqrt{1/(\bar{\alpha}L/\lambda)}}\bar{\omega}_{c}$. We can therefore
safely neglect multiple emissions at low opacities.

For dense media, $L\gg\lambda$, we discuss the three pertinent cases. As
usual, we will denote $\omega_{c}=\frac{1}{2}\hat{q}L^{2}$. Then, for
$\omega\gg\omega_{c}$ (rare hard scattering regime), we find

$N(\omega)\simeq\bar{\alpha}\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}\,,$
(80)

where we again neglected subleading terms. Hence hard emissions, described by
the OE expansion, can safely be considered to be rare.

Next, for $\omega_{\text{\tiny BH}}\ll\omega\ll\omega_{c}$ we get

$N(\omega)\simeq
2^{\frac{3}{2}}\bar{\alpha}\sqrt{\frac{\omega_{c}}{\omega}}\,,$ (81)

where we introduced the leading behavior of the spectrum and neglected the
multiplicity from hard emissions, following the discussion above. The
multiplicity becomes large $N(\omega)\gg 1$ at energies
$\omega\ll\bar{\alpha}^{2}\omega_{c}$ and thus multiple emission becomes
dominant. For large enough medium length $L\gg\lambda$, there is a significant
phase space allowing for multiple emission
($\bar{\alpha}_{s}^{2}\omega_{c}(t)>\omega_{\text{\tiny BH}}$), resulting in
the power enhancement.

Finally, for soft gluon energies $\omega<\omega_{\text{\tiny
BH}}\ll\omega_{c}$, the multiplicity is

$N(\omega)\simeq\bar{\alpha}\frac{L}{\lambda}\ln^{2}\frac{\omega_{\text{\tiny
BH}}}{\omega}+2^{\frac{3}{2}}\bar{\alpha}\sqrt{\frac{\omega_{c}}{\omega_{\text{\tiny
BH}}}}\,,$ (82)

where, again, only the leading terms from each regime were kept. The second
term in Eq. (82) scales as $\sim(L/\lambda)\sqrt{\hat{q}/\hat{q}_{0}}$. Based
on the discussion above, this term is already large and the multiplicity
continues to grow only logarithmically for small $\omega$, and therefore
multiple emissions are going to happen.

###  4.2 Resummation of multiple emissions

Considering multiple emissions in a medium poses a tremendous theoretical
challenge. The situation is quite analogous to the description of multiple
gluon emissions in vacuum. Similar to QCD jets, the main challenge when
considering medium effects lies in dealing with intricate interference effects
between subsequent emissions, see e.g. Refs. Arnold:2015qya ; Arnold:2020uzm ;
Arnold:2021pin . However, when considering strong ordering in the energies of
the subsequent emissions, these effects are suppressed Blaizot:2013vha .

In order to clarify the framework that we work in, let us briefly recall the
main arguments for neglecting interference effects for a set of multiple
induced emissions. For the time being, we stick to emissions in the HO region
which dominate the multiplicity. The typical time it takes an emission to
form, often referred to as a formation time (or in some works branching time),
of a soft gluon is ${t_{\text{f}}}\sim\sqrt{\omega/\hat{q}}$. This time is
much smaller than the extent of the medium ${t_{\text{f}}}\ll L$ as long as
$\omega\ll\omega_{c}$.

Another relevant quantity is the time between two subsequent emissions. This
is related to the no emission probability (or Sudakov factor). For a leading
particle with energy $E$, the first emission is produced at time
${t_{\text{f}}}_{1}$, with energy $\omega_{1}$. A second, strongly ordered
emission ($\omega_{2}\ll\omega_{1}$) of the original parton forms much quicker
${t_{\text{f}}}_{2}\ll{t_{\text{f}}}_{1}$. The time between the two emissions
$t_{\text{split}}$ can be estimated with
$\int_{{t_{\text{f}}}_{1}}^{t_{\text{split}}}{\rm
d}t\int_{\omega_{2}}^{\omega_{1}}{\rm d}\omega\,\frac{{\rm d}I}{{\rm
d}\omega{\rm d}t}\sim 1$, that is basically the probability of not having
emissions between $\omega_{2}<\omega<\omega_{1}$, resulting in

$t_{\text{split}}\sim{t_{\text{f}}}_{1}+\frac{1}{\bar{\alpha}}{t_{\text{f}}}_{2}\,.$
(83)

Hence, our rough estimate implies that typically
$t_{\text{split}}\gg{t_{\text{f}}}_{2}$ (for $\bar{\alpha}\ll 1$). Therefore,
the formation of emissions is short compared to the time that separates
emissions and thus emissions form independently. This motivates the
resummation of multiple independent emissions in terms of a rate equation.

Similar analysis can be done for hard emissions $\omega\gg\omega_{c}$, for
which $t_{\rm split}$ is very long due to the unlikeliness of hard
scatterings, and therefore the emissions are formed independently.

For soft emissions $\omega\ll\omega_{\text{\tiny BH}}$, $t_{\rm
split}\sim{t_{\text{f}}}_{1}+\frac{1}{\bar{\alpha}}\ln^{-2}\frac{\lambda}{{t_{\text{f}}}_{2}}$
and thus emissions form independently. Close to the boundary in cases where
$t_{\rm split}\approx{t_{\text{f}}}_{1}+{t_{\text{f}}}_{2}$, a more
complicated structure appears in terms of resummation, as emissions might
overlap. In this case, interference effects between the two emissions have to
be included. A similar thing happens in vacuum for wide angle soft emissions,
which result in angular ordering and in non-global effects for which the
resummation has been understood just recently Banfi:2021xzn . We will use the
rate equation to account for emissions with any $\omega$. However, it will not
necessarily account correctly for interference among them and further study is
needed in the future.

In this section we will focus on the single-inclusive energy distribution of
partons carrying energy $xE$ after travelling length $t$ in the medium, where
$E$ is the initial energy. It is defined as

$D(x,t)\equiv x\frac{{\rm d}N}{{\rm d}x}\,.$ (84)

The formalism can easily be extended to account for parton flavors, see e.g.
Mehtar-Tani:2018zba , but for now we restrict our attention to a pure gluon
cascade.

In Sec. 3 we focused on emissions of soft gluons with energies $\omega\ll E$.
Now we will consider generic splitting processes where a parton with flavor
index $a=q,g$ and initial energy $E$ shares its energy with two daughter
partons, with energies $zE$ and $(1-z)E$ and flavor indices $b$ and $c$,
respectively, for $0<z<1$. The spectrum of such splittings ${\rm d}I_{ba}/{\rm
d}z$ is given by Eq. (138). The general features of Fig. 3 remain the same
with the substitution $\omega\to z(1-z)E$.888Hence the upper limit of $\omega$
in Fig. 3 should now be $E/4$. For further details, see the discussion in App.
C.

As was shown in Ref. Blaizot:2012fh ; Blaizot:2013vha , for sufficiently soft
emissions, with formation times much smaller than the medium length,
interference effects are suppressed and one can consider multiple emissions as
occurring independently. The evolution equation for the energy distribution,
that accounts for an arbitrary number of induced emissions, is given by

$\displaystyle\frac{\partial}{\partial t}D(x,t)$
$\displaystyle=\int_{x}^{1}{\rm d}z\,{\cal
K}\left(z,\frac{x}{z}E,t\right)D\left(\frac{x}{z},t\right)-\int_{0}^{1}{\rm
d}z\,z\,{\cal K}\left(z,xE,t\right)D(x,t)\,.$ (85)

The initial condition is a single gluon carrying energy $E$, hence
$D(x,0)=\delta(1-x)$. The splitting kernel ${\cal K}(z,E,t)$ is the rate of
emissions off a particle with energy $E$,

${\cal K}(z,E,t)=2\left.\frac{{\rm d}I_{gg}}{{\rm d}z{\rm d}t}\right|_{E}\,,$
(86)

and the rate with full $z$-dependence for the $g\to gg$ splitting can be found
in Eq. (138).999In the soft limit, the kernel is closely related to the
spectrum discussed in the previous section or, more precisely, the rate ${\rm
d}I/({\rm d}\omega{\rm d}t)$, calculated in App. B. Importantly, for the gluon
splitting kernel, the divergences in $z\to 0$ and $z\to 1$ are folded together
in the limit $\omega\to 0$, hence the additional symmetry factor in Eq. (86).
The first term in Eq. (85) is a real emission describing an emitted gluon with
energy fraction $x$ (gain term), while the second is a virtual emission that
does not change the energy of the emitter (loss term). Both terms contribute
to cancelling out the apparent divergence at $z\to 1$. The evolution equation
conserves the total energy contained in the spectrum,

$\int_{0}^{1}{\rm d}x\,D(x,t)=1\,,$ (87)

which can be confirmed directly from (85).

The leading parametric behavior of the splitting kernels can be derived by
taking appropriate limits and is presented in Sec. 3, cf. Eqs. (27), (45), and
(69). This results in,

$\left.{\cal
K}(z,E,t)\right|_{t\ll\lambda}=\begin{cases}\frac{2{\rm\bar{\alpha}}}{z(1-z)}\frac{1}{\lambda}\ln\left(\frac{\bar{\omega}_{c}(t)}{z(1-z)E}\right)&\quad\text{for
}z(1-z)E\ll\bar{\omega}_{c}(t)\,,\\\
\frac{{\rm\bar{\alpha}}\pi}{2}\frac{{\hat{q}}_{0}t}{[z(1-z)]^{2}E}&\quad\text{for
}\bar{\omega}_{c}(t)\ll z(1-z)E\,,\end{cases}$ (88)

for $t\ll\lambda$, and

$\left.{\cal
K}(z,E,t)\right|_{t\gg\lambda}=\begin{cases}\frac{2{\rm\bar{\alpha}}}{z(1-z)}\frac{1}{\lambda}\ln\left(\frac{\omega_{\text{\tiny
BH}}}{z(1-z)E}\right)&\quad\text{for }z(1-z)E\ll\omega_{\text{\tiny BH}}\,,\\\
{\rm\bar{\alpha}}\sqrt{\frac{\hat{q}}{[z(1-z)]^{3}E}}&\quad\text{for
}\omega_{\text{\tiny BH}}\ll z(1-z)E\ll\omega_{c}(t)\,,\\\
\frac{{\rm\bar{\alpha}}\pi}{2}\frac{{\hat{q}}_{0}t}{[z(1-z)]^{2}E}&\quad\text{for
}\omega_{c}(t)\ll z(1-z)E\,,\end{cases}$ (89)

for $t\gg\lambda$, where $\bar{\omega}_{c}(t)=\frac{1}{2}\mu^{2}t$ and
$\omega_{c}=\frac{1}{2}\hat{q}t^{2}$. For the analytical estimates in this
section we neglect the running of $\hat{q}$, but this will be included in the
numerics presented in Sec. 5.

Figure 8: The medium-induced rate from Eqs. (86), and (75) for gluons for
different emitted energy and propagation time (black dashed lines). The black
solid lines include finite-$z$ corrections. The dotted lines are the
$N_{r}=1$, HO and $N=1$ lines according to Fig. 7. The color shading
corresponds to Fig. 3 i.e. denoting the dominant scattering processes. The
parameters correspond to $\lambda=0.06$ fm.

The rate ${\rm d}I/({\rm d}z{\rm d}t)$ is plotted in Fig. 8, in a similar
manner to the spectrum in Fig. 7, and rescaled by a factor $z^{3/2}$ to
highlight the behavior at small $z$. The panels in the upper part show the
$z$-dependence for three different $t$ (early, mid-, and late times), and has
a very similar structure to that of the spectrum. The color shading
corresponds to Fig. 3, i.e. the dominant scattering processes. The dotted
lines are the $N_{r}=1$, HO and $N=1$ lines from Sec. 3 as in Fig. 7, where
$\omega\to zE$ was used. Similarly, the dashed line is the full solution in
the soft limit with $\omega\to zE$. The solid black lines are from Eq. (107)
and they contain finite-$z$ corrections from App. C. The finite-$z$
corrections change the rate for hard emissions around $z\sim 1$. The panels in
the lower part show the time dependence of the rate for a fixed emitted
energy.

###  4.3 Analytic solutions of the evolution equation

The evolution equation (85) is readily solved by numerical evaluation, which
will be discussed in Sec. 5. Here we will discuss limiting cases where
analytical solutions can be found. We can find such solutions at early times
(considering only one emission) and at late times (considering many soft
emissions). These are by now well-known limiting cases. Finally, we also
consider the novel case of the evolution equation at intermediate times, where
both rare, hard emissions and multiple, soft emissions can occur in sequence
according to their respective allowed phase space of emissions, given in Fig.
3.

To simplify our discussion, in this section we will neglect the Bethe–Heitler
regime. We will nevertheless include it in the full numerical solutions
presented in Sec. 5.

#### Early time evolution:

At an early stage of medium propagation, the leading parton has little time to
interact with the medium which also translates into a small probability of
splitting. Given our previous discussion, the natural medium scale to compare
with is the mean free path $\lambda$. Hence, at $t\lesssim\lambda$ we consider
a single splitting, leading to

$\displaystyle D(x,t)$ $\displaystyle\simeq\delta(1-x)\left[1-\int_{0}^{t}{\rm
d}s\int_{0}^{1}{\rm d}z\,z{\cal K}(z,xE,s)\right]+\int_{0}^{t}{\rm
d}s\int_{x}^{1}{\rm d}z\,{\cal
K}\left(z,\frac{x}{z}E,s\right)\delta\left(1-\frac{x}{z}\right)$
$\displaystyle=\int_{0}^{t}{\rm d}s\,x{\cal K}(x,E,s)\,,$ (90)

where we dropped the term proportional to $\delta(1-x)$, which is only
important for energy conservation, cf. Eq. (87). Using the results in Eq.
(88), we find that

$D(x,t)\simeq\begin{cases}2{\rm\bar{\alpha}}\frac{t}{\lambda}\frac{1}{1-x}\ln\left(\frac{\bar{\omega}_{c}(t)}{x(1-x)E}\right)&\quad\text{for
}x\ll\bar{x}_{c}\,,\\\
\frac{\pi{\rm\bar{\alpha}}}{4}\frac{{\hat{q}}_{0}}{E}\frac{t^{2}}{x(1-x)^{2}}&\quad\text{for
}\bar{x}_{c}\ll x\ll 1-\bar{x}_{c}\,,\\\ \end{cases}$ (91)

where we have defined $\bar{x}_{c}=\bar{\omega}_{c}(t)/E$. For
$x\ll\bar{x}_{c}$, corresponding to the single soft or Bethe–Heitler
scattering regime, a characteristic $D\sim\ln 1/x$ structure appears that is
similar to the DGLAP energy distribution in vacuum. On the contrary, for
$x\gg\bar{x}_{c}$, corresponding to the single hard scattering regime, $D\sim
1/x$, and the two regimes are separated by $\bar{x}_{c}$. The limiting case of
Eq. (91) is shown in Fig. 9 with dashed lines for different times (different
panels). Since we, as we move onward, will largely neglect the description of
the infared regime, we only used the single hard scattering (green)
contribution that is valid $x>\bar{x}_{c}$.

Formally, the early-time expansion breaks down when $t>\lambda$ which is also
the characteristic time when multiple interactions with the medium become
important. Finally, we do not expect the early time solution to hold for for
$x<{\rm\bar{\alpha}}^{2}\omega_{c}/E$ where multiple emissions play an
important role.

#### Late time evolution:

The evolution equation can be solved exactly if one assumes coherent
scatterings dominate for all momentum fractions $x$ and all times $t$. This
approximation is most sound when $t\gg L_{c}$, as seen in Fig. 3. The
analytical solution neglects the Bethe–Heitler, which will make the solution
less reliable for very small $x<\omega_{\text{\tiny BH}}/E$. We will call this
solution $D_{0}$ as it serves as a baseline for subsequent calculations. It is
the solution to the evolution equation

$\displaystyle\frac{\partial}{\partial t}D_{0}(x,t)$
$\displaystyle=\int_{x}^{1}{\rm d}z\,{\cal K}_{\rm
coh}\left(z,\frac{x}{z}E,t\right)D_{0}\left(\frac{x}{z},t\right)-\int_{0}^{1}{\rm
d}z\,z{\cal K}_{\rm coh}\left(z,xE,t\right)D_{0}(x,t)$ (92)

where ${\cal K}_{\rm coh}$ is the soft limit of the harmonic oscillator
regime, see the middle term in Eq. (89). The solution from Ref. Blaizot:2013hx
is

$D_{0}(x,\tau)=\frac{\tau}{\sqrt{x}(1-x)^{3/2}}{\rm
e}^{-\pi\frac{\tau^{2}}{1-x}}\,,$ (93)

where we dropped the $\delta(x)$ zero mode term, which is only important for
energy conservation. We have also defined the re-scaled evolution variable
which absorbs the energy scale,

$\tau=\bar{\alpha}\sqrt{\frac{\hat{q}_{0}}{E}}t\,.$ (94)

The solution has the $D_{0}\sim\tau/\sqrt{x}$ shape characteristic of the
turbulent cascade. Equation (93) is shown in Fig. 9 with dotted lines for
different times.

Figure 9: Analytic solutions of the energy distribution in different
approximations. With these parameters $\lambda=0.06$ fm, and thus the
different panels are $t<\lambda$, $\lambda<t<L_{c}$, and $t\lesssim L_{c}$.

#### Intermediate time evolution:

To reach an approximate solution for intermediate times it is useful to cast
the evolution in different variables. Following Ref. Blaizot:2015jea , and
defining $\xi=x/z$ in the gain term and $\xi=xz$ in the loss term the
evolution can be rewritten as

$\frac{\partial}{\partial\tau}D(x,\tau)=\int_{x}^{1}{\rm
d}\xi\,P(x,\xi,\tau)D(\xi,\tau)-D(x,\tau)\int_{0}^{x}{\rm
d}\xi\,P(\xi,x,\tau)\,,$ (95)

where
$P(x,\xi,\tau)=\frac{1}{{\rm\bar{\alpha}}}\sqrt{\frac{E}{{\hat{q}}_{0}}}\frac{x}{\xi^{2}}{\cal
K}(\frac{x}{\xi},\xi E,t)$ and $\tau$ is defined in (94). This form of the
evolution equation clearly shows the cancellation of divergences between the
gain and loss terms at $\xi\to x$.

At early times, $t<\lambda$, the emissions are governed by single interactions
with the medium. In the regime $\lambda<t<L_{c}$ you also have to take into
account that coherent emissions play a role in the soft regime
$\omega<\omega_{c}(t)$ (neglecting Bethe–Heitler emissions). However, the
region above $\omega_{c}(t)$ will contain emissions from single hard
scatterings, see Eq. (89) (third line). In that case the kernel will be
modified to include both a soft and a hard component, divided by the
$\omega_{c}(t)$ separation line, as follows

$P(x,\xi,\tau)=\theta_{c}(x,\xi,\tau)P_{\rm
coh}(x,\xi,\tau)+\theta_{h}(x,\xi,\tau)P_{\rm hard}(x,\xi,\tau)\,,$ (96)

where the coherent and hard emission kernels were defined in Eq. (89). After
changing variables they become

$\displaystyle P_{\rm coh}(x,\xi,\tau)$
$\displaystyle=\sqrt{\frac{\xi}{x}}\frac{1}{(\xi-x)^{3/2}}\qquad\text{for
}x\ll x_{c}\text{ or }\xi-x_{c}\ll x\,,$ $\displaystyle P_{\rm
hard}(x,\xi,\tau)$
$\displaystyle=\frac{\pi}{2{\rm\bar{\alpha}}}\,\frac{\xi}{x}\frac{1}{(\xi-x)^{2}}\tau\qquad\text{for
}x_{c}\ll x\ll\xi-x_{c}\,,$ (97)

where $x_{c}\equiv\omega_{c}(t)/E=\tau^{2}/(2{\rm\bar{\alpha}}^{2})$ and we
have assumed that $x_{c}\ll 1$ in the limits. These conditions can be encoded
in a set of Heaviside theta-functions to make sure that each kernel is used
solely in its regime of validity, namely 101010In the case where $x_{c}$ is
not small the conditions are slightly more complicated, and the hard regime is
encoded in
$\theta_{h}(x,\xi,\tau)=\Theta\left(x-\frac{\xi}{2}\left(1-\sqrt{1-\frac{4x_{c}}{\xi}}\right)\right)\Theta\left(\frac{\xi}{2}\left(1+\sqrt{1-\frac{4x_{c}}{\xi}}\right)-x\right)$,
and $\theta_{c}=1-\theta_{h}$.

$\displaystyle\theta_{c}(x,\xi,\tau)$
$\displaystyle=\Theta(x_{c}-x)+\Theta(x_{c}+x-\xi)-\Theta(x_{c}-x)\Theta(x_{c}+x-\xi)\,,$
$\displaystyle\theta_{h}(x,\xi,\tau)$
$\displaystyle=\Theta(x-x_{c})\Theta(\xi-x_{c}-x)\,.$ (98)

The now $\tau$-dependent separation line $x_{c}(\tau)$ distinguishes the
different regimes.

The full solution can be written as

$D(x,\tau)=D_{0}(x,\tau)+\delta D(x,\tau)\,,$ (99)

where $D_{0}(x,\tau)$ is a solution to the coherent, soft kernel defined in
Eq. (93), and $\delta D(x,\tau)$ is a correcting factor.

Inserting this into the evolution equation (95), we get

$\displaystyle\frac{\partial}{\partial\tau}D(x,\tau)$
$\displaystyle=\int_{x}^{1}{\rm
d}\xi\,P(x,\xi,\tau)D_{0}(\xi,\tau)-D_{0}(x,\tau)\int_{0}^{x}{\rm
d}\xi\,P(\xi,x,\tau)$ $\displaystyle+\int_{x}^{1}{\rm
d}\xi\,P(x,\xi,\tau)\delta D(\xi,\tau)-\delta D(x,\tau)\int_{0}^{x}{\rm
d}\xi\,P(\xi,x,\tau)\,,$ (100)

Taking into account that $\theta_{c}+\theta_{h}=1$, one can rewrite the kernel
as

$\displaystyle P(x,\xi,\tau)$ $\displaystyle=P_{\rm coh}(x,\xi,\tau)+\delta
P(x,\xi,\tau)\,,$ (101)

where $\delta P\equiv(P_{\rm hard}-P_{\rm coh})\theta_{h}$. Inserting the new
kernel into Eq. (4.3), the term $\partial D_{0}/\partial\tau$ cancels, and we
are left with an iterative formula for $\delta D$,

$\displaystyle\delta D(x,\tau)$ $\displaystyle=D_{1}(x,\tau)$
$\displaystyle+\int^{\tau}_{0}{\rm d}\sigma\int_{x}^{1}{\rm
d}\xi\,P(x,\xi,\sigma)\delta D(\xi,\sigma)-\int^{\tau}_{0}{\rm d}\sigma\delta
D(x,\sigma)\int_{0}^{x}{\rm d}\xi\,P(\xi,x,\sigma)\,.$ (102)

Here we have defined the leading term in the correction as

$\displaystyle D_{1}(x,\tau)$ $\displaystyle=\int^{\tau}_{0}{\rm
d}\sigma\int_{x}^{1}{\rm d}\xi\,\delta
P(x,\xi,\sigma)D_{0}(\xi,\sigma)-\int^{\tau}_{0}{\rm d}\sigma
D_{0}(x,\sigma)\int_{0}^{x}{\rm d}\xi\,\delta P(\xi,x,\sigma)\,,$ (103)

which is given entirely in terms of known functions. To capture the main
modifications with respect to the purely coherent solution $D_{0}(x,\tau)$, it
is sufficient to keep only $D_{1}(x,\tau)$. This is sound because
$D_{0}(x,t=0)=\delta(1-x)$, implying that $\delta D(x,t=0)=0$. Therefore, all
terms going as $\sim\delta D$ start out small. At later times, single hard
emissions are rare, and thus $\delta D$ becomes less and less important.

Based on this discussion it is reasonable to assume that one can approximate
the intermediate time solution by the sum of the two leading terms
$D_{0}+D_{1}$. For this to be true it must be checked that it reproduces the
correct behavior at early and late times. The early times expansion is

$\displaystyle\lim_{\tau\to 0}(D_{0}(x,\tau)+D_{1}(x,\tau))$
$\displaystyle\simeq\int_{0}^{\tau}{\rm d}\sigma\int_{x}^{1}{\rm d}\xi\,P_{\rm
coh}(x,\xi,\sigma)\delta(1-\xi)$ $\displaystyle+\int_{0}^{\tau}{\rm
d}\sigma\int_{x}^{1}{\rm d}\xi\,\left[P_{\rm hard}(x,\xi,\tau)-P_{\rm
coh}(x,\xi,\tau)\right]\delta(1-\xi)$
$\displaystyle=\frac{\pi}{4{\rm\bar{\alpha}}}\frac{\tau^{2}}{x(1-x)^{2}}\,,$
(104)

where we kept the leading term and ignored virtual terms containing
$\delta(1-x)$. Therefore, the sum $D_{0}+D_{1}$ reproduces the hard part of
the early time expansion given in Eq. (91). Moreover, at late times
($t>L_{c}$), the phase space for hard ($\omega>\omega_{c}$) emissions
vanishes, and thus $D_{1}\to 0$. Hence, at late times the intermediate time
solution simply goes to the late time solution $D_{0}$. Consequently, we
expect that the sum of the two first terms $D_{0}+D_{1}$ to provide a decent
approximation of the true solution at all times. One can systematically
calculate corrections to this solution by iterating Eq. (4.3).

The first correction $D_{0}+D_{1}$ is shown in Fig. 9 with full lines for
different times. For short lengths (left panel), $D_{0}+D_{1}$ closely
resembles the early time solution, as expected. For late times $D_{0}+D_{1}$
reduces to $D_{0}$, as there is not any phase space for hard emissions left.
At intermediate times we see that $D_{0}+D_{1}$ goes to $D_{0}$ at low $x$,
while at high $x$ there is a suppression due to the lack of coherent
scatterings at early times. Qualitatively, Fig. 9 resembles the full numerical
solution shown in Fig. 10. In order to compare the two figures, we have marked
the value of $x_{\text{\tiny BH}}\equiv\omega_{\text{\tiny BH}}/E$ in Fig. 9.
We have also not included the color coding in this figure since it does not
include the physics from all the relevant regimes represented in Fig. 3.

##  5 Numerical evaluation of the medium cascade

In Sec. 3 we presented an effective framework that describes medium-induced
emissions up to arbitrary precision. By using this framework we showed in Sec.
4 how different scattering processes contribute to multiple induced emissions.
Based on the properties of the medium (e.g. length, mean free path), not all
induced emissions are necessary to resum (or to consider many of them). For
example, induced emissions from hard scatterings are not as important to resum
as emissions from multiple soft scatterings. We developed a simple analytic
model to include a single hard emission correction to the resummation of the
multiple soft ones. In Sec. 4, we made several simplifications e.g. we
simplified the emission kernel ${\cal K}$, and neglected Bethe–Heitler
emissions. It is important to understand how the energy distribution behaves
also without using these simplifications and to understand the error
introduced by employing them. That is the subject of this section.

In this section, we evaluate Eq. (85) numerically using the kernel in our new
framework. For more details about the numerical implementation, see App. D and
the complementary code MedKernels:2022 . The resulting energy distribution is
presented in Fig. 10, where the different lines stem from using emission
kernels ${\cal K}(z,E,t)$ at varying levels of approximation. The three panels
correspond to three different time stages in the evolution: $t<\lambda$,
$\lambda<t<L_{c}$, and $t>L_{c}$. As a reminder, there are no rare hard
splittings for $t>L_{c}$, which is evident on the rightmost panel of Fig. 10.

The evolution starts at $t=0$, with a single gluon of energy distribution
$D(x,0)=\delta(1-x)$ with $E=100$ GeV energy. We solely use gluons during the
evolution for simplicity, however our formalism is valid for other flavors
too. The kernels we use include finite-$z$ corrections, see App. C.3 for more
details. The vertical dashed lines in Fig. 10 separate the regions where
different scattering processes dominate, which are the same regions as in Fig.
3. The phase space for emissions is determined by comparing $z(1-z)E$ with
$\bar{\omega}_{c}(t),\omega_{\text{\tiny BH}}$, and $\omega_{c}(t)$ in the
relevant regions, as it is described in App. C. The turbulent cascade solution
$D_{0}(x,t)$ from Eq. (93) is also shown as a baseline, and was discussed in
more detail in Sec. 4.3.

Figure 10: The energy distribution solution of the medium-induced emission
evolution Eq. (85) at different times (different panels) and using different
scattering mechanisms to calculate the emission rate in Eq. (86) (different
line styles). The vertical dashed lines and the shaded areas separate
different scattering regions, where on the first panel $t<\lambda$ and on the
last panel $t>L_{c}$.

It is important to note that vacuum emissions are not included in the current
study, and thus many important effects (e.g. vacuum fragmentation, medium
resolution, color coherence), which are essential if one wishes to compare to
measurements, will not be discussed in this work. We refer the interested
readers to Refs. Mehtar-Tani:2012mfa ; Kurkela:2014tla ; Caucal:2018dla ;
Dominguez:2019ges ; Barata:2021byj  for further details.

#### Harmonic Oscillator:

A simple and much studied method of solving the evolution equation Eq. (85) is
by using the harmonic oscillator (HO) approximation of the emission rate (see
for example in Refs. Schenke:2009gb ; Blaizot:2013hx ; Blaizot:2015jea ;
Caucal:2018dla ; Caucal:2019uvr ; Caucal:2020xad ; Kutak:2018dim ;
Blanco:2020uzy ). The kernel is then simply given by

${\cal K}(z,E,t)\approx 2\left.\frac{{\rm d}^{2}I^{\text{\tiny HO}}}{{\rm
d}z{\rm d}t}\right|_{E}\,,$ (105)

where the HO approximation was discussed in Sec. 3.4 and in the numerics we
used Eq. (164). In this approximation, the induced emissions originate from
multiple soft elastic scatterings with the medium, and we expect this process
to dominate for $t\gg\lambda$ and $x_{\text{\tiny BH}}\ll x\ll x_{c}(t)$ (red
region in Fig. 3).

The energy distribution obtained by using the HO approximation is shown in
Fig. 10 with black dotted lines. The HO spectrum exhibits the well know
turbulent cascade behavior Blaizot:2013hx ; Blaizot:2015jea , resulting in a
characteristic tail $D(x\ll 1,t)\sim t/\sqrt{x}$, which can be seen on the
plot as horizontal lines at small $x$. This is also a feature of the
simplified analytic solution $D_{0}$. The turbulent cascade involves a
constant flux of energy propagating in time to $x\to 0$. The turbulence
appears below $x<x_{c}(t)$, and the energy increases with time due to the time
dependent emission phase space (see also $\omega_{c}(t)$ in Fig. 3).

Above $x>x_{c}(t)$ the HO kernel switches from ${\cal K}\sim 1/z^{3/2}$ to
$1/z^{3}$ and thus the energy distribution starts going as $D(x,t)\sim
1/x^{2}$. This is especially visible in the middle panel of Fig. 10.

It is interesting to note that due to the running of $\hat{q}$ (see Eq. (54)),
deviations from the pure $D\sim 1/\sqrt{x}$ are expected in the HO solution. A
consequence of this can be seen as a small kink on the dotted curves (most
visible in the left panel at $x\approx 10^{-4}$). Below the kink,
$\hat{q}(\omega)\to\hat{q}_{0}$ is used.

#### Improved Opacity Expansion:

The IOE from Sec. 3.4 makes it possible to extend the HO description to
include rare hard scatterings, covering both red and green regions in Fig. 3.
Here, the evolution equation is solved using the IOE kernel, which is

${\cal K}(z,E,t)\approx 2\left.\frac{{\rm d}^{2}I^{\rm IOE}}{{\rm d}z{\rm
d}t}\right|_{E}\,.$ (106)

This is valid from early to late times, for energies above the Bethe–Heitler
regime $x\gg x_{\text{\tiny BH}}$.

The energy distribution obtained using the IOE (HO+NHO is used in the
numerical implementation from Eqs. (164) - (171)) is shown in Fig. 10 with
black dashed lines. In the $x\ll x_{c}(t)$ region the HO dominates and
$D(x,t)$ qualitatively does not change by including hard emissions. The
turbulent tail $D\sim 1/\sqrt{x}$ is still present. The offset between the HO
and HO+NHO results originates from the effective contribution of the NHO term
to
$\hat{q}\approx\hat{q}_{0}\ln\frac{Q_{r}^{2}}{\mu_{\ast}^{2}}(1+1.016\ln^{-1}\frac{Q_{r}^{2}}{\mu_{\ast}^{2}})$.

In the region $x>x_{c}(t)$, the distribution function changes from the purely
HO result, due to the inclusion of the NHO corrections. The splitting function
here behaves as ${\cal K}\sim 1/z^{2}$, resulting in the distribution going as
$D(x)\sim 1/x$, which is visible in the green region. Our analytic result from
Sec. 4.3, given by $D_{0}+D_{1}$, includes a single hard emission in addition
to the HO cascade. Comparing the analytic results in Fig. 9 with the numerical
ones in Fig. 10 it is evident that the analytical result succeeds in capturing
qualitatively the behavior induced by using the IOE kernel.

Several models are based mainly on the $N=1$ term, see Refs. Idilbi:2008vm ;
Ovanesyan:2011xy ; Majumder:2013re , and also Refs. Guo:2000nz ; Wang:2013cia
. However, these models miss the soft scatterings that are present in a large
medium with $L>\lambda$. Soft scatterings are crucial to include to achieve an
accurate description of the energy distribution for $x<x_{c}(t)$. We would
like to also emphasize that by numerically solving Eq. (85), we also resum
hard emissions and thus we consider the possibility of emitting arbitrarily
many of them.

#### Full solution:

Finally, we present the energy distribution using our new picture. It is given
by solving the evolution equation (85) using the full kernel

${\cal K}(z,E,t)=2\frac{{\rm d}^{2}I^{\rm Full}}{{\rm d}z{\rm
d}t}=2\begin{cases}\frac{{\rm d}^{2}I^{\rm ROE}}{{\rm d}z{\rm d}t}\,,&{\rm
if}\quad z(1-z)E<\text{min}(\bar{\omega}_{c}(t),\omega_{\text{\tiny
BH}})\,,\\\ \frac{{\rm d}^{2}I^{\rm IOE}}{{\rm d}z{\rm d}t}\,,&{\rm
otherwise}\,.\end{cases}$ (107)

This covers the full emission phase space in ($z,t,E$) and was derived in
detail in Sec. 3. In practice we used $N_{r}=1$ and HO+NHO terms, derived from
Eqs. (42), (45), (69) and (72). This provides an excellent approximation to
the true result, as can be seen when comparing to the numerical evaluation of
the kernel in Fig. 7.

The resulting energy distribution is shown as full black lines in Fig. 10. One
can see that the results obtained by solely using the IOE (dashed lines)
agrees very well with the full kernel for $x\gg x_{\text{\tiny BH}}$, but they
start to differ when $x\lesssim x_{\text{\tiny BH}}$. It is clear that all the
difference between these two curves comes from Bethe–Heitler emissions. An
increasing tail appears at low $x$, which becomes more and more important for
later times, due to the logarithmic soft limit of the ROE kernel. In this
region the ROE kernel is ${\cal K}\sim\ln(z)/z$, resulting in the energy
distribution going as $D\sim\ln\frac{1}{x}$, which is similar to the DGLAP
evolution in vacuum.

A new, interesting bump also appears in the full solution close to
$x_{\text{\tiny BH}}$, which warrants an explanation. First, the energy flux
that brings quanta from $x=1$ to $x=0$ is not the same on the two sides of
$x_{\text{\tiny BH}}$. If the energy transport is more efficient from the
right (red region), that will result in a slowing of the flux when going from
higher to lower $x$. The bump can then be understood as a sediment of energy
building up around $x_{\text{\tiny BH}}$. Secondly, the emission kernel
defined in Eq. (107) is not smooth around the transition point
$\text{min}(\bar{\omega}_{c}(t),\omega_{\text{\tiny BH}})$. This discontinuity
introduces additional uncertainties in the true behavior of the energy
distribution around $x_{\text{\tiny BH}}$. This transition would smoothen by
including more orders of the expansion, which should be studied in the future.
This uncertainty could also be connected with the observed bump. In our
numerical implementation we introduced a smoothing function to minimize this
uncertainty, but it is still present. In addition, other effects like $2\to 2$
scattering process and thermal masses are also important here
Schlichting:2020lef , which we do not presently discuss. A thorough study of
this region is needed where all these effects are included, which motivates
future work.

Bethe–Heitler emissions have a soft divergence and thus an IR regulator
$\omega_{\text{min}}$ has to be introduced. A similar regulator was introduced
in the numerical differential equation solver. For more details, see App. D.

##  6 Conclusions and outlook

In this paper, we studied emissions induced by elastic scatterings on the
quark gluon plasma. We considered the interplay of the relevant length scales
of the problem, which are the propagation length, the mean free path between
scatterings, and the formation time. From these considerations, we derived the
emergent energy scales that separate the induced emission spectrum into
regimes governed by different scattering processes. The emerging hierarchy of
scales and related processes is illustrated in Fig. 3.

We presented a new theoretical framework consisting of different perturbative
expansions, namely the opacity expansion (OE), the improved OE (IOE), and the
resummed OE (ROE), which we derived rigorously. Together these are suitable to
describe the induced emission spectrum in different regimes. We showed that at
least one of the expansions is valid in every phase space ($\omega,t$) point,
where $\omega$ is the energy of the emitted gluon and $t$ is the propagation
time (length) in the medium.111111For the general case valid beyond the
strictly soft limit $\omega$ refers to the reduced energy of the three-body
evolution $\omega=z(1-z)E$, where $E$ is the initial energy of the emitter and
$z$ is the momentum sharing fraction. While these expansions are formally
simply reorganizations of the multiple scattering series, it is also important
to note that none of the expansions when truncated at any fixed order is valid
everywhere. Relying on the multiple approaches to obtain the spectrum, our
composite framework can account for an arbitrary number of medium interactions
that can be soft or hard. Most importantly, it is systematically improvable
and can, in principle, describe the emission spectrum (and rate) up to
arbitrary precision, which has only been achieved numerically before
CaronHuot:2010bp ; Feal:2018sml ; Andres:2020vxs .

In the current work, we have elucidated the convergence properties up to
second order in the studied resummations. Identifying the expansion structure
in the different regimes opens for the possibility of studying the accuracy of
resummations in the medium. Finally, our new description provides a quick and
efficient way to evaluate the induced emission spectrum (and rate), which is
an essential ingredient of medium-induced cascades and jet quenching study.
Our implementation is available in an online repository MedKernels:2022 .

In order to tackle multiple emissions, we also studied the in-medium energy
distribution $D(x,t)$ within our new formalism. This is an essential
ingredient of jet quenching phenomenology as it describes how the energy of a
leading particle gets distributed within a cascade. Having first identified
the conditions to generate multiple emissions in the medium, we demonstrated
the separation between early, rare emissions – generated mainly by a single
momentum exchange with the medium – and a following cascade of soft
splittings. We also showed that emissions are formed independently up to
suppressed terms, justifying a posteriori the formulation of the cascade via a
rate equation. We developed new analytic tools to combine the resummation of
multiple soft with rare hard emissions, and thus we showed how different
scattering processes appear on the level of the energy distribution. Finally,
using numerical evaluation based on the previously derived full splitting
kernels, we showed more rigorously the transition effects between different
scattering regions. We identified the importance of the time dependent phase
space separation, multiple Bethe–Heitler emissions, and the running of
$\hat{q}$, to mention some.

Even though this work can immediately be applied to jet quenching
phenomenology one should also tackle other challenges, such as vacuum
emissions and coherence effects, which we have neglected here and left for
future studies. We also acknowledge the importance of soft emissions for which
one should study thermal masses, energy changing $2\to 2$ scatterings, and
thermalization Ghiglieri:2015ala ; Schlichting:2020lef  to correctly describe
the very infrared regime close to the thermal scale and below. Moreover, other
non-perturbative effects such as expanding, inhomogeneous medium are also
important, in which direction our framework is extendable Casalderrey-
Solana:2012evi ; Adhya:2019qse ; Caucal:2020xad ; Adhya:2021kws .

###### Acknowledgements.

We thank Paul Caucal, Alexandre Falcão, Aleksas Mazeliauskas, and Sören
Schlichting for helpful discussions. We also thank Carlota Andres, Fabio
Dominguez, and Marcos Gonzalez Martinez for providing us the numerical results
for the spectrum computed in Ref. Andres:2020kfg . The work is supported by a
Starting Grant from Trond Mohn Foundation (BFS2018REK01) and the University of
Bergen.

##  Appendix A All order formulas for medium-induced spectrum in the soft
limit

Here we present some general formulas for the different expansions, that aided
us in calculating the spectrum in Sec. 3. In all cases the full spectrum is
given as the sum of all the terms $\omega\frac{{\rm d}I}{{\rm
d}\omega}=\sum_{n}\omega\frac{{\rm d}I^{N=n}}{{\rm d}\omega}$. The expansion
therefore only converges if the terms decrease sufficiently fast order by
order $\omega\frac{{\rm d}I^{N=n+1}}{{\rm d}\omega}<\omega\frac{{\rm
d}I^{N=n}}{{\rm d}\omega}$.

###  A.1 Opacity expansion

Combining the formulas for the opacity expansion Eq. (3.1) and the spectrum
Eq. (23) one can write the formula for the $n^{\rm th}$ term of the opacity
expansion as

$\displaystyle\omega\frac{{\rm d}I^{N=n}}{{\rm
d}\omega}=(-1)^{n-1}\frac{4\alpha_{s}C_{R}}{\omega}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\Sigma({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}v({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots
v({\bm{p}}_{2}-{\bm{p}}_{1})$ (108) $\displaystyle\times{\rm
Re}\,i\int_{0}^{L}{\rm d}t_{n}\int_{0}^{t_{n}}{\rm
d}t_{0}\int_{t_{0}}^{t_{n}}{\rm d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm
d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm d}t_{1}\,{\rm
e}^{-i\frac{{\bm{p}}_{n}^{2}}{2\omega}(t_{n}-t_{n-1})}\dots\rm
e^{-i\frac{{\bm{p}}_{1}^{2}}{2\omega}(t_{1}-t_{0})}\,,$

where $v({\bm{p}},s)$ is given in Eq. (18) and $\Sigma({\bm{p}}^{2})$ in given
in Eq. (24) in the GW model. It is useful to go to unitless integration
variables, by defining $\sqrt{L/2\omega}{\bm{p}}_{k}\to{\bm{p}}_{k}$ and
$\frac{t_{k}}{L}\to t_{k}$,

$\displaystyle\omega\frac{{\rm d}I^{N=n}}{{\rm
d}\omega}=(-1)^{n-1}8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{n}\frac{\bar{\omega}_{c}}{\omega}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{v}({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\tilde{v}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times{\rm Re}\,i\int_{0}^{1}{\rm d}t_{n}\int_{0}^{t_{n}}{\rm
d}t_{0}\int_{t_{0}}^{t_{n}}{\rm d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm
d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm d}t_{1}\,{\rm
e}^{-i{\bm{p}}_{n}^{2}(t_{n}-t_{n-1})}\dots\rm
e^{-i{\bm{p}}_{1}^{2}(t_{1}-t_{0})}\,.$ (109)

Here we have used the unitless function
$\tilde{v}({\bm{p}})=(2\pi)^{2}\delta({\bm{p}})-\frac{\bar{\omega}_{c}}{\omega}\tilde{\sigma}({\bm{p}})$,
where in the GW model

$\displaystyle\tilde{\sigma}({\bm{p}})$
$\displaystyle=\frac{4\pi}{\left({\bm{p}}^{2}+\frac{\bar{\omega}_{c}}{\omega}\right)^{2}}\,,$
$\displaystyle\tilde{\Sigma}({\bm{p}}^{2})$
$\displaystyle=\frac{1}{{\bm{p}}^{2}+\frac{\bar{\omega}_{c}}{\omega}}\,.$
(110)

The momentum and time integrals only depend on the unitless combination
$\bar{\omega}_{c}/\omega$. This means that the spectrum will take the form
$\omega\frac{{\rm d}I^{N=n}}{{\rm
d}\omega}={\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{n}h_{n}(\frac{\omega}{\bar{\omega}_{c}})$,
for some function $h_{n}$. This would naively imply convergence for
$L/\lambda<1$, but the exact for of the function $h_{n}$ must also be taken
into account. In Sec. 3.2 we showed for $N=1,2$ at low energy
$\omega\ll\bar{\omega}_{c}$ that $h_{n}(\frac{\omega}{\bar{\omega}_{c}})$ is a
finite function, meaning the series converges when $L/\lambda<1$. However, at
high energy $\omega\gg\bar{\omega}_{c}$ the function takes the form
$h_{n}(\frac{\bar{\omega}_{c}}{\omega})\sim(\frac{\bar{\omega}_{c}}{\omega})^{n}\tilde{h}_{n}(\frac{\bar{\omega}_{c}}{\omega})$,
where $\tilde{h}_{n}$ is finite, implying convergence when
$\frac{L}{\lambda}\frac{\bar{\omega}_{c}}{\omega}<1$. Proving this is true to
all orders is deferred to future work.

One can also derive a corresponding formula for the rate $\omega\frac{{\rm
d}I^{N=n}}{{\rm d}\omega{\rm d}L}$, by taking the length derivative of Eq.
(A.1)

$\displaystyle\omega\frac{{\rm d}I^{N=n}}{{\rm d}\omega{\rm
d}L}=(-1)^{n-1}8\pi{\rm\bar{\alpha}}\frac{1}{L}\left(\frac{L}{\lambda}\right)^{n}\frac{\bar{\omega}_{c}}{\omega}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{v}({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\tilde{v}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times{\rm Re}\,i\int_{0}^{1}{\rm d}t_{0}\int_{t_{0}}^{1}{\rm
d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm
d}t_{1}\,{\rm e}^{-i{\bm{p}}_{n}^{2}(1-t_{n-1})}{\rm
e}^{-i{\bm{p}}_{n-1}^{2}(t_{n-1}-t_{n-2})}\dots\rm
e^{-i{\bm{p}}_{1}^{2}(t_{1}-t_{0})}\,.$ (111)

###  A.2 Resummed opacity expansion

One can also derive all order formulas for the resummed opacity expansion, by
starting with Eq. (3.3) and following the same procedure,

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm
d}\omega}=\frac{4\alpha_{s}C_{R}}{\omega}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\Sigma({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\sigma({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\sigma({\bm{p}}_{2}-{\bm{p}}_{1})$
(112) $\displaystyle\times{\rm Re}\,i\int_{0}^{L}{\rm
d}t_{n}\int_{0}^{t_{n}}{\rm d}t_{0}\int_{t_{0}}^{t_{n}}{\rm
d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm
d}t_{1}\,\Delta(t_{n},t_{0}){\rm
e}^{-i\frac{{\bm{p}}_{n}^{2}}{2\omega}(t_{n}-t_{n-1})}\dots\rm
e^{-i\frac{{\bm{p}}_{1}^{2}}{2\omega}(t_{1}-t_{0})}\,.$

The ROE and OE expansions are equivalent at infinite order at low opacity
$L/\lambda\ll 1$, but at finite order the terms are mixed up. The terms
containing a delta function in the potential $v$ in the OE are included in the
Sudakov factor $\Delta$ in the ROE. One can see this by taking $\Delta\to 1$
and $v\to-\sigma$ in Eq. (112) and Eq. (108) respectively, in which case the
expansions become exactly the same. This is however not a good approximation
at any order.

After changing to unitless variables we arrive at

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm
d}\omega}=8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{n}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{n}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{\sigma}({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\tilde{\sigma}({\bm{p}}_{2}-{\bm{p}}_{1})$
(113) $\displaystyle\times{\rm Re}\,i\int_{0}^{1}{\rm
d}t_{n}\int_{0}^{t_{n}}{\rm d}t_{0}\int_{t_{0}}^{t_{n}}{\rm
d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm
d}t_{1}\,{\rm e}^{-\frac{L}{\lambda}(t_{n}-t_{0})}{\rm
e}^{-i{\bm{p}}_{n}^{2}(t_{n}-t_{n-1})}\dots\rm
e^{-i{\bm{p}}_{1}^{2}(t_{1}-t_{0})}\,,$

where we have used that in the static medium
$\Delta(t_{n},t_{0})=\exp(-\frac{t_{n}-t_{0}}{\lambda})$. Again we use the
unitless functions $\tilde{\Sigma}({\bm{p}}^{2})$ and
$\tilde{\sigma}({\bm{p}})$ defined in Eq. (A.1) for the GW model. In this case
the time integrals can be done analytically

$\displaystyle i\int_{0}^{1}{\rm d}t_{n}\int_{0}^{t_{n}}{\rm
d}t_{0}\int_{t_{0}}^{t_{n}}{\rm d}t_{n-1}\int_{t_{0}}^{t_{n-1}}{\rm
d}t_{n-2}\dots\int_{t_{0}}^{t_{2}}{\rm d}t_{1}\,{\rm
e}^{-\frac{L}{\lambda}(t_{n}-t_{0})}{\rm
e}^{-i{\bm{p}}_{n}^{2}(t_{n}-t_{n-1})}\dots\rm
e^{-i{\bm{p}}_{1}^{2}(t_{1}-t_{0})}$ $\displaystyle=i^{n}\int_{0}^{1}{\rm
d}t_{n}\int_{0}^{t_{n}}{\rm d}t_{0}\sum_{k=1}^{n}\frac{{\rm
e}^{-i\left({\bm{p}}_{k}^{2}-i\frac{L}{\lambda}\right)(t_{n}-t_{0})}}{\prod_{l\neq
k}^{n}({\bm{p}}_{k}^{2}-{\bm{p}}_{l}^{2})}=\left(\cos{\frac{n\pi}{2}}+i\sin{\frac{n\pi}{2}}\right)\sum_{k=1}^{n}\frac{T\left({\bm{p}}_{k}^{2}-i\frac{L}{\lambda}\right)}{\prod_{l\neq
k}^{n}({\bm{p}}_{k}^{2}-{\bm{p}}_{l}^{2})}\,,$ (114)

where we have defined the function

$\displaystyle T(x)$ $\displaystyle=\int_{0}^{1}{\rm
d}t_{n}\int_{0}^{t_{n}}{\rm d}t_{0}\,{\rm
e}^{-ix(t_{n}-t_{0})}=\frac{1-ix-{\rm e}^{-x}}{x^{2}}\,.$ (115)

The reason this simplification was not possible in the OE is that the formula
for the time integrals is only valid if ${\bm{p}}_{k}\neq{\bm{p}}_{l}$, while
the OE has terms containing $\delta({\bm{p}}_{k}-{\bm{p}}_{k-1})$. However, in
the ROE there are no such deltas, so the formula is valid. In the end the
spectrum is

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm d}\omega}$
$\displaystyle=8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{n}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{n}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{\sigma}({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\tilde{\sigma}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times\sum_{k=1}^{n}\frac{\cos{\frac{n\pi}{2}}{\rm
Re}\,T({\bm{p}}_{k}^{2}-i\chi)-\sin{\frac{n\pi}{2}}{\rm
Im}\,T({\bm{p}}_{k}^{2}-i\chi)}{\prod_{l\neq
k}^{n}({\bm{p}}_{k}^{2}-{\bm{p}}_{l}^{2})}\,.$ (116)

The remaining momentum integrals will in the end be some function of
$\bar{\omega}_{c}/\omega$ and opacity $\chi=L/\lambda$. As a reference, the
real and imaginary parts of $T({\bm{p}}^{2}-i\chi)$ are

$\displaystyle{\rm Re}\,T({\bm{p}}^{2}-i\chi)$
$\displaystyle=\frac{\chi({\bm{p}}^{4}+\chi^{2})+{\bm{p}}^{4}-\chi^{2}-{\rm
e}^{-\chi}\left(({\bm{p}}^{4}-\chi^{2})\cos{{\bm{p}}^{2}}+2\chi{\bm{p}}^{2}\sin{{\bm{p}}^{2}}\right)}{({\bm{p}}^{4}+\chi^{2})^{2}}$
$\displaystyle-{\rm Im}\,T({\bm{p}}^{2}-i\chi)$
$\displaystyle=\frac{{\bm{p}}^{2}({\bm{p}}^{4}+\chi^{2})-2\chi{\bm{p}}^{2}-{\rm
e}^{-\chi}\left(({\bm{p}}^{4}-\chi^{2})\sin{{\bm{p}}^{2}}-2\chi{\bm{p}}^{2}\cos{{\bm{p}}^{2}}\right)}{({\bm{p}}^{4}+\chi^{2})^{2}}\,.$
(117)

It is also possible to look at the limits of this function when the opacity
$\chi$ is low or high. In the low opacity case $\chi\ll 1$ we have

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm d}\omega}$
$\displaystyle\simeq
8\pi{\rm\bar{\alpha}}\left(\frac{L}{\lambda}\right)^{n}\left(\frac{\bar{\omega}_{c}}{\omega}\right)^{n}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{\sigma}({\bm{p}}_{n}-{\bm{p}}_{n-1})\dots\tilde{\sigma}({\bm{p}}_{2}-{\bm{p}}_{1})$
$\displaystyle\times\sum_{k=1}^{n}\frac{\cos{\frac{n\pi}{2}}(1-\cos{{\bm{p}}_{k}^{2}})+\sin{\frac{n\pi}{2}}({\bm{p}}_{k}^{2}-\sin{{\bm{p}}_{k}^{2}})}{{\bm{p}}_{k}^{4}\prod_{l\neq
k}^{n}({\bm{p}}_{k}^{2}-{\bm{p}}_{l}^{2})}\,.$ (118)

In this case the momentum integrals will give a function of
$\bar{\omega}_{c}/\omega$. As mentioned in Sec. 3.3 is it preferable to use
the OE (A.1) at low opacity, as the order of opacity is mixed up in the ROE.

In the high opacity case $\chi\gg 1$ you get

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm d}\omega}$
$\displaystyle\simeq
8\pi{\rm\bar{\alpha}}\frac{L}{\lambda}\left(\frac{\omega_{\text{\tiny
BH}}}{\omega}\right)^{n}\int_{{\bm{p}}_{n},\dots,{\bm{p}}_{1}}\,\tilde{\Sigma}({\bm{p}}_{n}^{2},\omega_{\text{\tiny
BH}})\frac{{\bm{p}}_{n}\cdot{\bm{p}}_{1}}{{\bm{p}}_{n}^{2}}\tilde{\sigma}({\bm{p}}_{n}-{\bm{p}}_{n-1},\omega_{\text{\tiny
BH}})\dots\tilde{\sigma}({\bm{p}}_{2}-{\bm{p}}_{1},\omega_{\text{\tiny BH}})$
$\displaystyle\times\sum_{k=1}^{n}\frac{\cos{\frac{n\pi}{2}}+\sin{\frac{n\pi}{2}}{\bm{p}}_{k}^{2}}{(1+{\bm{p}}_{k}^{4})\prod_{l\neq
k}^{n}({\bm{p}}_{k}^{2}-{\bm{p}}_{l}^{2})}\,.$ (119)

Here we see the emergence of the scale $\omega_{\text{\tiny
BH}}=\frac{1}{2}\mu^{2}\lambda$, which takes the place of $\bar{\omega}_{c}$
in the functions $\tilde{\Sigma}$ and $\tilde{\sigma}$, as emphasized in the
above formula. The momentum integrals then only become some function of
$\omega_{\text{\tiny BH}}/\omega$. Notice that this scales as $\sim L/\lambda$
at every order, meaning it can converge also when $L/\lambda>1$.

In Sec. 3.4 we calculated the high opacity limit to orders $N=1,2$. There it
was clear that for low $\omega\ll\omega_{\text{\tiny BH}}$ the $N=2$ limit is
subleading compared to the $N=1$ limit, while at high
$\omega\gg\omega_{\text{\tiny BH}}$ the $N=2$ and $N=1$ limits are of the same
order. This seems to imply that the expansion converges at low $\omega$, while
it breaks down at high $\omega$. Again, proving this for all orders is
deferred to future work.

Again the calculation for the rate $\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm
d}\omega{\rm d}L}$ is similar

$\displaystyle\omega\frac{{\rm d}I^{N_{r}=n}}{{\rm d}\omega{\rm