# A no expanding breather theorem for noncompact Ricci flows

Liang Cheng and Yongjia Zhang  School of Mathematics and Statistics $\&$ Hubei
Key Laboratory of Mathematical Sciences, Central China Normal University,
Wuhan, 430079, P.R.China<EMAIL_ADDRESS>School of Mathematics,
University of Minnesota, Twin Cities, MN, 55414, USA<EMAIL_ADDRESS>
###### Abstract.

In this note we show that, under certain curvature positivity conditions (the
weak $\operatorname{PIC}-2$ condition or the nonnegative bisectional curvature
condition), a complete and noncompact expanding breather of the Ricci flow is
also an expanding gradient Ricci soliton. This is the first no expanding
breather theorem in the noncompact category.

###### Key words and phrases:

expanding breathers, gradient expanding Ricci solitons, no breather theorem,
noncompact Ricci flows

Liang Cheng’s Research partially supported by Natural Science Foundation of
Hubei 2019CFB511

##  1\. Introduction

Perelman [13] proved the no shrinking, steady, and expanding breather theorems
on compact manifolds, by applying the monotonicity formulas of his
$\mathcal{W}$-functional, $\mathcal{F}$-functional, and normalized
$\mathcal{F}$-functional, respectively: these types of breathers must also be
gradient Ricci solitons of the corresponding types, which evolve by self-
diffeomorphism and scaling. In fact, in the compact category, the steady and
expanding gradient Ricci solitons must also be Einstein manifolds. Let us
first of all recall the definitions of breathers and solitons.

######  Definition 1.1.

Let $(M,g(t))$ be a complete Ricci flow. If there exist two time instances
$t_{1}<t_{2}$, a constant $\alpha>0$, and a self-diffeomorphism
$\phi:M\rightarrow M$, such that

$\displaystyle g(t_{1})=\alpha\phi^{*}g(t_{2}),$

then $(M,g(t))$ is called a breather. If $\alpha=1$, $\alpha>1$, or
$\alpha<1$, then the breather is called steady, shrinking, or expanding,
respectively.

######  Definition 1.2.

A Ricci soliton is a tuple $(M,g,X,\lambda)$, where $(M,g)$ is a complete
Riemannian manifold, $X$ is a smooth vector field on $M$, and $\lambda>0$ is a
constant, satisfying

$\displaystyle Rc+\frac{1}{2}\mathcal{L}_{X}g=\frac{\lambda}{2}g.$

If $\lambda=0$, $\lambda>0$, or $\lambda<0$, then the soliton is called
steady, shrinking, or expanding, respectively. The soliton is called
_complete_ if the vector field $X$ is complete. If there exists a smooth
function $f$ on $M$ such that $X=\nabla f$, then $(M,g,f,\lambda)$ is called a
gradient Ricci soliton.

A complete Ricci soliton $(M,g,X,\lambda)$ always generates a canonical form,
that is, a Ricci flow $g(t)=\tau(t)\phi_{t}^{*}g$ which moves by self-
diffeomorphism and scaling, where $\tau(t)=\lambda t+1$ and
$\frac{\partial}{\partial t}\phi_{t}(x)=\frac{1}{\tau(t)}X(\phi_{t}(x))$. As
indicated by Perelman [13], if one views a Ricci flow as an orbit in the space
$\operatorname{Met}(M)/\operatorname{Diff}(M)$, where $\operatorname{Met}(M)$
is the space of all smooth Riemannian metrics and $\operatorname{Diff}(M)$
stands for the group of self-diffeomorphisms and scalings, then a breather is
a periodic orbit and a soliton is a static one. Therefore, the no breather
theorem is tantamount to saying that the periodic orbits must also be static.

Perelman’s proofs of the no breather theorems require the existence of
minimizers for the $\mathcal{W}$-functional, the $\mathcal{F}$-functional, and
the normalized $\mathcal{F}$-functional, respectively. When the manifold is
compact, such existence results are straightforward applications of
variational problems. A natural question to ask is, under what condition can
the no breather theorems be established for noncompact manifolds, and how to
carry our the proof. One may certainly attempt to find the minimizers of these
functionals, and this is possible under certain geometric conditions. For
results of this type, the reader may refer to [14] and [15]. We remark here
that, in the noncomapct category, this approach is almost impossible for the
no steady or expanding breather theorems. The reason is that the
$\mathcal{F}$-functional or the normalized $\mathcal{F}$-functional are
generally not finite on noncompact steady or expanding gradient Ricci
solitons, respectively, when these functionals are evaluated using the
potential functions of the corresponding solitons. (One may think of a Bryant
soliton for example.)

Another approach was initiated by the result of Lu-Zheng [12], where they
constructed a Type I ancient solution using a shrinking breather, and proved,
under certain geometric conditions, that the backward blow-down limit of this
ancient solution must be the breather itself, which must also be a shrinking
gradient Ricci soliton by Naber [9]. This method was refined, and the
conclusion improved, by the second author [16], where the condition for the no
shrinking breather theorem is reduced to bounded curvature alone. The authors
[4] recently further reduced the condition for the no shrinking breather
theorem to a lower bound of the Ricci curvature alone.

In this paper, we continue our study in [4] and apply our method to noncompact
expanding breathers. Recall that Feldman-Ilmanen-Ni [6] established some
forward monotonicity formulas for the Ricci flow as the dual version of
Perelman’s [13], whose equalities are fulfilled on expanding gradient Ricci
solitons. We will be implementing Feldman-Ilmanen-Ni’s forward reduced
geometry in this paper in the same way as we have applied Perelman’s reduced
geometry in [4]. However, since the forward reduced geometry does not behave
as nicely as Perelman’s reduced geometry (The reason, intuitively speaking, is
this, that on steady or expanding solitons, as it is in the case of a
shrinking soliton, the forward reduced volumes should coincide with the
$\mathcal{F}$-functional or the normalized $\mathcal{F}$-functional,
respectively, evaluated using the corresponding potential functions, whereas
the latter two are generally infinite in the noncompact case), we will need to
impose some strong curvature conditions.

######  Theorem 1.3.

Let $(M,g(t))$ be a complete and noncompact expanding breather of the Ricci
flow. Assume $g(t)$ has bounded curvature on each time-slice and either one of
the following is true.

  1. (1)

$g(t)$ satisfies the weak $\operatorname{PIC}-2$ condition.

  2. (2)

$(M,g(t))$ is Kähler with nonnegative bisectional curvature.

Then $(M,g(t))$ is the canonical form of an expanding gradient Ricci soliton.

######  Remark 1.4.

Our proof depends heavily on Hamilton’s Harnack estimate ([1], [2], and [11]),
and this is why we assumed the above curvature conditions. The reader may
easily verify that, if Hamilton’s trace Harnack is assumed to be valid, then a
nonnegative Ricci curvature assumption is sufficient for the proof. Lott [8]
gave an example of a complete but nongradient expanding soliton on noncompact
manifold (see page 635 in therein). This soliton has bounded curvature but
does not have nonnegative Ricci curvature. Therefore, Theorem 1.3 cannot be
proved in general without any curvature positivity condition.

Before we conclude our introduction, a word is to be said about our method.
Similar to [4], we constructed a Type III immortal solution starting from the
given expanding breather and considered the monotonicity of the quantity

$\displaystyle\tilde{\theta}_{+}(t)=\frac{1}{(4\pi
t)^{\frac{n}{2}}}\int_{M}e^{-\ell_{+}}d\mu_{t},$

where $\ell_{+}$ is the forward reduced distance constructed in [6]. Note that
the monotonicity of this quantity relies on Hamilton’s Harnack, as proved in
Theorem 4.3 of [10]. However, the forward reduced volume defined in [6]

$\displaystyle\theta_{+}(t)=\frac{1}{(4\pi
t)^{\frac{n}{2}}}\int_{M}e^{\ell_{+}}d\mu_{t},$

being automatically decreasing along all Ricci flows on compact manifolds, is
yet generally infinite on noncompact manifolds.

The organization of the paper is as follows. In section 2, we review basic
$\mathcal{L}_{+}$-geometry introduced by Feldman-Ilmanen-Ni [6]. In section 3,
we estimate the $\ell_{+}$-distance on Type III immortal Ricci flows. In
section 4, we prove the main theorem.

##  2\. Preliminaries

Let $g(t)$ be a metric evolving by the Ricci flow equation on $M\times[0,T]$.
We always assume that either $M$ is compact or $g(t)$ has bounded curvature at
each time. Fixing a point $x$, Feldman-Ilmanen-Ni [6] defined the following
dual version of Perelman’s reduced distance, called the _forward reduced
distance function_.

$\displaystyle\ell_{+}(y,t)=\frac{1}{2\sqrt{t}}\inf\limits_{\gamma}\int^{t}_{0}\sqrt{\eta}\left(R(\gamma(\eta),\eta)+|\gamma^{\prime}(\eta)|_{g(\eta)}^{2}\right)d\eta,$
(2.1)

where $(y,t)\in M\times(0,T]$, and the infimum is taken among all piecewise
smooth curves $\gamma:[0,t]\rightarrow M$ satisfying $\gamma(0)=x$ and
$\gamma(t)=y$. $(x,0)$ is called the _base point_ of $\ell_{+}$. We remark
that since the minimizing curve in (2.1) also satisfies an
$\mathcal{L}_{+}$-geodesic equation, whose form is very similar to that of an
$\mathcal{L}$-geodesic equation, one may easily modify the arguments in, say,
Chapter 7 in [5], to verify that $\ell_{+}$ is locally Lipschitz under our
assumption. We then summarize some equations and inequalities satisfied by
$\ell_{+}$. The reader may note their similarity to the case of Perelman’s
$\mathcal{L}$-geometry.

######  Theorem 2.1 (Corollary 2.1 in [6]).

The $\ell_{+}$ function satisfies the following equalities for almost every
$(y,t)\in M\times(0,T]$

$\displaystyle\frac{\partial\ell_{+}}{\partial
t}=R-\frac{\ell_{+}}{t}-\frac{K}{2t^{\frac{3}{2}}},$ (2.2)
$\displaystyle|\nabla\ell_{+}|^{2}=\frac{\ell_{+}}{t}-R+\frac{K}{t^{\frac{3}{2}}},$
(2.3)

Moreover, $\ell_{+}$ satisfies the following inequalities in the barrier sense
or in the sense of distribution.

$\displaystyle\Delta\ell_{+}\leq R+\frac{n}{2t}-\frac{K}{2t^{\frac{3}{2}}},$
(2.4) $\displaystyle\frac{\partial\ell_{+}}{\partial
t}+\Delta\ell_{+}+|\nabla\ell_{+}|^{2}-R-\frac{n}{2t}\leq 0,$ (2.5)
$\displaystyle 2\Delta\ell_{+}+|\nabla\ell_{+}|^{2}-R-\frac{\ell_{+}+n}{t}\leq
0,$ (2.6)

where

$\displaystyle K=\int^{t}_{0}\eta^{\frac{3}{2}}H(X)d\eta,$ (2.7)

$X$ is the velocity of the minimizing $\mathcal{L}_{+}$-geodesic connecting
$(x,0)$ and $(y,t)$, and

$H(X)=\frac{\partial R}{\partial t}+2<\nabla R,X>+2Rc(X,X)+\frac{R}{t}$

is Hamilton’s trace Harnack. Furthermore, $\nabla\ell_{+}(y,t)=X(t)$ whenever
the minimizing $\mathcal{L}_{+}$-geodesic connecting $(x,0)$ and $(y,t)$ is
unique.

######  Theorem 2.2 (Theorem 4.3 in [10]).

Let $(M,g(t))$ be a Ricci flow with bounded curvature at each time slice.
Assume Hamilton’s trace Harnack is nonnegative, then the quantity

$\displaystyle\tilde{\theta}_{+}(t)=\frac{1}{(4\pi
t)^{\frac{n}{2}}}\int_{M}e^{-\ell_{+}}d\mu_{t}$

is monotonically decreasing in $t$, where $\mu_{t}$ is the Riemannian measure
of $g(t)$ and $\ell_{+}$ is the forward reduced distance based at some fixed
point on $M\times\\{0\\}$.

###### Sketch of proof.

Combing (2.2), (2.3), and (2.4), we have

$\displaystyle\frac{\partial\ell_{+}}{\partial
t}-\Delta\ell_{+}+|\nabla\ell_{+}|^{2}+R+\frac{n}{2t}\geq\frac{K}{t^{\frac{3}{2}}}\geq
0$ (2.8)

in the barrier sense or in the sense of distribution. Then, taking for granted
the integration by parts at infinity, we compute

$\displaystyle\frac{d}{dt}\int_{M}(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}(x,t)}d\mu_{t}$ $\displaystyle=$
$\displaystyle-\int_{M}\left(\frac{\partial\ell_{+}}{\partial
t}-\Delta\ell_{+}+|\nabla\ell_{+}|^{2}+R+\frac{n}{2t}\right)(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}(x,t)}d\mu_{t}$ $\displaystyle=$
$\displaystyle-\int_{M}\frac{K}{t^{\frac{3}{2}}}(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}(x,t)}d\mu_{t}\leq 0.$

Note that the above computation is valid in our case according to the
estimates in the next section. ∎

##  3\. Estimates for $\ell_{+}$ function on Type III Ricci flows

In this section, we use similar methods as in [9] to derive some estimates for
the forward reduced distance on a Type III immortal Ricci flow. Since Theorem
3.1(2)—(4) are already well established in [9], we shall be relatively brief
in their proofs. Recall that an immortal solution $(M,g(t))_{t\in[0,\infty)}$
is called Type III, if there exists a constant $C_{0}>0$, such that

$\displaystyle|Rm|\leq\frac{C_{0}}{t}$ everywhere on $\displaystyle
M\times(0,\infty).$ (3.1)

######  Theorem 3.1.

Let $(M^{n},g(t))_{t\in[0,\infty)}$ be an $n$-dimensional Type III immortal
Ricci flow with nonnegative Ricci curvature everywhere. Let $\ell_{+}$ be the
forward reduced distance function based a fixed point $(x,0)$. Furthermore,
assume that there exists a sequence of space-time points
$\\{(x_{j},t_{j})\\}_{j=1}^{\infty}$ such that $t_{j}\nearrow\infty$ and

$\ell_{+}(x_{j},t_{j})\leq A<\infty\ \text{ for all }\ j\geq 1.$ (3.2)

Then, there exists a positive constant $Q$ depending only on $A$,
$\alpha\in(0,1)$, $n$, and the Type III constant $C_{0}$ in (3.1), such that
the following inequalities hold for all $(y,t)\in
M^{n}\times[\alpha^{-1},\alpha]$ and for all $j\geq 1$, understood in the
barrier sense if any differentiation is involved.

  1. (1)

$\displaystyle|K^{j}|(y,t)\leq
Q\sqrt{t}\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}$,
$\displaystyle|\nabla K^{j}|(y,t)\leq
Q\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}$,
$\displaystyle\left|\frac{\partial}{\partial
t}K^{j}\right|(y,t)\leq\frac{Q}{\sqrt{t}}\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}$,

  2. (2)

$\displaystyle\frac{1}{Q}\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}-Q\leq\ell^{j}_{+}(y,t)\leq
Q\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}$,

  3. (3)

$\displaystyle|\nabla\ell^{j}_{+}|(y,t)\leq\frac{Q}{\sqrt{t}}\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)$
,

  4. (4)

$\displaystyle\left|\frac{\partial\ell^{j}_{+}}{\partial
t}\right|(y,t)\leq\frac{Q}{t}\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}$,

where $g_{j}(t):=t_{j}^{-1}g(tt_{j})$ is the Ricci flow sequence obtained by
Type III scaling, $K^{j}$ is as defined in (2.7) for the scaled Ricci flow
$g_{j}$, and $\ell_{+}^{j}(\cdot,t):=\ell_{+}(\cdot,tt_{j})$ is the forward
reduced distance based at $(x,0)$ and with respect to the Ricci flow
$g_{j}(t)$.

###### Proof.

In this proof, we will use the capital letter $C$ to denote a general
estimation constant, which depends on $\alpha$, $n$, $A$, and $C_{0}$ as
indicated in the statement of this theorem, and could vary from line to line.
Since the Ricci flow is Type III, by Shi’s gradient estimates we have

$\displaystyle|R|(y,t)\leq\frac{C}{t},\ \ \ |\nabla
Rm|(y,t)\leq\frac{C}{t^{\frac{3}{2}}},\ \ \ \left|\frac{\partial R}{\partial
t}\right|(y,t)\leq\frac{C}{t^{2}},\ \ \ \left|\nabla\frac{\partial}{\partial
t}R\right|\leq\frac{C}{t^{\frac{5}{2}}},$

for all $(y,t)\in M\times(0,\infty)$. Let $(y,t)\in M\times(0,\infty)$ be such
that the minimizing $\mathcal{L}_{+}$-geodesic $\gamma$ connecting $(x,0)$ and
$(y,t)$ is unique. We denote $X:=\gamma^{\prime}$ and calculate that

$\displaystyle|K|$
$\displaystyle=\left|\int^{t}_{0}\eta^{\frac{3}{2}}\left(\frac{\partial
R}{\partial\eta}+\frac{R}{\eta}+2\langle\nabla
R,X\rangle+2Rc(X,X)\right)d\eta\right|\leq\int^{t}_{0}\eta^{\frac{3}{2}}\left(\frac{C}{\eta^{2}}+\frac{C}{\eta^{\frac{3}{2}}}|X|+\frac{C}{\eta}|X|^{2}\right)d\eta$
$\displaystyle\leq\int_{0}^{t}\frac{C}{\eta^{\frac{1}{2}}}d\eta+\int_{0}^{t}\sqrt{\eta}|X|^{2}d\eta\leq
2C\sqrt{t}+C\int_{0}^{t}\sqrt{\eta}(|X|^{2}+R)d\eta=2C\sqrt{t}(1+\ell_{+}(y,t)).$
(3.3)

Since the Type III scaling process does not alter the Type III constant in
(3.1), we obtain from (2.2), (2.3), and (3) that

$\displaystyle\left|\frac{\partial{\ell^{j}_{+}}}{\partial t}\right|(y,t)$
$\displaystyle\leq\frac{C}{t}(1+\ell^{j}_{+}(y,t)),$ (3.4)
$\displaystyle|\nabla\ell^{j}_{+}|^{2}(y,t)$
$\displaystyle\leq\frac{C}{t}(1+\ell^{j}_{+}(y,t)).$ (3.5)

In view of the fact

$\displaystyle\ell^{j}_{+}(x_{j},1)=\ell_{+}(x_{j},t_{j})\leq A,$

we may integrate (3.4) to obtain

$\displaystyle\ell^{j}_{+}(x_{j},t)\leq C\ \ \text{ for all }\ \
t\in[\alpha^{-1},\alpha].$ (3.6)

Integrating (3.5) in space and applying (3.6), we have

$\displaystyle\ell^{j}_{+}(y,t)\leq
Q\left(1+\frac{d_{g_{j}(t)}(x_{j},y)}{\sqrt{t}}\right)^{2}\ \ \text{ for all
}\ \ (y,t)\in M\times[\alpha,\alpha^{-1}].$ (3.7)

Conclusion (3), (4), the second inequality of conclusion (2), and the first
inequality of conclusion (1) now follow from (3), (3.4), (3.5), and (3.7).

To obtain the first inequality of conclusion (2), we fix $(y,t)\in
M\times[\alpha,\alpha^{-1}]$ and let $\gamma_{1}(s)$ and $\gamma_{2}(s)$ be
minimizing $\mathcal{L}_{+}$-geodesics from $(x,0)$ to $(x^{j},t)$ and to
$(y,t)$, respectively, all with respect to $g_{j}(t)$. We denote
$f(s):=d_{g_{j}(s)}(\gamma_{1}(s),\gamma_{2}(s))$. Then we have

$\displaystyle\frac{d^{-}}{ds}f(s)$ $\displaystyle=\langle\nabla
d_{g_{j}(s)},\gamma_{1}^{\prime}(s)\rangle+\langle\nabla
d_{g_{j}(s)},\gamma_{2}^{\prime}(s)\rangle+\left(\frac{d^{-}}{d\tau}d_{g_{j}(\tau)}(\gamma_{1}(s),\gamma_{2}(s))\right)\Bigg{|}_{\tau=s}$
$\displaystyle=\langle\nabla
d_{g_{j}(s)},\nabla\ell^{j}_{+}(\gamma_{1}(s),s)\rangle+\langle\nabla
d_{g_{j}(s)},\nabla\ell^{j}_{+}(\gamma_{2}(s),s)\rangle+\left(\frac{d^{-}}{d\tau}d_{g_{j}(\tau)}(\gamma_{1}(s),\gamma_{2}(s))\right)\Bigg{|}_{\tau=s}$
$\displaystyle\leq|\nabla\ell^{j}_{+}(\gamma_{1}(s),s)|+|\nabla\ell^{j}_{+}(\gamma_{2}(s),s)|+\left(\frac{d^{-}}{d\tau}\int_{\sigma}\sqrt{g_{j}(\tau)(\sigma^{\prime},\sigma^{\prime})}\right)\Bigg{|}_{\tau=s},$

where $\sigma$ is a unit speed minimizing geodesic from $\gamma_{1}(s)$ to
$\gamma_{2}(s)$ with respect to metric $g_{j}(s)$, and
$\frac{d^{-}f}{ds}=\liminf\limits_{h\rightarrow 0^{+}}\frac{f(s+h)-f(s)}{h}$
is the lower forward Dini derivative. By Lemma 18.1 in [5], we have

$\displaystyle\left(\frac{d^{-}}{d\tau}\int_{\sigma}\sqrt{g_{j}(\tau)(\sigma^{\prime},\sigma^{\prime})}\right)\Bigg{|}_{\tau=s}=-\min\limits_{\eta\in\mathcal{Z}\left(s\right)}\int_{\eta}Rc_{g_{j}(s)}(\eta^{\prime},\eta^{\prime})\leq
0,$

where $\mathcal{Z}\left(s\right)$ denotes the set of all unit speed minimizing
geodesics from $\gamma_{1}(s)$ to $\gamma_{2}(s)$ with respect to metric
$g_{j}(s)$. Then, by (3.5), we get

$\displaystyle\frac{d^{-}}{ds}f(s)\leq\sqrt{\frac{C_{2}}{s}(1+\ell^{j}_{+}(\gamma_{1}(s),s))}+\sqrt{\frac{C_{2}}{s}(1+\ell^{j}_{+}(\gamma_{2}(s),s))}.$
(3.8)

On the other hand, since for all $s\in(0,t]$, we have

$\displaystyle\ell^{j}_{+}(\gamma_{2}(s),s)$
$\displaystyle=\frac{1}{2\sqrt{s}}\int^{s}_{0}\sqrt{\eta}(R+|X|^{2})d\eta\leq\frac{1}{2\sqrt{s}}\int^{t}_{0}\sqrt{\eta}(R+|X|^{2})d\eta=\frac{\sqrt{t}}{\sqrt{s}}\ell^{j}_{+}(y,t),$
$\displaystyle\ell^{j}_{+}(\gamma_{1}(s),s)$
$\displaystyle\leq\frac{\sqrt{t}}{\sqrt{s}}\ell^{j}_{+}(x_{j},t)\leq
C\frac{\sqrt{t}}{\sqrt{s}},$

where in the latter formula we have applied (3.6). Then, (3.8) becomes

$\displaystyle\frac{d^{-}}{ds}f(s)\leq
C\frac{t^{\frac{1}{4}}}{s^{\frac{3}{4}}}\left(1+\sqrt{1+l^{j}_{+}(y,t)}\right)\
\ \text{ for all }\ \ s\in(0,t].$

Integrating the above inequality from $0$ to $t$, we obtain the first
inequality of conclusion (2).

Finally, to obtain the last two inequalities of conclusion (1), we let
$\gamma$ be a minimizing $\mathcal{L}_{+}$-geodesic with respect to $g$, which
connects $(x,0)$ and $(y,t)$, and $Y$ an $\mathcal{L}_{+}$-Jacobi field along
$\gamma$, satisfying $[X,Y]=0$ and

$\displaystyle\nabla_{X}Y=Ric(Y)+\frac{1}{2\eta}Y,\ \ \
|Y(\eta)|^{2}=\frac{\eta}{t}|Y(t)|^{2}=\frac{\eta}{t},\ \ \
|\nabla_{X}Y|\leq\frac{C}{t^{\frac{1}{2}}\eta^{\frac{1}{2}}}.$ (3.9)

Then, we may compute

$\displaystyle|\delta_{Y}K|$
$\displaystyle=\Bigg{|}\int_{0}^{t}\eta^{\frac{3}{2}}\bigg{(}\Big{\langle}\nabla\frac{\partial}{\partial\eta}R,Y\Big{\rangle}+\frac{1}{\eta}\langle\nabla
R,Y\rangle+2\langle\nabla^{2}R,X\otimes Y\rangle+2\langle\nabla
R,\nabla_{X}Y\rangle$ (3.10) $\displaystyle\quad+2\nabla
Ric(Y,X,X)+2Ric(X,\nabla_{X}Y)\bigg{)}d\eta\Bigg{|}\leq
C\int_{0}^{t}\eta^{\frac{3}{2}}\left(\frac{1}{t^{\frac{1}{2}}\eta^{2}}+\frac{1}{t^{\frac{1}{2}}\eta^{\frac{3}{2}}}|X|+\frac{1}{t^{\frac{1}{2}}\eta}|X|^{2}\right)d\eta$
$\displaystyle\leq\frac{C}{\sqrt{t}}\int_{0}^{t}\frac{1}{\eta^{\frac{1}{2}}}d\eta+\frac{C}{\sqrt{t}}\int_{0}^{t}\sqrt{\eta}(|X|^{2}+R)d\eta\leq
C(1+\ell_{+}(y,t)).$

Since the Type III constant in (3.1) is not affected by the Type III scaling,
we obtain the second inequality of conclusion (1) by (3.10). Next, we observe
that

$\displaystyle\left|\frac{\partial}{\partial t}K\right|(y,t)$
$\displaystyle=\left|\frac{d}{d\eta}\Big{|}_{\eta=t}K-\langle\nabla
K,X\rangle(y,t)\right|\leq t^{\frac{3}{2}}\left|\frac{\partial R}{\partial
t}+\frac{R}{t}+2\langle\nabla R,X\rangle+2Rc(X,X)\right|+|\nabla K||X|$
$\displaystyle\leq\frac{C}{\sqrt{t}}+C|X|+C\sqrt{t}|X|^{2}+C(1+\ell_{+}(y,t))|X|\leq\frac{C}{\sqrt{t}}(1+\ell_{+}(y,t))+C\sqrt{t}|\nabla\ell_{+}|^{2}(y,t)$
$\displaystyle\leq\frac{C}{\sqrt{t}}(1+\ell_{+}(y,t)).$

Here we have used $X(t)=\nabla\ell_{+}(y,t)$ and formula (3.5). Again, by the
scaling invariance of the Type III constant, we obtain the third inequality of
conclusion (1).

∎

##  4\. the proof of the main theorem

In this section, we prove the main theorem by implementing the method of, [4],
[12], and [16]. Since our techniques and arguments are very similar to that of
[4], we are not including obvious modifications, and the readers are referred
to this paper for more detailed treatment.

Let $(M,g_{0}(t))$ be an expanding breather as described in the statement of
Theorem 1.3. Note that the assumptions therein guarantees the validity of
Hamilton’s Harnack estimate ([1], [2] and [11]). After rescaling and
translating in time, we may assume that there exists $\alpha\in(0,1)$ and a
diffeomorphism $\phi:M\to M$, such that

$\alpha g_{0}(1)=\phi^{*}g_{0}(0).$ (4.1)

Similar to [4], for each $j\geq 0$ we define

$\displaystyle\displaystyle t_{j}$ $\displaystyle=$
$\displaystyle\sum^{j}_{k=0}\alpha^{-k},\ \ t_{0}=1,$ $\displaystyle g_{j}(t)$
$\displaystyle=$
$\displaystyle\alpha^{-j}(\phi^{j})^{*}g_{0}(\alpha^{j}(t-t_{j-1})),\ \
t\in[t_{j-1},t_{j}].$

Obviously, $t_{j}\nearrow\infty$. We may then define the spliced immortal
solution

$\displaystyle g(t)=\left\\{\begin{array}[]{ll}g_{0}(t),&t\in[0,1],\\\
g_{j}(t),&t\in[t_{j-1},t_{j}].\end{array}\right.$

It is straightforward to check that $g_{j}(t_{j-1})=g_{j-1}(t_{j-1})$ and
$|Rm_{g_{j}(t)}|\leq\frac{C}{t}$ for all $j\geq 1$, where the constant $C$
depends only on $\alpha$ and the curvature bound of the original breather.
Then the immortal solution $g(t)$ is of Type III and is smooth by the
uniqueness of the Ricci flow (c.f. [3] and [7]).

Next, we fix a point $p_{0}\in M$ and define $x_{j}=\phi^{-{(j+1)}}(p_{0})$
for $j\geq 0$. Let $\ell_{+}$ be the reduced distance from $(p_{0},0)$, and we
shall proceed to show

$\displaystyle\limsup_{j\rightarrow\infty}\ell_{+}(x_{j},t_{j})<\infty.$ (4.4)

Let $\sigma:[0,1]\to M^{n}$ be a smooth curve satisfying $\sigma(0)=p_{0}$ and
$\sigma(1)=x_{0}=\phi^{-1}(p_{0})$. We define
$\sigma_{j}:[t_{j},t_{j+1}]\rightarrow M$ and $\gamma_{j}:[0,t_{j+1}]\to M$ as

$\displaystyle\sigma_{j}(t)=\phi^{-{(j+1)}}\circ\sigma(\alpha^{j+1}(t-t_{j})),\quad
t\in[t_{j},t_{j+1}],$ (4.5)
$\displaystyle\gamma_{j}(t)=\left\\{\begin{array}[]{ll}\sigma(t),&t\in[0,1],\\\
\sigma_{i}(t),&t\in[t_{i},t_{i+1}],0\leq i\leq j.\end{array}\right.$ (4.8)

Since

$\sigma_{j}(t_{j})=\phi^{-(j+1)}\circ\sigma(0)=\phi^{-j}\circ\sigma(1)=\phi^{-j}\circ\sigma(\alpha^{j}(t_{j}-t_{j-1}))=\sigma_{j-1}(t_{j}),$

we have that $\gamma_{j}(t)$ defined as (4.8) is a piecewise smooth $C^{0}$
curve satisfying $\gamma_{j}(0)=p_{0}$, $\gamma_{j}(t_{j})=x_{j}$. We then
compute

$\displaystyle 2\sqrt{t_{j+1}}\ell_{+}(x_{j+1},t_{j+1})$
$\displaystyle\leq\mathcal{L}_{+}(\sigma)+\sum\limits_{i=1}^{j}\int\limits_{t_{i}}^{t_{i+1}}\sqrt{t}\left(R(\sigma_{i}(t),t)+|\sigma_{i}^{\prime}(t)|^{2}_{g(t)}\right)dt$
$\displaystyle\leq
D+\sum\limits_{i=1}^{j}\int\limits_{t_{i}}^{t_{i+1}}\sqrt{t}\left(\frac{C}{t}+A\alpha^{i+1}\right)dt\leq
D+C\sum\limits_{i=1}^{j}\alpha^{-\frac{i+1}{2}}$ $\displaystyle\leq
D+C\alpha^{-\frac{j+1}{2}},$

for all $j\geq 0$, where
$A:=\max\limits_{t\in[0,1]}|\sigma^{\prime}(t)|_{g_{0}(t)}$, $C$, and $D$ are
all constants independent of $j$. Since $t_{j+1}\geq\alpha^{-(j+1)}$, we
obtain (4.4) from the above computation. Theorem 3.1 is now applicable to
$(M,g(t))_{t\in[0,\infty)}$ along the space-time sequence
$\\{(x_{j},t_{j})\\}_{j=1}^{\infty}$.

From the construction of $g(t)$, we easily observe that

$\Big{(}M,t_{j}^{-1}g(t_{j}t),x_{j}\Big{)}_{t\in\left[1,\frac{t_{j+1}}{t_{j}}\right]}\rightarrow\Big{(}M,g_{\infty}(t),p_{0}\Big{)}_{t\in[1,\alpha^{-1}]},$

where $g_{\infty}$ and $g_{0}$ differ only by a scaling and a time-shifting.
Furthermore, by Theorem 3.1, we have that there exists a function
$\ell_{+}^{\infty}:M\times[1,\alpha^{-1}]\rightarrow\mathbb{R}$, such that

$\displaystyle\ell_{+}^{j}\rightarrow\ell_{+}^{\infty}$

in the $C^{0,\alpha}_{\operatorname{loc}}$ sense and the weak
$*W^{1,2}_{\operatorname{loc}}$ sense, where
$\ell_{+}^{j}(\cdot,t)=\ell_{+}(\cdot,tt_{j})$. Our next goal is to show that
$\ell_{+}^{\infty}$ gives rise to an expander structure on
$(M,g_{\infty}(t))$.

Arguing as in section 4 of [9] or section 6 of [4] and applying the
monotonicity in theorem 2.2 in the same way as one have applied Perelman’s
reduced geometry, we have that
$\ell_{+}^{\infty}:M\times[1,\alpha^{-1}]\rightarrow\mathbb{R}$ is a smooth
function, satisfying

$\displaystyle\frac{\partial\ell_{+}^{\infty}}{\partial
t}-\Delta\ell_{+}^{\infty}+|\nabla\ell_{+}^{\infty}|^{2}+R_{\infty}+\frac{n}{2t}=0.$
(4.9)

Furthermore, by Theorem 3.1(1), we may find a function
$K^{\infty}:M\times[1,\alpha^{-1}]\rightarrow\mathbb{R}$, such that
$K^{j}\rightarrow K^{\infty}$ in the $C^{0,\alpha}_{\operatorname{loc}}$ sense
and the weak $*W^{1,2}_{\operatorname{loc}}$ sense. Fixing arbitrary
$1<s_{1}<s_{2}<\alpha^{-1}$ and an arbitrary nonnegative and smooth time-
independent cut-off function $\varphi$ compactly supported on $M$, we obtain
the following by (2.8)

$\displaystyle 0$
$\displaystyle\leq\int_{s_{1}}^{s_{2}}\int_{M}\frac{1}{t^{\frac{3}{2}}}K^{j}\varphi(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}^{j}}d\mu_{t}^{j}dt$
$\displaystyle\leq\int_{s_{1}}^{s_{2}}\int_{M}\left(\varphi\frac{\partial}{\partial
t}\ell_{+}^{j}-\langle\nabla\varphi,\nabla\ell_{+}^{j}\rangle+R_{j}\varphi+\frac{n}{2t}\varphi\right)(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}^{j}}d\mu_{t}^{j}dt$
$\displaystyle\rightarrow\int_{s_{1}}^{s_{2}}\int_{M}\left(\varphi\frac{\partial}{\partial
t}\ell_{+}^{\infty}-\langle\nabla\varphi,\nabla\ell_{+}^{\infty}\rangle+R_{\infty}\varphi+\frac{n}{2t}\varphi\right)(4\pi
t)^{-\frac{n}{2}}e^{-\ell_{+}^{\infty}}d\mu_{t}^{\infty}dt$ $\displaystyle=0.$

This shows that

$\displaystyle K^{\infty}\equiv 0$ everywhere on $\displaystyle
M\times[1,\alpha^{-1}].$ (4.10)

Since equations (2.2) and (2.3) are both carried to the limit in the sense of
distribution, and since $\ell_{+}^{\infty}$ is smooth, the following hold on
$M\times[1,\alpha^{-1}]$

$\displaystyle\frac{\partial}{\partial
t}\ell_{+}^{\infty}=R_{\infty}-\frac{\ell_{+}^{\infty}}{t},\ \ \
|\nabla\ell_{+}^{\infty}|^{2}=\frac{\ell_{+}^{\infty}}{t}-R_{\infty}.$

In combination with (4.9), we then have

$\displaystyle\frac{\partial\ell_{+}^{\infty}}{\partial
t}+\Delta\ell_{+}^{\infty}+|\nabla\ell_{+}^{\infty}|^{2}-R_{{\infty}}-\frac{n}{2t}=0,$
$\displaystyle
2\Delta\ell_{+}^{\infty}+|\nabla\ell_{+}^{\infty}|^{2}-R_{{\infty}}-\frac{\ell_{+}^{\infty}+n}{t}=0.$

Then, by Theorem 1.2 in [6],

$\displaystyle\left(\frac{\partial}{\partial t}+\Delta-
R_{\infty}\right)\left(t(2\Delta\ell_{+}^{\infty}+|\nabla\ell_{+}^{\infty}|^{2}-R_{\infty})-\ell_{+}^{\infty}-n\right)(4\pi
t)^{-\frac{n}{2}}e^{\ell_{+}^{\infty}}$ $\displaystyle=$
$\displaystyle-2t\left|Rc_{\infty}-\nabla\nabla\ell_{+}^{\infty}+\frac{g_{\infty}}{2t}\right|^{2}(4\pi
t)^{-\frac{n}{2}}e^{\ell_{+}^{\infty}}=0.$

Hence $\ell_{+}^{\infty}$ satisfies the following gradient expanding soliton
equation

$\displaystyle
Rc_{{\infty}}-\nabla\nabla\ell_{+}^{\infty}=-\frac{1}{2t}g_{\infty}.$

Note that Theorem 3.1(2) guarantees that $\displaystyle(4\pi
t)^{-\frac{n}{2}}e^{\ell_{+}^{\infty}}>0$ everywhere. This finishes the proof.

## References

  * [1] Simon Brendle. _A generalization of Hamilton’s differential Harnack inequality for the Ricci flow._ J. Differential Geom., 82(1):207-227, 2009
  * [2] H.D.Cao. _On Harnack’s inequalities for the Kaehler-Ricci flow_ , Invent. Math. 109 (1992), no. 2, 247-263
  * [3] B.L.Chen, X.P. Zhu, Uniqueness of the Ricci flow on complete noncompact manifolds. J. Differ.Geom. 74(1), 119–154 (2006)
  * [4] Liang Cheng, Yongjia Zhang._Perelman-type no breather theorem for noncompact Ricci flows._ Trans. Amer. math. Soc., to appear
  * [5] B.Chow, S.C.Chu, D.Glickenstein, C. Guenther, J. Isenberg, T. Ivey, D.Knopf, P.Lu, F.Luo, and L.Ni, _The Ricci flow: techniques and applications._ Parts I—IV. Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI.
  * [6] M.Feldman, T.Ilmanen, L.Ni, _Entropy and reduced distance for Ricci expanders_ , J. Geom. Anal. 15 (2005), 49-62.
  * [7] Kotschwar, B. _Backwards uniqueness for the ricci flow._ Int. Math. Res. Not. 2010(21), 4064–4097
  * [8] J.Lott. _On the long-time behavior of type-III Ricci flow solutions._ Mathematische Annalen, 2007, 339(3):627-666.
  * [9] Naber. A. _Noncompact Shrinking 4-Solitons with Nonnegative Curvature._ J.reine Angew.math, 2007, 645(2):125-153.
  * [10] Lei Ni,_Monotonicity and Li-Yau-Hamilton inequality_ , Surveys in Differential Geometry, Vol 12 (2008), pages 251–302.
  * [11] R.Hamilton, _The Harnack estimate for the Ricci flow_ , Journal of Differential Geometry, 1993, 37 (1): 225–243.
  * [12] Peng Lu and Yu Zheng. _New proofs of Perelman’s theorem on shrinking breathers in Ricci flow._ The Journal of Geometric Analysis, pages 1-7, 2017.
  * [13] G.Perelman, _The entropy formula for the Ricci flow and its geometric applications._ http://arxiv.org/abs/math/0211159.
  * [14] M. Rimoldi, G. Veronelli, _Extremals of Log Sobolev inequality on non-compact manifolds and Ricci soliton structures._ Calc. Var. 58, 66 (2019). https://doi.org/10.1007/s00526-019-1518-z
  * [15] Qi S. Zhang, _A no breathers theorem for some noncompact Ricci flows_. Asian Journal of Mathematics, 2012, 18(4).
  * [16] Yongjia, Zhang._A Note on Perelman’s No Shrinking Breather Theorem._ The Journal of Geometric Analysis (2018).


# Langevin dynamics of heavy quarks in a soft-hard factorized approach

Shuang Li<EMAIL_ADDRESS>College of Science, China Three Gorges University,
Yichang 443002, China  Key Laboratory of Quark and Lepton Physics (MOE),
Central China Normal University, Wuhan 430079, China  Fei Sun
<EMAIL_ADDRESS>College of Science, China Three Gorges University, Yichang
443002, China  Wei Xie<EMAIL_ADDRESS>College of Science, China Three
Gorges University, Yichang 443002, China  Wei Xiong<EMAIL_ADDRESS>College of Science, China Three Gorges University, Yichang 443002, China

###### Abstract

By utilizing a soft-hard factorized model, which combines a thermal
perturbative description of soft scatterings and a perturbative QCD-based
calculation for hard collisions, we study the energy and temperature
dependence of the heavy quark diffusion coefficients in Langevin dynamics. The
adjustable parameters are fixed from the comprehensive model-data comparison.
We find that a small value of the spatial diffusion coefficient at transition
temperature is preferred by data $2\pi TD_{s}(T_{c})\simeq 6$. With the
parameter-optimized model, we are able to describe simultaneously the prompt
$D^{0}$ $R_{\rm AA}$ and $v_{\rm 2}$ data at $p_{\rm T}\leq 8$ GeV in Pb–Pb
collisions at $\sqrt{s_{\rm NN}}=2.76$ and $\sqrt{s_{\rm NN}}=5.02$ TeV. We
also make predictions for non-prompt $D^{0}$ meson for future experimental
tests down to the low momentum region.

##  I INTRODUCTION

Ultrarelativistic heavy-ion collisions provide a unique opportunity to create
and investigate the properties of strongly interacting matter in extreme
conditions of temperature and energy density, where the normal matter turns
into a new form of nuclear matter, consisting of deconfined quarks and gluons,
namely quark-gluon plasma (QGP Bzdak  _et al._ (2020)). Such collisions allow
us to study the properties of the produced hot and dense partonic medium,
which are important for our understanding of the properties of the universe in
the first few milliseconds and the composition of the inner core of neutron
stars E. Shuryak (1980); Braun-Munzinger and Wambach (2009); Braun-Munzinger
_et al._ (2016). Over the past two decades, the measurements with heavy-ion
collisions have been carried at the Relativistic Heavy Ion Collider (RHIC) at
BNL and the Large Hadron Collider (LHC) at CERN Gyulassy and L. McLerran
(2005); E. Shuryak (2005); B. Muller, J. Schukraft, and B. Wyslouch (2012), to
search and explore the fundamental properties of QGP, notably its transport
coefficients related to the medium interaction of hard probes.

Heavy quarks (HQs), including charm and bottom, provide a unique insight into
the microscopic properties of QGP Zhou  _et al._ (2014); Tang  _et al._
(2014); Andronic  _et al._ (2016); X. Dong and V. Greco (2019); Dong  _et al._
(2019); Zhao  _et al._ (2020). Due to large mass, they are mainly produced in
initial hard scatterings and then traverse the QGP and experience elastic and
inelastic scatterings with its thermalized constituents Rapp and van Hees
(2010); F. Prino and R. Rapp (2016). The transport properties of HQ inside QGP
are encoded in the HQ transport coefficients, which are expected to affect the
distributions of the corresponding open heavy-flavor hadrons. The resulting
experimental observables like the nuclear modification factor $R_{\rm AA}$ and
elliptic anisotropy $v_{\rm 2}$ of various $D$ and $B$-mesons are therefore
sensitive to the HQ transport coefficients, in particular their energy and
temperature dependence. There now exist an extensive set of such measurements,
which allow a data-based extraction of these coefficients. In this work, we
make such an attempt by using a soft-hard factorized model (see Sec. III) to
calculate the diffusion and drag coefficients relevant for heavy quark
Langevin dynamics.

A particularly important feature of the QGP transport coefficients is their
momentum and temperature dependence, especially how they change within the
temperature region accessed by the RHIC and LHC experiments. For instance the
normalized jet transport coefficient $\hat{q}/T^{3}$ was predicted to present
a rapidly increasing behavior with decreasing temperature and develop a
near-$T_{c}$ peak structure J. F. Liao and E. Shuryak (2009). The subsequent
studies J. C. Xu, J. F. Liao and M. Gyulassy (2015, 2016); JET Collaboration
(2014); S. Z. Shi, J. F. Liao, and M. Gyulassy (2018, 2019); Ramamurti and
Shuryak (2018) seem to confirm this scenario. Another important transport
property, shear viscosity over entropy density ratio $\eta/s$, also presents a
visible $T$-dependence with a considerable increase above $T_{c}$ J. E.
Bernhard, J. S. Moreland and S. A. Bass (2019). Concerning the HQ diffusion
and drag coefficients, there are also indications of nontrivial temperature
dependence from both the phenomenological extractions S. S. Cao, G. Y. Qin,
and S. A. Bass (2015); Cao  _et al._ (2016); M. Greif, F. Senzel, H. Kremer,
K. Zhou, C. Greiner, and Z. Xu (2017); S. Li, C. W. Wang, X. B. Yuan, and S.
Q. Feng (2018); S. Li and C. W. Wang (2018); Prado  _et al._ (2020); Wang  _et
al._ (2020) and theoretical calculations M. He, R. J. Fries, and R. Rapp
(2013); T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya
(2016); W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M.
Nardi, and F. Prino (2011); O. Fochler, Z. Xu, and C. Greiner (2010); S. Z.
Shi, J. F. Liao, and M. Gyulassy (2018, 2019). The difference among the
derived hybrid models is mainly induced by the treatment of the scale of QCD
strong coupling, hadronization and non-perturbative effects Gossiaux (2019);
Cao (2021). See Ref. Beraudo  _et al._ (2018); Xu  _et al._ (2019); Cao  _et
al._ (2019) for the recent comparisons.

The paper is organized as follows. In Sec. II we introduce the general setup
of the employed Langevin dynamics. Section III is dedicated to the detailed
calculation of the heavy quark transport coefficients with the factorization
model. In Sec. IV we show systematic comparisons between modeling results and
data and optimize the model parameters based on global $\chi^{2}$ analysis.
With the parameter-optimized model, the energy and temperature dependence of
the heavy quark transport coefficients are presented in Sec. V, as well as the
comparisons with data and other theoretical calculations. Section VI contains
the summary and discussion.

##  II Langevin dynamics

The classicial Langevin Transport Equation (LTE) of a single HQ reads Rapp and
van Hees (2010)

$\displaystyle dx^{i}=\frac{p^{i}}{E}dt$ (1a) $\displaystyle
dp^{i}=-\eta_{D}p^{i}dt+C^{ik}\rho^{k}\sqrt{dt}.$ (1b)

with $i,k=1,2,3$. The first term on the right hand side of Eq. 1b represents
the deterministic drag force, $F^{i}_{drag}=-\eta_{D}p^{i}$, which is given by
the drag coefficient $\eta_{D}(E,T)$ with the HQ energy
$E=\sqrt{\vec{p}^{\;2}+m^{2}_{Q}}$ and the underlying medium temperature $T$.
The second term denotes the stochastic thermal force,
$F^{i}_{thermal}=C^{ik}\rho^{k}/\sqrt{dt}$, which is described by the momentum
argument of the covariance matrix $C^{ik}$,

$C^{ik}\equiv\sqrt{\kappa_{T}}(\delta^{ik}-\frac{p^{i}p^{k}}{\vec{p}^{\;2}})+\sqrt{\kappa_{L}}\frac{p^{i}p^{k}}{\vec{p}^{\;2}},$
(2)

together with a Gaussian-normal distributed random variable $\rho^{k}$,
resulting in the uncorrelated random momentum kicks between two different time
scales

$\displaystyle<\vec{F}^{i}_{thermal}(t)\cdot\vec{F}^{j}_{thermal}(t^{\prime})>_{\rho}$
(3) $\displaystyle=C^{ik}C^{jk}\delta(t-t^{\prime})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:LTEtensor})}}{{=}}\bigr{[}\kappa_{T}(\delta^{ij}-\frac{p^{i}p^{j}}{\vec{p}^{\;2}})+\kappa_{L}\frac{p^{i}p^{j}}{\vec{p}^{\;2}}\bigr{]}\delta(t-t^{\prime}).$

$\kappa_{T}$ and $\kappa_{L}$ are the transverse and longitudinal momentum
diffusion coefficients, respectively, which describe the momentum fluctuations
in the direction that perpendicular (i.e. transverse) and parallel (i.e.
longitudinal) to the propagation. Considering the Einstein relationship, which
enforces the drag coefficient starting from the momentum diffusion
coefficients as S. Li and J. F. Liao (2020)

$\displaystyle\eta_{D}=$
$\displaystyle\frac{\kappa_{L}}{2TE}+(\xi-1)\frac{\partial\kappa_{L}}{\partial\vec{p}^{\;2}}+\frac{d-1}{2\vec{p}^{\;2}}\bigr{[}\xi(\sqrt{\kappa_{T}}+\sqrt{\kappa_{L}})^{2}$
(4) $\displaystyle-(3\xi-1)\kappa_{T}-(\xi+1)\kappa_{L}\bigr{]}.$

The parameter $\xi$ denotes the discretization scheme of the stochastic
integral, which typically takes the values $\xi=0,0.5,1$, representing the
pre-point Ito, the mid-point Stratonovic, and the post-point discretization
schemes, respectively; $d=3$ indicates the spatial dimension. In the framework
of LTE the HQ-medium interactions are conveniently encoded into three
transport coefficients, i.e. $\eta_{D}$, $\kappa_{T}$ and $\kappa_{L}$ (Eq.
4). All the problems are therefore reduced to the evaluation of
$\kappa_{T/L}(E,T)$, which will be mainly discussed in this work.

Finally, we introduce a few detailed setups of the numerical implementation.
The space-time evolution of the temperature field and the velocity field are
needed to solve LTE (Eq. 1a and 1b). Following our previous analysis S. Li, C.
W. Wang, X. B. Yuan, and S. Q. Feng (2018); S. Li, C. W. Wang, R. Z. Wan, and
J. F. Liao (2019), they are obtained in a 3+1 dimensional viscous hydrodynamic
calculation I. Karpenko, P. Huoviven and M. Bleicher (2014), with the local
thermalization started at $\tau_{0}=0.6~{}{\rm fm}/{\it c}$, the shear
viscosity ${\eta/s=1/(4\pi)}$ and the critical temperature $T_{c}=165~{}{\rm
MeV}$ (see details in Ref. I. Karpenko, P. Huoviven and M. Bleicher (2014)).
When the medium temperature drops below $T_{c}$, heavy quark will hadronize
into the heavy-flavor hadrons via a fragmentation-coalescence approach. The
Braaten-like fragmentation functions are employed for both charm and bottom
quarks E. Braaten, K. Cheung, and T. C. Yuan (1993); M. Cacciari, P. Nason and
R. Vogt (2005). An instantaneous approach is utilized to characterize the
coalescence process for the formation of heavy-flavor mesons from the heavy
and (anti-)light quark pairs. The relevant coalescence probability is
quantified by the overlap integral of the Winger functions for the meson and
partons, which are defined through a harmonic oscillator and the Gaussian
wave-function C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi (1991),
respectively. See Ref. S. Li and J. F. Liao (2020) for more details.

##  III Momentum diffusion coefficients in a soft-hard factorized approach

When propagating throughout the QGP, the HQ scattering off the gluons and
(anti-)partons of the thermal deconfined medium, can be characterized as the
two-body elementary processes,

$Q~{}(p_{1})+i~{}(p_{2})\rightarrow Q~{}(p_{3})+i~{}(p_{4}),$ (5)

with $p_{1}=(E_{1},\vec{p}_{1})$ and $p_{2}$ are the four-momentum of the
injected HQ ($Q$) and the incident medium partons $i=q,g$, respectively, while
$p_{3}$ and $p_{4}$ are for the ones after scattering. Note that the medium
partons are massless ($m_{2}=m_{4}\sim 0$) in particular comparing with the
massive HQ ($m_{1}=m_{3}=m_{Q}$ in a few times $\rm GeV$). The corresponding
four-momentum transfer is
$(\omega,~{}\vec{q}^{\;})=(\omega,~{}\vec{q}_{T},~{}q_{L})$. The Mandelstam
invariants read

$\displaystyle s\equiv(p_{1}+p_{2})^{2}$ (6) $\displaystyle
t\equiv(p_{1}-p_{3})^{2}=\omega^{2}-q^{2}$ $\displaystyle
u\equiv(p_{1}-p_{4})^{2}$

with $q\equiv|\vec{q}\;|\ll E_{1}$ for small momentum exchange. The transverse
and longitudinal momentum diffusion coefficients can be determined by
weighting the differential interaction rate with the squared transverse and
longitudinal momentum transfer, respectively. It yields

$\displaystyle\kappa_{T}$ $\displaystyle=\frac{1}{2}\int
d{\Gamma}\;\vec{q}_{T}^{\;2}=\frac{1}{2}\int
d{\Gamma}\bigr{[}\omega^{2}-t-\bigr{(}\frac{2\omega
E_{1}-t}{2|\vec{p}_{1}|}\bigr{)}^{2}\bigr{]}$ (7)

and

$\displaystyle\kappa_{L}$ $\displaystyle=\int
d{\Gamma}\;q_{L}^{2}=\frac{1}{4\vec{p}_{1}^{\;2}}\int d{\Gamma}\;(2\omega
E_{1}-t)^{2}.$ (8)

As the momentum transfer vanishes ($|t|\rightarrow 0$), the gluon propagator
in the $t$-channel of the elastic process causes an infrared divergence in the
squared amplitude $\overline{|\mathcal{M}^{2}}|_{t-channel}\propto 1/t^{2}$,
which is usually regulated by a Debye screening mass, i.e. $t\rightarrow
t-\lambda m^{2}_{D}$ with an adjustable parameter $\lambda$ B. Svetitsky
(1988); P. B. Gossiaux and J. Aichelin (2008). Alternatively, it can be
overcome by utilizing a soft-hard factorized approach S. Peign$\acute{\rm e}$
and A. Peshier (2008a, b), which starts with the assumptions that the medium
is thermal and weakly coupled, and then the interactions between the heavy
quarks and the medium can be computed in thermal perturbation theory. Finally,
this approach allows to decompose the soft HQ-medium interactions with
$t>t^{\ast}$, from the hard ones with $t<t^{\ast}$. For soft collisions the
gluon propagator should be replaced by the hard-thermal loop (HTL) propagator
E. Braaten and T. C. Yuan (1991); J. P. Blaizot and E. Iancu (2002), while for
hard collisions the hard gluon exchange is considered and the Born
approximation is appropriate. Therefore, the final results of $\kappa_{T/L}$
include the contributions from both soft and hard components.

As discussed in the QED case S. Peign$\acute{\rm e}$ and A. Peshier (2008a),
$\mu+\gamma\rightarrow\mu+\gamma$, the complete calculation for the energy
loss of the energic incident heavy-fermion is independent of the intermediate
scale $t^{\ast}$ in high energy limit. However, in the QCD case, there is the
complication that the challenge of the validity of the HTL scenario,
$m_{D}^{2}\ll T^{2}$ S. Peign$\acute{\rm e}$ and A. Peshier (2008b), due to
the temperatures reached at RHIC and LHC energies. Consequently, in the QCD
case, the soft-hard approach is in fact not independent of the intermediate
scale $t^{\ast}$ P. B. Gossiaux and J. Aichelin (2008); W. M. Alberico, A.
Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011).
In this analysis we have checked that the calculations for $\kappa_{T/L}$ are
not sensitive to the choice of the artificial cutoff $t^{\ast}\sim m_{D}^{2}$.

In the next parts of this section, we will focus on the energy and temperature
dependence of the interaction rate ${\Gamma}$ at leading order in $g$ for the
elastic process, as well as the momentum diffusion coefficients (Eq. 7 and 8)
in soft and hard collisions, respectively.

###  III.1 $\kappa_{T/L}$ in soft region $t^{\ast}<t<0$

In soft collisions the exchanged four-momentum is soft, $\sqrt{-t}\sim gT$
($\lambda_{mfp}\sim 1/g^{2}T$ J. Ghiglieri, G. Moore, and D. Teaney (2016)),
and the $t$-channel long-wavelength gluons are screened by the mediums, thus,
they feel the presence of the medium and require the resummation. Here we just
show the final results, and the details are relegated to A. The transverse and
longitudinal momentum diffusion coefficients can be expressed as

$\displaystyle\kappa_{T}(E_{1},T)=$
$\displaystyle\frac{C_{F}g^{2}}{16\pi^{2}v_{1}^{3}}\int^{0}_{t^{\ast}}dt\;(-t)^{3/2}\int_{0}^{v_{1}}dx\frac{v_{1}^{2}-x^{2}}{(1-x^{2})^{5/2}}$
(9)
$\displaystyle\bigr{[}\rho_{L}(t,x)+(v_{1}^{2}-x^{2})\rho_{T}(t,x)\bigr{]}coth\bigr{(}\frac{x}{2T}\sqrt{\frac{-t}{1-x^{2}}}\bigr{)}$

and

$\displaystyle\kappa_{L}(E_{1},T)=$
$\displaystyle\frac{C_{F}g^{2}}{8\pi^{2}v_{1}^{3}}\int^{0}_{t^{\ast}}dt\;(-t)^{3/2}\int_{0}^{v_{1}}dx\frac{x^{2}}{(1-x^{2})^{5/2}}$
(10)
$\displaystyle\bigr{[}\rho_{L}(t,x)+(v_{1}^{2}-x^{2})\rho_{T}(t,x)\bigr{]}coth\bigr{(}\frac{x}{2T}\sqrt{\frac{-t}{1-x^{2}}}\bigr{)},$

respectively, with the HQ velocity $v_{1}=|\vec{p}|/E_{1}$ and the transverse
and longitudinal parts of the HTL gluon spectral functions444 The spectral
function involves only the low frequency excitations, namely the Landau cut,
while the quasiparticle excitations is irrelevant in this regime. See Eq.25
for more details. are given by

$\displaystyle\rho_{T}(t,x)=$ $\displaystyle\frac{\pi
m_{D}^{2}}{2}x(1-x^{2})\biggr{\\{}\bigr{[}-t+\frac{m_{D}^{2}}{2}x^{2}(1+\frac{1-x^{2}}{2x}ln\frac{1+x}{1-x})\bigr{]}^{2}$
(11) $\displaystyle+\bigr{[}\frac{\pi
m_{D}^{2}}{4}x(1-x^{2})\bigr{]}^{2}\biggr{\\}}^{-1}$

and

$\displaystyle\rho_{L}(t,x)=$ $\displaystyle\pi
m_{D}^{2}x\biggr{\\{}\bigr{[}\frac{-t}{1-x^{2}}+m_{D}^{2}(1-\frac{x}{2}ln\frac{1+x}{1-x})\bigr{]}^{2}$
(12) $\displaystyle+\bigr{(}\frac{\pi
m_{D}^{2}}{2}x\bigr{)}^{2}\biggr{\\}}^{-1}.$

###  III.2 $\kappa_{T/L}$ in hard region $t_{min}<t<t^{\ast}$

In hard collisions the exchanged four-momentum is hard, $\sqrt{-t}\gtrsim T$
($\lambda_{mfp}\sim 1/g^{4}T$ J. Ghiglieri, G. Moore, and D. Teaney (2016)),
and the pQCD Born approximation is valid in this regime. In analogy with the
previous part we give the $\kappa_{T/L}$ results directly, and the detailed
aspects of the calculations can be found in B. The momentum diffusion
coefficients reads

$\displaystyle\kappa^{Qi}_{T}(E_{1},T)=$
$\displaystyle\frac{1}{256\pi^{3}|\vec{p}_{1}|E_{1}}\int_{|\vec{p}_{2}|_{min}}^{\infty}d|\vec{p}_{2}|E_{2}n_{2}(E_{2})$
(13)
$\displaystyle\int_{-1}^{cos\psi|_{max}}d(cos\psi)\int_{t_{min}}^{t^{\ast}}dt\frac{1}{a}\bigr{[}-\frac{m_{1}^{2}(D+2b^{2})}{8\vec{p}_{1}^{\;2}a^{4}}$
$\displaystyle+\frac{E_{1}tb}{2\vec{p}_{1}^{\;2}a^{2}}-t(1+\frac{t}{4\vec{p}_{1}^{\;2}})\bigr{]}\;\overline{|\mathcal{M}^{2}}|^{Qi}$

and

$\displaystyle\kappa^{Qi}_{L}(E_{1},T)=$
$\displaystyle\frac{1}{256\pi^{3}|\vec{p}_{1}|^{3}E_{1}}\int_{|\vec{p}_{2}|_{min}}^{\infty}d|\vec{p}_{2}|E_{2}n_{2}(E_{2})$
(14)
$\displaystyle\int_{-1}^{cos\psi|_{max}}d(cos\psi)\int_{t_{min}}^{t^{\ast}}dt\frac{1}{a}\bigr{[}\frac{E_{1}^{2}(D+2b^{2})}{4a^{4}}$
$\displaystyle-\frac{E_{1}tb}{a^{2}}+\frac{t^{2}}{2}\bigr{]}\;\overline{|\mathcal{M}^{2}}|^{Qi}.$

The integrations limits and the short notations are shown in Eq. 45-51.

###  III.3 Complete results in soft-hard scenario

Combining the soft and hard contributions to the momentum diffusion
coefficients via

$\displaystyle\kappa_{T/L}(E,T)$
$\displaystyle=\kappa_{T/L}^{soft}(E,T)+\kappa_{T/L}^{hard}(E,T)$ (15)
$\displaystyle=\kappa_{T/L}^{soft}(E,T)+\sum_{i=q,g}\kappa_{T/L}^{hard-Qi}$

while $\kappa_{T/L}^{soft}(E\equiv E_{1},T)$ is given by Eq. 9 and 10, and
$\kappa_{T/L}^{hard-Qi}$ is expressed in Eq. 13 and 14 for a given incident
medium parton $i=q,g$. Adopting the post-point discretization scheme of the
stochastic integral, i.e. $\xi=1$ in Eq. 4, the drag coefficient
$\eta_{D}(E,T)$ can be obtained by inserting Eq. 15 into Eq. 4.

##  IV Data-based parameter optimization

Following the strategies utilized in our previous work S. Li and J. F. Liao
(2020); S. Li, W. Xiong, and R. Z. Wan (2020), the two key parameters in this
study, the intermediate cutoff $t^{\ast}$ and the scale $\mu$ of running
coupling (Eq. 31a), are tested within a wide range of possibility and drawn
constrains by comparing the relevant charm meson data with model results. We
calculate the corresponding final observable y for the desired species of
D-meson. Then, a $\chi^{2}$ analysis can be performed by comparing the model
predictions with experimental data

$\displaystyle\chi^{2}=\sum_{i=1}^{N}\biggr{(}\frac{\bf y^{\rm
Data}_{i}-y^{\rm Model}_{i}}{\sigma_{i}}\biggr{)}^{2}.$ (16)

In the above $\sigma_{i}$ is the total uncertainty in data points, including
the statistic and systematic components which are added in quadrature. $n=N-1$
denotes the degree of freedom ($d.o.f$) when there are $N$ data points used in
the comparison. In this study, we use an extensive set of LHC data in the
range $p_{\rm T}\leq 8~{}\rm{GeV}$: $D^{0}$, $D^{+}$, $D^{\ast+}$ and
$D_{s}^{+}$ $R_{\rm AA}$ data collected at mid-rapidity ($|y|<0.5$) in the
most central ($0-10\%$) and semi-central ($30-50\%$) Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=2.76~{}{\rm TeV}$ ALICE Collaboration (2016a, b) and
$\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$ ALICE Collaboration (2018a), as well as
the $v_{\rm 2}$ data in semi-central ($30-50\%$) collisions ALICE
Collaboration (2014); ALICE Collaboration (2018b); CMS Collaboration (2018a).

We scan a wide range of valus for ($t^{\ast}$, $\mu$):
$1\leq\frac{|t^{\ast}|}{m_{D}^{2}}\leq 3$ and $1\leq\frac{\mu}{\pi T}\leq 3$.
A total of 20 different combinations were computed and compared with the
experimental data. The obtained results are summarized in Tab. 1. The
$\chi^{2}$ values are computed separately for $R_{\rm AA}$ and $v_{\rm 2}$ as
well as for all data combined. To better visualize the results, we also show
them in Fig. 1, with left panels for $R_{\rm AA}$ analysis and right panels
for $v_{2}$ analysis. In both panels, the y-axis labels the desired
parameters, $\frac{|t^{\ast}|}{m_{D}^{2}}$ (upper) and $\frac{\mu}{\pi T}$
(lower), and x-axis labels the $\chi^{2}/d.o.f$ within the selected
ranges555$\chi^{2}/d.o.f$ is shown in the range $0.5<\chi^{2}/d.o.f<2.5$ for
better visualization.. The different points (filled gry circles) represent the
different combinations of parameters ($\frac{|t^{\ast}|}{m_{D}^{2}}$,
$\frac{\mu}{\pi T}$) in Tab. 1, with the number on top of each point to
display the relevant “$Model~{}ID$” for that model. A number of observations
can be drawn from the comprehensive model-data comparison. For the $R_{\rm
AA}$, several models achieve $\chi^{2}/d.o.f\simeq 1$ with
$\frac{|t^{\ast}|}{m_{D}^{2}}\lesssim 1.5$ and widespread values of $\mu$:
$1\leq\frac{\mu}{\pi T}\leq 3$. This suggests that $R_{\rm AA}$ appears to be
more sensitive to the intermediate cutoff while insensitive to the scale of
coupling constant. For the $v_{\rm 2}$, it clearly shows a stronger
sensitivity to $\mu$, which seems to give a better description
($\chi^{2}/d.o.f\simeq 1.5-2.0$) of the data with $\frac{\mu}{\pi T}=1$. It is
interesting to see that $R_{\rm AA}$ data is more powerful to constrain
$\frac{|t^{\ast}|}{m_{D}^{2}}$, while $v_{\rm 2}$ data is more efficient to
nail down $\frac{\mu}{\pi T}$. Taken all together, we can identify a
particular model that outperforms others in describing both $R_{\rm AA}$ and
$v_{\rm 2}$ data simultaneously with $\chi^{2}/d.o.f=1.3$. This one will be
the parameter-optimized model in this work: $|t^{\ast}|=1.5m_{D}^{2}$ and
$\mu=\pi T$.

$\it Model~{}ID$ | $|t^{\ast}|/m_{D}^{2}$ | $\mu/\pi T$ | $\chi^{2}/d.o.f$ | $\chi^{2}/d.o.f$ | Total
---|---|---|---|---|---
| (Cutoff) | (Scale) | ($R_{\rm AA}$) | ($v_{\rm 2}$) |
1 | 1.00 | 1.00 | 1.02 | 5.08 | 1.61
2 | 1.00 | 1.50 | 2.46 | 6.08 | 2.98
3 | 1.00 | 2.00 | 1.71 | 6.11 | 2.35
4 | 1.00 | 3.00 | 0.83 | 4.20 | 1.32
5 | 1.50 | 1.00 | 1.16 | 2.10 | 1.30
6 | 1.50 | 1.50 | 3.90 | 9.82 | 4.76
7 | 1.50 | 2.00 | 5.38 | 10.90 | 6.18
8 | 1.50 | 3.00 | 4.85 | 9.12 | 5.47
9 | 2.00 | 1.00 | 2.44 | 1.54 | 2.31
10 | 2.00 | 1.50 | 2.99 | 7.11 | 3.58
11 | 2.00 | 2.00 | 6.81 | 10.19 | 7.30
12 | 2.00 | 3.00 | 8.91 | 8.19 | 8.80
13 | 2.50 | 1.00 | 4.04 | 1.61 | 3.68
14 | 2.50 | 1.50 | 2.08 | 6.78 | 2.77
15 | 2.50 | 2.00 | 6.59 | 8.51 | 6.87
16 | 2.50 | 3.00 | 11.80 | 12.04 | 11.83
17 | 3.00 | 1.00 | 5.46 | 1.46 | 4.88
18 | 3.00 | 1.50 | 1.52 | 8.05 | 2.47
19 | 3.00 | 2.00 | 6.27 | 9.88 | 6.79
20 | 3.00 | 3.00 | 12.94 | 15.21 | 13.27
Table 1: Summary of the adjustable parameters in this work, together with the
relevant $\chi^{2}/d.o.f$ obtained for $R_{\rm AA}$ and $v_{\rm 2}$. Figure 1:
Comparision of $\chi^{2}/d.o.f$ based on the experimental data of $R_{\rm AA}$
(left) and $v_{\rm 2}$ (right). The model predictions are calculated by using
various combinations of parameter $|t^{\ast}|/m_{D}^{2}$ (upper) and $\mu/\pi
T$ (lower), which are represented as $\chi^{2}/d.o.f$ (x-axis) and the desired
parameter (y-axis), respectively. See the legend and text for details.

##  V Results

In this section we will first examine the $t^{\ast}$ dependence of
$\kappa_{T/L}$ (Eq. 15) for both charm and bottom quark. Then, the relevant
energy and temperature dependence of $\kappa_{T/L}(E,T)$ will be discussed
with the optimized parameters. For the desired observables we will perform the
comparisons with the results from lattice QCD at zero momentum limit, as well
as the ones from experimental data in the low to intermediate $p_{\rm T}$
region.

###  V.1 Energy and temperature dependence of the transport coefficients

In Fig. 2, charm (left) and bottom quark (right) $\kappa_{T}$ (thin curves)
and $\kappa_{L}$ (thick curves) are calculated, with
$m_{D}^{2}\leq{|t^{\ast}|}\leq 2.5m_{D}^{2}$ and $\mu=\pi T$, at the
temperature $T=0.40$ GeV (upper) and the energy $E=10.0$ GeV (lower).
$\kappa_{T}$ and $\kappa_{L}$ are, as expected, identical at zero momentum
limit ($E=m_{Q}$), while the latter one has a much stronger energy dependence
at larger momentum. Furthermore, $\kappa_{T/L}(E,T)$ behave a mild sensitivity
to the intermediate cutoff $t^{\ast}$ W. M. Alberico, A. Beraudo, A. De Pace,
A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011). Because the soft-hard
approach is strictly speaking valid when the coupling is small $m_{D}^{2}\ll
T^{2}$ S. Peign$\acute{\rm e}$ and A. Peshier (2008b), thus, the above
observations support the validity of this approach within the temperature
regions even though the coupling is not small.

Figure 2: (Color online) Comparison of the momentum diffusion coefficients
$\kappa_{T}$ (thin black curves) and $\kappa_{L}$ (thick red curves), for
charm (left) and bottom quarks (right), displaying separately the results
based various testing parameters: $m_{D}^{2}\leq{|t^{\ast}|}\leq 2.5m_{D}^{2}$
and $\mu=\pi T$.

In Fig. 3, charm quark $\kappa_{T}$ (left) and $\kappa_{L}$ (middle) are
evaluated with the optimized parameters, $|t^{\ast}|=1.5m_{D}^{2}$ and
$\mu=\pi T$, including both the soft (dotted blue curves) and hard
contributions (dashed black curves), at fixed temperature $T=0.40$ GeV (upper)
and at fixed energy $E=10.0$ GeV (lower). It is found that the soft components
are significant at low energy/temperature, while they are compatible at larger
values. The combined results (solid red curves) are presented as well for
comparison. With the post-point scheme ($\xi=1$ in Eq. 15), the drag
coefficients (right) behave (1) a nonmonotonic dependence, in particular on
the temperature, which is in part due to the improved treatment of the
screening in soft collisions, and in part due to the procedure of inferring
$\eta_{D}$ from $\kappa_{T/L}$ F. Prino and R. Rapp (2016); (2) a weak energy
dependence at $E\gtrsim 2m_{c}=3$ GeV. Similar conclusions can be drawn for
bottom quark, as shown in Fig. 4.

Figure 3: (Color online) Charm quark $\kappa_{T}$ (left) and $\kappa_{L}$
(middle) are shown at fixed temperature $T=0.40$ GeV (upper) and at fixed
energy $E=10.0$ GeV (lower), contributed by the soft (dotted blue curves) and
hard collisions (dashed black curves). The combined results (solid red curves)
are shown as well for comparison. The derived drag coefficient $\eta_{D}$
(right; Eq. 4) are obtained with the post-point scenario.

Figure 4: Same as Fig. 3 but for bottom quark.

Figure 5 presents the transport coefficient of charm quark,
$\hat{q}=2\kappa_{T}$, at fix momentum $p=10$ GeV (solid red curve). It can be
seen that $\hat{q}/T^{3}$ reaches the maximum near the critical temperature,
and then followed by a decreasing trend with $T$, providing a good description
of the light quark transport parameter (black circle points). The results from
various phenomenological extractions and theoretical calculations, including a
phenomenological fitting analysis with the Langevin-transport with Gluon
Radiation (LGR; dotted blue curve S. Li and J. F. Liao (2020)), a LO
calculation with a Linearized Boltzmann Diffusion Model (LIDO333LIDO results
are shown with only elastic scattering channels.; dashed black curve W. Y. Ke,
Y. R. Xu, and S. A. Bass (2018)), a nonperturbative treatment with Quasi-
Particle Model (QPM in Catania; dot-dashed green curve F. Scardina, S. K. Das,
V. Minissale, S. Plumari, and V. Greco (2017)), a novel confinement with semi-
quark-gluon-monopole plasma approach (CUJET3; shadowed red band J. C. Xu, J.
F. Liao and M. Gyulassy (2016); S. Z. Shi, J. F. Liao, and M. Gyulassy (2019,
2018)), are displayed as well for comparison. Similar temperature dependence
can be observed except the CUJET3 approach, which shows a strong enhancement
near $T_{c}$ regime.

Figure 5: (Color online) Transport coefficient, $\hat{q}/T^{3}(T)$, of charm
quark from the various calculations, including: the soft-hard factorized
approach (solid red curve), the LGR model with data optimized parameters
(dotted blue curve S. Li and J. F. Liao (2020)), a Bayesian ayalysis from LIDO
(dashed black curve W. Y. Ke, Y. R. Xu, and S. A. Bass (2018)), a quasi-
particle model from Catania (dot-dashed green curve G. Coci (2018)), CUJET3
(shadowed red band J. C. Xu, J. F. Liao and M. Gyulassy (2016); S. Z. Shi, J.
F. Liao, and M. Gyulassy (2019, 2018)) and JET Collaboration (black circle
points JET Collaboration (2014)) at $p=10$ GeV.

The scaled spatial diffusion coefficient describes the low energy interaction
strength of HQ in medium G. D. Moore and D. Teaney (2005),

$2\pi TD_{s}=\lim_{E\rightarrow m_{Q}}\frac{2\pi
T^{2}}{m_{Q}\cdot\eta_{D}(E,T)},$ (17)

and it can be calculated by substituting Eq. 4 into Eq. 17. The obtained
result for charm quark ($m_{c}=1.5$ GeV) is displayed as the solid red curve
in Fig. 6. It is found that a relatively strong increase of $2\pi TD_{s}(T)$
from crossover temperature $T_{c}$ toward high temperature. Meanwhile, the
$v_{\rm 2}$ data prefers a small value of $2\pi TD_{s}$ near $T_{c}$, $2\pi
TD_{s}(T_{c})\simeq 3-6$444This range is estimated from the testing models as
displayed in the right panels of Fig. 1., which is close to the lattice QCD
calculations D. Banerjee, S. Datta, R. Gavai, and P. Majumdar (2012); H. T.
Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner (2012); O.
Kaczmarek (2014); N. Brambilla, V. Leino, P. Petreczky, and A. Vairo (2020);
L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur, G. D.Moore, and H. T. Shu
(2020). The relevant results from other theoretical analyses, such as LGR S.
Li and J. F. Liao (2020), LIDO W. Y. Ke, Y. R. Xu, and S. A. Bass (2018),
Catania F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco (2017)
and CUJET3555CUJET3 results are obtained by performing the energy
interpolation down to $E=m_{Q}$. J. C. Xu, J. F. Liao and M. Gyulassy (2016);
S. Z. Shi, J. F. Liao, and M. Gyulassy (2019, 2018), show a similar trend but
with much weaker temperature dependence.

Figure 6: (Color online) Spatial diffusion constant $2\pi TD_{s}(T)$ of charm
quark from the sof-hard scenario (solid red curve), LGR (dotted blue curve S.
Li and J. F. Liao (2020), LIDO (dashed black curve W. Y. Ke, Y. R. Xu, and S.
A. Bass (2018)), Catania (dot-dashed green curve F. Scardina, S. K. Das, V.
Minissale, S. Plumari, and V. Greco (2017)), CUJET3 (shadowed red band J. C.
Xu, J. F. Liao and M. Gyulassy (2016); S. Z. Shi, J. F. Liao, and M. Gyulassy
(2019, 2018)) and lattice QCD calculations (black circle D. Banerjee, S.
Datta, R. Gavai, and P. Majumdar (2012), blue triangle H. T. Ding, A. Francis,
O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner (2012), pink square O.
Kaczmarek (2014), red inverted triangle N. Brambilla, V. Leino, P. Petreczky,
and A. Vairo (2020) and green plus L. Altenkort, A. M. Eller, O. Kaczmarek, L.
Mazur, G. D.Moore, and H. T. Shu (2020)). The result for D-meson (long dashed
pink curve L. Tolos and J. M. Torres-Rincon (2013)) in the hadronic phase is
shown for comparison.

###  V.2 Comparison with experimental data: $R_{\rm AA}$ and $v_{\rm 2}$

Figure 7 shows the $R_{\rm AA}$ of $D^{0}$ (a), $D^{+}$ (b), $D^{\ast+}$ (c)
and $D_{s}^{+}$ (d) in the most central ($0-10\%$) Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=2.76$ TeV, respectively. The calculations are done with
FONLL initial charm quark spectra and EPS09 NLO parametrization for the nPDF
in Pb S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018), and the pink band
reflects the theoretical uncertainties coming from these inputs. It can be
seen that, within the experimental uncertainties, the model calculations
provide a very good description of the measured $p_{\rm T}$-dependent $R_{\rm
AA}$ data for various charm mesons. Concerning the results in Pb–Pb collisions
at $\sqrt{s_{\rm NN}}=5.02$ TeV, as shown in Fig. 8, a good agreement is found
between the model and the measurement at $p_{\rm T}\lesssim 6~{}{\rm
GeV/\it{c}}$, while a slightly larger discrepancy observed at larger $p_{\rm
T}$.

Figure 7: (Color online) Comparison between experimental data (red box ALICE
Collaboration (2016a, b)) and soft-hard factorized model calculations (solid
black curve with pink uncertainty band) for the nuclear modification factor
$R_{\rm AA}$, of $D^{0}$ (a), $D^{+}$ (b), $D^{\ast+}$ (c) and $D_{s}^{+}$ (d)
at mid-rapidity ($|y|<0.5$) in central ($0-10\%$) Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=2.76$ TeV. Figure 8: Same as Fig. 7 but for Pb–Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. The data (solid ALICE
Collaboration (2018a), open Vertesi (2019)) are shown for comparison.

Figure 9 presents the elliptic flow coefficient $v_{\rm 2}$ of non-strange
D-meson (averaged $D^{0}$, $D^{+}$, and $D^{\ast+}$) in semi-central
($30-50\%$) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV (a) and
$\sqrt{s_{\rm NN}}=5.02$ TeV (b). Within the uncertainties of the experimental
data, our model calculations describe well the anisotropy of the transverse
momentum distribution of the non-strange D-meson. The sizable $v_{\rm 2}$ of
these charm mesons, in particular at intermediate $p_{\rm T}\simeq 3-5$ GeV,
suggests that charm quarks actively participate in the collective expansion of
the fireball.

Figure 9: (Color online) Comparison between experimental data (red ALICE
Collaboration (2014), black ALICE Collaboration (2018b) and blue boxes CMS
Collaboration (2018b)) and model calculations (solid black curve with pink
uncertainty band) for the elliptic flow $v_{\rm 2}$ of non-strange D-meson at
mid-rapidity ($|y|<0.5$) in semi-central ($30-50\%$) Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=2.76$ TeV (a) and $\sqrt{s_{\rm NN}}=5.02$ TeV (b).

In Fig. 10, the $p_{\rm T}$-differential $R_{\rm AA}$ of non-prompt $D^{0}$
mesons are predicted with the parameter-optimized model, and shown as a
function of $p_{\rm T}$ in central ($0-10\%$) Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV. Comparing with prompt $D^{0}$ $R_{\rm AA}$
(solid curve), at moderate $p_{\rm T}$ ($p_{\rm T}\sim 5-7$), a less
suppression behavior is found for $D^{0}$ from $B$-hadron decays (dashed
curve), reflecting a weaker in-medium energy loss effect of bottom quark
($m_{b}=4.75$ GeV), which has larger mass with respect to that of charm quark
($m_{c}=1.5$ GeV). We note that the future measurements performed for the non-
prompt $D^{0}$ $R_{\rm AA}$ and $v_{\rm 2}$, are powerful in nailing down the
varying ranges of the model parameters (see Sec. IV), which will largely
improve and extend the current understanding of the in-medium effects.

Figure 10: Comparison of non-prompt (dashed blue curve) and prompt $D^{0}$
(solid black curve) $R_{\rm AA}$ calculations as a function of $p_{\rm T}$
with the measured values (point ALICE Collaboration (2018a)), in the central
$0-10\%$ Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. See legend and text
for details.

##  VI Summary

In this work we have used a soft-hard factorized model to investigate the
heavy quark momentum diffusion coefficients $\kappa_{T/L}$ in the quark-gluon
plasma in a data-driven approach. In particular we’ve examined the validity of
this scenario by systematically scanning a wide range of possibilities. The
global $\chi^{2}$ analysis using an extensive set of LHC data on charm meson
$R_{\rm AA}$ and $v_{\rm 2}$ has allowed us to constrain the preferred range
of the two parameters: soft-hard intermediate cutoff $|t^{\ast}|=1.5m_{D}^{2}$
and the scale of QCD coupling constant $\mu=\pi T$. It is found that
$\kappa_{T/L}$ have a mild sensitivity to $t^{\ast}$, supporting the validity
of the soft-hard approach when the coupling is not small. With this
factorization model, we have calculated the transport coefficient $\hat{q}$,
drag coefficient $\eta_{D}$, spatial diffusion coefficient $2\pi TD_{s}$, and
then compared with other theoretical calculations and phenomenological
extractions. Our analysis suggests that a small value $2\pi TD_{s}\simeq 6$
appears to be much preferred near $T_{c}$. Finally we’ve demonstrated a
simultaneous description of charm meson $R_{\rm AA}$ and $v_{\rm 2}$
observables in the range $p_{\rm T}\leq 8$ GeV. We’ve further made predictions
for bottom meson observables in the same model.

We end with discussions on a few important caveats in the present study that
call for future studies:

  * $\bullet$

In the framework of Langevin approach, the heavy flavor dynamics is encoded
into three coefficient, $\kappa_{T}$, $\kappa_{L}$ and $\eta_{D}$, satisfying
Einstein’s relationship. It means that two of them are independent while the
third one can be obtained accordingly. The final results therefore depend on
the arbitrary choice of which of the two coefficients are calculated with the
employed model F. Prino and R. Rapp (2016); Xu  _et al._ (2019). For
consistency, we calculate independently the momentum diffusion coefficients
$\kappa_{T/L}$ with the factorization approach, and obtain the drag
coefficient via $\eta_{D}=\eta_{D}(\kappa_{T},\kappa_{L})$ (see Eq. 4). A
systematic study among the different options may help remedy this situation.

  * $\bullet$

According to the present modeling, the resulting spatial diffusion coefficient
exhibits a relatively strong increase of temperature comparing with lattice
QCD calculations, in particular at large $T$, as shown in Tab. 2. It
corresponds to a larger relaxation time and thus weaker HQ-medium coupling
strength. This comparison can be improved by (1) determining the key
parameters based on the upcoming measurements on high precision observables
(such as $R_{\rm AA}$, $v_{\rm 2}$ and $v_{\rm 3}$) of both $D$ and $B$-mesons
at low $p_{\rm T}$; (2) replacing the current $\chi^{2}$ analysis with a
state-of-the-art deep learning technology.

| $T=T_{c}$ | $T=2T_{c}$ | $T=3T_{c}$
---|---|---|---
This work |  $\sim$6 |  $\sim$22 |  $\sim$31
LQCD (median value) D. Banerjee, S. Datta, R. Gavai, and P. Majumdar (2012); N. Brambilla, V. Leino, P. Petreczky, and A. Vairo (2020) |  $\sim$5 |  $\sim$10 |  $\sim$13
Table 2: Summary of the different models for $2\pi TD_{s}$ at desired
temperature values.

  * $\bullet$

It is realized Cao  _et al._ (2020) that the heavy quark hadrochemistry, the
abundance of various heavy flavor hadrons, provides special sensitivity to the
heavy-light coalescence mechanism and thus plays an important role to
understand the observables like the baryon production and the baryon-to-meson
ratio. A systematic comparison including the charmed baryons over a broad
momentum region is therefore crucial for a better constraining of the model
parameters, as well as a better extraction of the heavy quark transport
coefficients in a model-to-data approach.

  * $\bullet$

The elastic scattering ($2\rightarrow 2$) processes between heavy quark and
QGP constituents are dominated for heavy quark with low to moderate transverse
momentum S. Li, C. W. Wang, R. Z. Wan, and J. F. Liao (2019). Thus, here we
consider only the elastic energy loss mechanisms to study the observables at
$p_{\rm T}\leq 8$ GeV. The missing radiative ($2\leftrightarrow 3$) effects
may help to reduce the discrepancy with $R_{\rm AA}$ data in the vicinity of
$p_{\rm T}=8$ GeV, as mentioned above (see Fig. 8). With the soft-hard
factorized approach, it would be interesting to include both elastic and
radiative contributions in a simultaneous best fit to data in the whole
$p_{\rm T}$ region. We also plan to explore this idea in the future.

###### Acknowledgements.

The authors are grateful to Andrea Beraudo and Jinfeng Liao for helpful
discussions and communications. We also thank Weiyao Ke and Gabriele Coci for
providing the inputs as shown in Fig. 5 and 6. This work is supported by the
Hubei Provincial Natural Science Foundation under Grant No.2020CFB163, the
National Science Foundation of China (NSFC) under Grant Nos.12005114, 11847014
and 11875178, and the Key Laboratory of Quark and Lepton Physics Contracts
Nos.QLPL2018P01 and QLPL201905.

##  Appendix A Derivation of the interaction rate and momentum diffusion
coefficients in soft collisions

In the QED case, $\mu+\gamma\rightarrow\mu+\gamma$, one can calculate relevant
interaction rate with small momentum transfer by using the imaginary part of
the muon self-energy $\Sigma(p_{1})$ H. A. Weldon (1983)

$\displaystyle{\Gamma}(E_{1},T)$
$\displaystyle=-\frac{1}{2E_{1}}\bar{n}_{F}(E_{1}){\rm
Tr}\bigr{[}({p\\!\\!\\!/}_{1}+m_{1}){\rm Im}\Sigma(p_{1})\bigr{]}.$ (18)

where, $p_{1}=(E_{1},\vec{p}_{1})$ and $m_{1}$ are the four-momentum and mass
of the injected muon, respectively 444The notations for the injected muon are
same with the ones for the injected HQ in the elastic process.. The trace term
in Eq. 18 was calculated with a resummed photon propagator, which is very
similar with the one in QCD 555The structure of the QED HTL-propagator follows
Eq. LABEL:eq:PropTL but with opposite signs S. Peign$\acute{\rm e}$ and A.
Peshier (2008a).. It is realized E. Braaten and M. H. Thoma (1991a, b) that,
in the QCD case, the contributions to ${\Gamma}$ can be obtained from the
corresponding QED calculations by simply substitution: $e^{2}\rightarrow
C_{F}g^{2}$, where, $e$ ($g$) is the QED (QCD) coupling constant and
$C_{F}=4/3$ is the quark Casimir factor. It yields E. Braaten and M. H. Thoma
(1991a)

$\displaystyle{\rm Tr}\bigr{[}({p\\!\\!\\!/}_{1}+m_{1}){\rm
Im}\Sigma(p_{1})\bigr{]}$ (19)
$\displaystyle=-2C_{F}g^{2}(1+e^{-E_{1}/T})\int_{q}\int
d\omega\;\bar{n}_{B}(\omega)\frac{E_{1}^{2}}{E_{3}}\mathcal{A}\mathcal{B},$

with

$\displaystyle\mathcal{A}\equiv(1-\frac{\omega+\vec{v}_{1}\cdot\vec{q}}{2E_{1}})\rho_{L}(\omega,q)+\bigr{[}\vec{v}_{1}^{\;2}(1-(\hat{v}_{1}\cdot\hat{q})^{2})$
(20)
$\displaystyle\qquad-\frac{\omega-\vec{v}_{1}\cdot\vec{q}^{\;}}{E_{1}}\bigr{]}\rho_{T}(\omega,q)$
$\displaystyle\mathcal{B}\equiv\bar{n}_{F}(E_{3})\delta(E_{1}-E_{3}-\omega)-n_{F}(E_{3})\delta(E_{1}+E_{3}-\omega)$
(21)

where, $\vec{v}_{1}=\vec{p}_{1}/E_{1}$ denotes the HQ velocity;
$n_{B/F}(E)=(e^{E/T}\mp 1)^{-1}$ indicates the thermal distributions for
Bosons/Fermions and $\bar{n}_{B/F}\equiv 1\pm n_{B/F}$ accounts for the Bose-
enhancement or Pauli-blocking effect. Here we have used the short notation
$\int_{q}\equiv\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}$ for phase space integrals.
Taking $|\vec{p}_{1}|$ and $m_{1/3}=m_{Q}$ are both much greater than the
underlying medium temperature, i.e. $|\vec{p}_{1}|,m_{1/3}\gg T$, thus,
$E_{1/3}\gg T$ and $n_{F}(E_{3})$ is exponentially suppressed and can be
dropped. Moreover, the first $\delta$ funciton in Eq. LABEL:eq:Trace_Soft3 can
be simplified since $E_{3}=\sqrt{(\vec{p}_{1}-\vec{q})^{2}+m_{1}^{2}}\approx
E_{1}-\vec{v}_{1}\cdot\vec{q}$. Concerning the second $\delta$ funciton, it
cannot contribute for $\omega$ less than or on the order of $T$ E. Braaten and
M. H. Thoma (1991a), which will be deleted in this work. Finally, Eq.
LABEL:eq:Trace_Soft2 and LABEL:eq:Trace_Soft3 can be reduced to

$\displaystyle\mathcal{A}\equiv\rho_{L}(\omega,q)+\vec{v}_{1}^{\;2}\bigr{[}1-(\hat{v}_{1}\cdot\hat{q})^{2}\bigr{]}\rho_{T}(\omega,q)$
(22a) $\displaystyle\mathcal{B}\equiv\delta(\omega-\vec{v}_{1}\cdot\vec{q}).$
(22b)

By substituting Eq. 22a and 22b into Eq. LABEL:eq:Trace_Soft1, one gets

$\displaystyle{\rm Tr}\bigr{[}({p\\!\\!\\!/}_{1}+m_{1}){\rm
Im}\Sigma(p_{1})\bigr{]}$ (23) $\displaystyle=-2C_{F}g^{2}E_{1}\int_{q}\int
d\omega\;\bar{n}_{B}(\omega)\delta(\omega-\vec{v}_{1}\cdot\vec{q}\;)$
$\displaystyle\qquad\biggr{\\{}\rho_{L}(\omega,q)+\vec{v}_{1}^{\;2}\bigr{[}1-(\hat{v}_{1}\cdot\hat{q})^{2}\bigr{]}\rho_{T}(\omega,q)\biggr{\\}}.$

Equation 18 can be rewritten as

$\displaystyle{\Gamma}(E_{1},T)$ (24) $\displaystyle=C_{F}g^{2}\int_{q}\int
d\omega\;\bar{n}_{B}(\omega)\delta(\omega-\vec{v}_{1}\cdot\vec{q}\;)$
$\displaystyle\qquad\biggr{\\{}\rho_{L}(\omega,q)+\vec{v}_{1}^{\;2}\bigr{[}1-(\hat{v}_{1}\cdot\hat{q})^{2}\bigr{]}\rho_{T}(\omega,q)\biggr{\\}}$

with the transverse and longitudinal spectral functions are given by the
imaginary part of the retarded propagator

$\displaystyle\rho_{T/L}(\omega,q)\equiv 2\cdot ImD_{T/L}^{R}(\omega,q).$ (25)

We note that, in the weak coupling limit, a consistent method is to use the
HTL resummed propagators, which is contributed by the quasiparticle poles and
the Landau damping cuts J. P. Blaizot and E. Iancu (2002). In this analysis,
we mainly focus on the low frequency excitation ($|\omega|<q$), where the
Landau damping is dominant and the quasiparticle excitation is irrelevant. The
resulting spectral function is therefore denoted by $\rho$ in Eq. 25.

The retarded propagator in Eq.25 reads

$\displaystyle D_{T/L}^{R}(\omega,q)\equiv{\Delta}_{T/L}(\omega+i\eta,q),$
(26)

which is defined by setting $q^{0}=\omega+i\eta$ $(\eta\rightarrow 0_{+})$,
i.e. the real energy, for the dressed gluon propagator
${\Delta}_{T/L}(q^{0},q)$ J. P. Blaizot and E. Iancu (2002)

$\displaystyle{\Delta}_{T}(q^{0},q)=\frac{-1}{(q^{0})^{2}-q^{2}-{\Pi}_{T}(x)}$
(27) $\displaystyle{\Delta}_{L}(q^{0},q)=\frac{-1}{q^{2}+{\Pi}_{L}(x)}.$

The medium effects are embedded in the HTL gluon self-energy

$\displaystyle{\Pi}_{T}(x)=\frac{m_{D}^{2}}{2}\bigr{[}x^{2}+(1-x^{2})Q(x)\bigr{]}$
(28) $\displaystyle{\Pi}_{L}(x)=m_{D}^{2}\bigr{[}1-Q(x)\bigr{]}$

where, $x=q^{0}/q$; $Q(x)$ is the Legendre polynomial of second kind

$\displaystyle Q(x)=\frac{x}{2}ln\frac{x+1}{x-1}$ (29)

and $m_{D}^{2}$ is the Debye screening mass squared for gluon

$\displaystyle m_{D}^{2}=g^{2}T^{2}\bigr{(}1+\frac{N_{f}}{6}\bigr{)}.$ (30)

The coupling constant, $g$, is quantified by the two-loop QCD beta-function O.
Kaczmarek and F. Zantow (2005)

$\displaystyle
g^{-2}(\mu)=2\beta_{0}ln(\frac{\mu}{\Lambda_{QCD}})+\frac{\beta_{1}}{\beta_{0}}ln\bigr{[}2ln(\frac{\mu}{\Lambda_{QCD}})\bigr{]}$
(31a) $\displaystyle\qquad\beta_{0}=\frac{1}{16\pi^{2}}(11-\frac{2}{3}N_{f})$
(31b)
$\displaystyle\qquad\beta_{1}=\frac{1}{(16\pi^{2})^{2}}(102-\frac{38}{3}N_{f})$
(31c)

where, $\pi T\leq\mu\leq 3\pi T$ and $\Lambda_{QCD}=261~{}{\rm MeV}$. $N_{f}$
is the number of active flavors in the QGP. Finally, for space-like momentum,
Eq. 25 can be expressed as

$\displaystyle\rho_{T}(\omega,q)=$ $\displaystyle\frac{\pi\omega
m_{D}^{2}}{2q^{3}}(q^{2}-\omega^{2})\biggr{\\{}\bigr{[}q^{2}-\omega^{2}$ (32)
$\displaystyle+\frac{\omega^{2}m_{D}^{2}}{2q^{2}}(1+\frac{q^{2}-\omega^{2}}{2\omega
q}ln\frac{q+\omega}{q-\omega})\bigr{]}^{2}$
$\displaystyle+\bigr{[}\frac{\pi\omega
m_{D}^{2}}{4q^{3}}(q^{2}-\omega^{2})\bigr{]}^{2}\biggr{\\}}^{-1}$

$\displaystyle\rho_{L}(\omega,q)=$ $\displaystyle\frac{\pi\omega
m_{D}^{2}}{q}\biggr{\\{}\bigr{[}q^{2}+m_{D}^{2}(1-\frac{\omega}{2q}ln\frac{q+\omega}{q-\omega})\bigr{]}^{2}$
(33) $\displaystyle+\bigr{(}\frac{\pi\omega
m_{D}^{2}}{2q}\bigr{)}^{2}\biggr{\\}}^{-1},$

with which Eq. LABEL:eq:Trace_Soft7 is computable.

The momentum diffusion coefficients can be calculated by subsituting Eqs.
LABEL:eq:Trace_Soft7, 32 and 33 back into Eq. 7 and 8, respectively, and then
performing the angular integral, yielding555 Note that the spectral functions
(Eq. 32 and 33) are odd, and the resulting
$\rho_{T/L}(-\omega,q)=-\rho_{T/L}(\omega,q)$ are used to obtain Eq. 34 and
35.

$\displaystyle\kappa_{T}(E_{1},T)=$
$\displaystyle\frac{C_{F}g^{2}}{8\pi^{2}v_{1}}\int_{0}^{q_{max}}dq\;q^{3}\int_{0}^{v_{1}q}d\omega\;(1-\frac{\omega^{2}}{v_{1}^{2}q^{2}})$
(34)
$\displaystyle\bigr{[}\rho_{L}(\omega,q)+(v_{1}^{2}-\frac{\omega^{2}}{q^{2}})\rho_{T}(\omega,q)\bigr{]}coth\frac{\omega}{2T}$

and

$\displaystyle\kappa_{L}(E_{1},T)=$
$\displaystyle\frac{C_{F}g^{2}}{4\pi^{2}v_{1}}\int_{0}^{q_{max}}dq\;q\int_{0}^{v_{1}q}d\omega\;\frac{\omega^{2}}{v_{1}^{2}}$
(35)
$\displaystyle\bigr{[}\rho_{L}(\omega,q)+(v_{1}^{2}-\frac{\omega^{2}}{q^{2}})\rho_{T}(\omega,q)\bigr{]}coth\frac{\omega}{2T}$

with $v_{1}\equiv|\vec{v}_{1}|$. The maximum momentum exchange is
$q_{max}=\sqrt{4E_{1}T}$ in the high-energy limit J. D. Bjorken (1982).

Next, we implement the further calculations by performing a simple change of
variables,

$\displaystyle t=\omega^{2}-q^{2}<0\quad\quad q^{2}=\frac{-t}{1-x^{2}}$ (36)
$\displaystyle
x=\frac{\omega}{q}<v_{1}~{}~{}~{}\qquad\quad\omega=x\sqrt{\frac{-t}{1-x^{2}}},$

resulting in

$\displaystyle
dtdx=\biggr{|}\frac{\partial(t,x)}{\partial(q,\omega)}\biggr{|}dqd\omega=2(1-x^{2})dqd\omega.$
(37)

Using Eq. LABEL:eq:OmegaQ2TX and 37, one arrives at Eq. 9-12. Similar results
can be found in Ref. W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M.
Monteno, M. Nardi, and F. Prino (2011); A. Beraduo, A. De Pace, W. M.
Alberico, and A. Molinari (2009).

##  Appendix B Derivation of interaction rate and momentum diffusion
coefficients in hard collisions

For two-boday scatterings, the transition rate is defined as the rate of
collisions with medium parton $i$, which changes the momentum of the HQ
(parton $i$) from $\vec{p}_{1}$ ($\vec{p}_{2}$) to
$\vec{p}_{3}=\vec{p}_{1}-\vec{q}$ ($\vec{p}_{4}=\vec{p}_{2}+\vec{q}$),

$\displaystyle\omega^{Qi}(\vec{p}_{1},\vec{q},T)=$
$\displaystyle\int_{p_{2}}n(E_{2})\bar{n}(E_{3})\bar{n}(E_{4})v_{rel}d\sigma^{Qi}(\vec{p}_{1},\vec{p}_{2}\rightarrow\vec{p}_{3},\vec{p}_{4}).$
(38)

Typically, one can assume $\bar{n}(E_{3})=1$ by neglecting the thermal effects
on the HQ after scattering. The differential cross section summed over the
spin/polarization and color of the final partons and averaged over those of
incident partons,

$\displaystyle
v_{rel}d\sigma^{Qi}(\vec{p}_{1},\vec{p}_{2}\rightarrow\vec{p}_{3},\vec{p}_{4})$
(39)
$\displaystyle=\frac{1}{2E_{1}}\frac{1}{2E_{2}}\frac{d^{3}\vec{p}_{3}}{(2\pi)^{3}2E_{3}}\frac{d^{3}\vec{p}_{4}}{(2\pi)^{3}2E_{4}}\overline{|\mathcal{M}^{2}}|^{Qi}(2\pi)^{4}$
$\displaystyle\qquad\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$

where, $v_{rel}=(\sqrt{(p_{1}\cdot p_{2})^{2}-(m_{1}m_{2})^{2}})/(E_{1}E_{2})$
is the relative velocity between the projectile HQ and the target parton. The
interaction rate for a given elastic process reads

$\displaystyle{\Gamma}^{Qi}(E_{1},T)=$ $\displaystyle\int
d^{3}\vec{q}~{}\omega^{Qi}(\vec{p}_{1},\vec{q},T)$ (40)
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:IndiTraRate},\ref{eq:DiffCroSec})}}{{=}}\frac{1}{2E_{1}}\int_{p_{2}}\frac{n(E_{2})}{2E_{2}}\int_{p_{3}}\frac{1}{2E_{3}}\int_{p_{4}}\frac{\bar{n}(E_{4})}{2E_{4}}$
$\displaystyle\quad\overline{|\mathcal{M}^{2}}|^{Qi}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$

The momentum diffusion coefficients can be obtained by inserting Eq. 40 into
Eq. 7 and 8, yielding

$\displaystyle\kappa^{Qi}_{T}(E_{1},T)$ (41)
$\displaystyle=\frac{1}{2E_{1}}\int_{p_{2}}\frac{n(E_{2})}{2E_{2}}\int_{p_{3}}\frac{1}{2E_{3}}\int_{p_{4}}\frac{\bar{n}(E_{4})}{2E_{4}}\frac{\vec{q}_{T}^{\;2}}{2}$
$\displaystyle\quad\theta(|t|-|t^{\ast}|)\overline{|\mathcal{M}^{2}}|^{Qi}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle=\frac{1}{256\pi^{4}|\vec{p}_{1}|E_{1}}\int_{|\vec{p}_{2}|_{min}}^{\infty}d|\vec{p}_{2}|E_{2}n(E_{2})\int_{-1}^{cos\psi|_{max}}d(cos\psi)$
$\displaystyle\quad\int_{t_{min}}^{t^{\ast}}dt\int_{\omega_{min}}^{\omega_{max}}d\omega\frac{\bar{n}({\omega+E_{2}})}{\sqrt{G(\omega)}}\bigr{[}\omega^{2}-t-\frac{(2E_{1}\omega-t)^{2}}{4\vec{p}_{1}^{2}}\bigr{]}\overline{|\mathcal{M}^{2}}|^{Qi}$

and

$\displaystyle\kappa^{Qi}_{L}(E_{1},T)$ (42)
$\displaystyle=\frac{1}{2E_{1}}\int_{p_{2}}\frac{n(\vec{p}_{2})}{2E_{2}}\int_{p_{3}}\frac{1}{2E_{3}}\int_{p_{4}}\frac{\bar{n}(\vec{p}_{4})}{2E_{4}}\bigr{(}\frac{2E_{1}\omega-t}{2|\vec{p}_{1}|}\bigr{)}^{2}$
$\displaystyle\quad\theta(|t|-|t^{\ast}|)\overline{|\mathcal{M}^{2}}|^{Qi}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle=\frac{1}{512\pi^{4}|\vec{p}_{1}|^{3}E_{1}}\int_{|\vec{p}_{2}|_{min}}^{\infty}d|\vec{p}_{2}|E_{2}n(E_{2})\int_{-1}^{cos\psi|_{max}}d(cos\psi)$
$\displaystyle\quad\int_{t_{min}}^{t^{\ast}}dt\int_{\omega_{min}}^{\omega_{max}}d\omega\frac{\bar{n}({\omega+E_{2}})}{\sqrt{G(\omega)}}(2E_{1}\omega-t)^{2}\overline{|\mathcal{M}^{2}}|^{Qi}$

Note that,

  * (1)

for hard collisions the momentum exchange is constrained by imposing
$\theta(|t|-|t^{\ast}|)$ in the first equality of Eq. LABEL:eq:KappaT3 and
LABEL:eq:KappaL3;

  * (2)

the tree level matrix elements squared includes the contributions from the
various channels, which are given in Ref. B. L. Combridge (1979):

    * $\bullet$

for $Q+q\rightarrow Q+q$ ($t$-channel only)

$\displaystyle\overline{|\mathcal{M}^{2}}|^{Qq}(s,t)$ (43)
$\displaystyle=2N_{f}\cdot 2N_{c}\cdot
g^{4}\frac{4}{9}\frac{(m_{1}^{2}-u)^{2}+(s-m_{1}^{2})^{2}+2m_{1}^{2}t}{t^{2}}$

    * $\bullet$

for $Q+g\rightarrow Q+g$ ($t$, $s$ and $u$-channel combined)

$\displaystyle\overline{|\mathcal{M}^{2}}|^{Qg}(s,t)$ (44)
$\displaystyle=2(N_{c}^{2}-1)g^{4}\biggr{[}2\frac{(s-m_{1}^{2})(m_{1}^{2}-u)}{t^{2}}$
$\displaystyle\quad+\frac{4}{9}\frac{(s-m_{1}^{2})(m_{1}^{2}-u)+2m_{1}^{2}(s+m_{1}^{2})}{(s-m_{1}^{2})^{2}}$
$\displaystyle\quad+\frac{4}{9}\frac{(s-m_{1}^{2})(m_{1}^{2}-u)+2m_{1}^{2}(m_{1}^{2}+u)}{(m_{1}^{2}-u)^{2}}$
$\displaystyle\quad+\frac{1}{9}\frac{m_{1}^{2}(4m_{1}^{2}-t)}{(s-m_{1}^{2})(m_{1}^{2}-u)}+\frac{(s-m_{1}^{2})(m_{1}^{2}-u)+m_{1}^{2}(s-u)}{t(s-m_{1}^{2})}$
$\displaystyle\quad-\frac{(s-m_{1}^{2})(m_{1}^{2}-u)-m_{1}^{2}(s-u)}{t(m_{1}^{2}-u)}\biggr{]}$

Concerning the degeneracy factors, in Eq. LABEL:eq:MatrixHQq, $2N_{f}$
reflects the identical contribution from all light quark and anti-quark
flavors, and $2N_{c}$ indicates the summing, rather than averaging, over the
helicities and colors of the incident light quark, while in Eq.
LABEL:eq:MatrixHQg, the factor $2(N_{c}^{2}-1)$ denotes the summing over the
polarization and colors of the incident gluon. The running coupling constant
takes $g(\mu)^{2}=4\pi\alpha_{s}(\mu)$, which is given by Eq. 31a with the
scale $\mu=\sqrt{-t}$.

  * (3)

with the help of $\delta$-function, we can reduce the integral in Eq.
LABEL:eq:KappaT3 and LABEL:eq:KappaL3 from 9-dimension (9D) to 4D in the
numerical calculations, by transforming the integration variables from
$(\vec{p}_{2},\vec{p}_{3},\vec{p}_{4})$ to $(|\vec{p}_{2}|,cos\psi,t)$, where
$\psi$ is the polar angle of $\vec{p}_{2}$; it yields the resutls as shown in
the second equality of Eq. LABEL:eq:KappaT3 and LABEL:eq:KappaL3; the relevant
limits of integration together with the additional notations are summarized
below:

$\displaystyle|\vec{p}_{2}|_{min}=\frac{|t^{\ast}|+\sqrt{(t^{\ast})^{2}+4m_{1}^{2}|t^{\ast}|}}{4(E_{1}+|\vec{p}_{1}|)}$
(45)

$\displaystyle
cos\psi|_{max}=min\biggr{\\{}1,\frac{E_{1}}{|\vec{p}_{1}|}-\frac{|t^{\ast}|+\sqrt{(t^{\ast})^{2}+4m_{1}^{2}|t^{\ast}|}}{4|\vec{p}_{1}|\cdot|\vec{p}_{2}|}\biggr{\\}}$
(46)

$\displaystyle t_{min}=-\frac{(s-m_{1}^{2})^{2}}{s}$ (47)

$\displaystyle\omega_{max/min}=\frac{b\pm\sqrt{D}}{2a^{2}}~{}with$ (48)

$\displaystyle\quad\quad a=\frac{s-m_{1}^{2}}{|\vec{p}_{1}|}$ (49)

$\displaystyle\quad\quad
b=-\frac{2t}{\vec{p}_{1}^{\;2}}\bigr{[}E_{1}(s-m_{1}^{2})-E_{2}(s+m_{1}^{2})\bigr{]}$
(50)

$\displaystyle\quad\quad
c=-\frac{t}{\vec{p}_{1}^{\;2}}\biggr{\\{}t\bigr{[}(E_{1}+E_{2})^{2}-s\bigr{]}+4\vec{p}_{1}^{\;2}\vec{p}_{2}^{\;2}sin^{2}\psi\biggr{\\}}$
(51)

$\displaystyle\quad\quad
D=b^{2}+4a^{2}c=-t\biggr{[}ts+(s-m_{1}^{2})^{2}\biggr{]}\cdot\biggr{(}\frac{4E_{2}sin\psi}{|\vec{p}_{1}|}\biggr{)}^{2}$
(52)

$\displaystyle G(\omega)=-a^{2}\omega+b\omega+c$ (53)

  * (4)

then, one can follow the procedure of Ref. S. Peign$\acute{\rm e}$ and A.
Peshier (2008a) for the analytical evaluation of $\omega$ integral. The
obtained results for $\kappa_{T}$ and $\kappa_{L}$ are shown in Eq. 13 and 14,
respectively.

## References

  * Bzdak  _et al._ (2020) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * E. Shuryak (1980) E. Shuryak, Phys. Rep. 61, 71 (1980).
  * Braun-Munzinger and Wambach (2009) P. Braun-Munzinger and J. Wambach, Rev. Mod. Phys. 81, 1031 (2009), arXiv:0801.4256 [hep-ph] .
  * Braun-Munzinger  _et al._ (2016) P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Phys. Rept. 621, 76 (2016), arXiv:1510.00442 [nucl-th] .
  * Gyulassy and L. McLerran (2005) M. Gyulassy and L. L. McLerran, Nucl. Phys. A 750, 30 (2005).
  * E. Shuryak (2005) E. Shuryak, Nucl. Phys. A 750, 64 (2005).
  * B. Muller, J. Schukraft, and B. Wyslouch (2012) B. Muller, J. Schukraft, and B. Wyslouch, Ann. Rev. Nucl. Part. Sci. 62, 361 (2012).
  * Zhou  _et al._ (2014) K. Zhou, N. Xu, Z. Xu, and P. Zhuang, Phys. Rev. C 89, 054911 (2014), arXiv:1401.5845 [nucl-th] .
  * Tang  _et al._ (2014) Z. Tang, N. Xu, K. Zhou, and P. Zhuang, J. Phys. G 41, 124006 (2014), arXiv:1409.5559 [nucl-th] .
  * Andronic  _et al._ (2016) A. Andronic _et al._ , Eur. Phys. J. C 76, 107 (2016), arXiv:1506.03981 [nucl-ex] .
  * X. Dong and V. Greco (2019) X. Dong and V. Greco, Prog. Part. Nucl. Phs. 104, 97 (2019).
  * Dong  _et al._ (2019) X. Dong, Y. Lee, and R. Rapp, Ann. Rev. Nucl. Part. Sci. 69, 417 (2019), arXiv:1903.07709 [nucl-ex] .
  * Zhao  _et al._ (2020) J. Zhao, K. Zhou, S. Chen, and P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020), arXiv:2005.08277 [nucl-th] .
  * Rapp and van Hees (2010) R. Rapp and H. van Hees, in  _Quark-Gluon Plasma 4 (World Sciencefic, Singapore, 2010)_ (2010) pp. 111–206.
  * F. Prino and R. Rapp (2016) F. Prino and R. Rapp, J. Phys. G 43, 093002 (2016).
  * J. F. Liao and E. Shuryak (2009) J. F. Liao and E. Shuryak, Phys. Rev. Lett. 102, 202302 (2009).
  * J. C. Xu, J. F. Liao and M. Gyulassy (2015) J. C. Xu, J. F. Liao and M. Gyulassy, Chin. Phys. Lett. 32, 092501 (2015).
  * J. C. Xu, J. F. Liao and M. Gyulassy (2016) J. C. Xu, J. F. Liao and M. Gyulassy, Journal of High Energy Physics 1602, 169 (2016).
  * JET Collaboration (2014) JET Collaboration, Phys. Rev. C 90, 014909 (2014).
  * S. Z. Shi, J. F. Liao, and M. Gyulassy (2018) S. Z. Shi, J. F. Liao, and M. Gyulassy, Chin. Phys. C 42, 104104 (2018).
  * S. Z. Shi, J. F. Liao, and M. Gyulassy (2019) S. Z. Shi, J. F. Liao, and M. Gyulassy, Chin. Phys. C 43, 044101 (2019).
  * Ramamurti and Shuryak (2018) A. Ramamurti and E. Shuryak, Phys. Rev. D 97, 016010 (2018), arXiv:1708.04254 [hep-ph] .
  * J. E. Bernhard, J. S. Moreland and S. A. Bass (2019) J. E. Bernhard, J. S. Moreland and S. A. Bass, Nat. Phys. 15, 1113 (2019).
  * S. S. Cao, G. Y. Qin, and S. A. Bass (2015) S. S. Cao, G. Y. Qin, and S. A. Bass, Phys. Rev. C 92, 024907 (2015).
  * Cao  _et al._ (2016) S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 94, 014909 (2016), arXiv:1605.06447 [nucl-th] .
  * M. Greif, F. Senzel, H. Kremer, K. Zhou, C. Greiner, and Z. Xu (2017) M. Greif, F. Senzel, H. Kremer, K. Zhou, C. Greiner, and Z. Xu, Phys. Rev. C 95, 054903 (2017).
  * S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng, Phys. Rev. C 98, 014909 (2018).
  * S. Li and C. W. Wang (2018) S. Li and C. W. Wang, Phys. Rev. C 98, 034914 (2018).
  * Prado  _et al._ (2020) C. A. G. Prado, W.-J. Xing, S. Cao, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 101, 064907 (2020), arXiv:1911.06527 [nucl-th] .
  * Wang  _et al._ (2020) S. Wang, W. Dai, B.-W. Zhang, and E. Wang,  (2020), 10.1088/1674-1137/abf4f5, arXiv:2012.13935 [nucl-th] .
  * M. He, R. J. Fries, and R. Rapp (2013) M. He, R. J. Fries, and R. Rapp, Phys. Rev. Lett. 110, 112301 (2013).
  * T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya (2016) T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 93, 034906 (2016).
  * W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011) W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino, Eur. Phys. J. C 71, 1666 (2011).
  * O. Fochler, Z. Xu, and C. Greiner (2010) O. Fochler, Z. Xu, and C. Greiner, Phys. Rev. C 82, 024907 (2010).
  * Gossiaux (2019) P. B. Gossiaux, Nucl. Phys. A 982, 113 (2019), arXiv:1901.01606 [nucl-th] .
  * Cao (2021) S. Cao, Nucl. Phys. A 1005, 121984 (2021).
  * Beraudo  _et al._ (2018) A. Beraudo _et al._ , Nucl. Phys. A 979, 21 (2018), arXiv:1803.03824 [nucl-th] .
  * Xu  _et al._ (2019) Y. Xu _et al._ , Phys. Rev. C 99, 014902 (2019), arXiv:1809.10734 [nucl-th] .
  * Cao  _et al._ (2019) S. Cao _et al._ , Phys. Rev. C C99, 054907 (2019), arXiv:1809.07894 [nucl-th] .
  * S. Li and J. F. Liao (2020) S. Li and J. F. Liao, Eur. Phys. J. C 80, 671 (2020).
  * S. Li, C. W. Wang, R. Z. Wan, and J. F. Liao (2019) S. Li, C. W. Wang, R. Z. Wan, and J. F. Liao, Phys. Rev. C 99, 054909 (2019).
  * I. Karpenko, P. Huoviven and M. Bleicher (2014) I. Karpenko, P. Huoviven and M. Bleicher, Comput. Phys. Commun. 185, 3016 (2014).
  * E. Braaten, K. Cheung, and T. C. Yuan (1993) E. Braaten, K. Cheung, and T. C. Yuan, Phys. Rev. D 48, R5049 (1993).
  * M. Cacciari, P. Nason and R. Vogt (2005) M. Cacciari, P. Nason and R. Vogt, Phys. Rev. Lett. 95, 122001 (2005).
  * C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi (1991) C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi, Phys. Rev. C 44, 1636 (1991).
  * B. Svetitsky (1988) B. Svetitsky, Phys. Rev. D 37, 2484 (1988).
  * P. B. Gossiaux and J. Aichelin (2008) P. B. Gossiaux and J. Aichelin, Phys. Rev. C 78, 014904 (2008).
  * S. Peign$\acute{\rm e}$ and A. Peshier (2008a) S. Peign$\acute{\rm e}$ and A. Peshier, Phys. Rev. D 77, 014015 (2008a).
  * S. Peign$\acute{\rm e}$ and A. Peshier (2008b) S. Peign$\acute{\rm e}$ and A. Peshier, Phys. Rev. D 77, 114017 (2008b).
  * E. Braaten and T. C. Yuan (1991) E. Braaten and T. C. Yuan, Phys. Rev. Lett. 2183, 66 (1991).
  * J. P. Blaizot and E. Iancu (2002) J. P. Blaizot and E. Iancu, Phys. Rep. 359, 355 (2002).
  * J. Ghiglieri, G. Moore, and D. Teaney (2016) J. Ghiglieri, G. Moore, and D. Teaney, J. High Energy Phys. 03, 095 (2016).
  * S. Li, W. Xiong, and R. Z. Wan (2020) S. Li, W. Xiong, and R. Z. Wan, Eur. Phys. J. C 80, 1113 (2020).
  * ALICE Collaboration (2016a) ALICE Collaboration, Journal of High Energy Physics 03, 081 (2016a).
  * ALICE Collaboration (2016b) ALICE Collaboration, Journal of High Energy Physics 03, 082 (2016b).
  * ALICE Collaboration (2018a) ALICE Collaboration, Journal of High Energy Physics 10, 174 (2018a).
  * ALICE Collaboration (2014) ALICE Collaboration, Phys. Rev. C 90, 034904 (2014).
  * ALICE Collaboration (2018b) ALICE Collaboration, Phys. Rev. Lett. 120, 102301 (2018b).
  * CMS Collaboration (2018a) CMS Collaboration, Phys. Rev. Lett. 120, 202301 (2018a).
  * W. Y. Ke, Y. R. Xu, and S. A. Bass (2018) W. Y. Ke, Y. R. Xu, and S. A. Bass, Phys. Rev. C 98, 064901 (2018).
  * F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco (2017) F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco, Phys. Rev. C 96, 044905 (2017).
  * G. Coci (2018) G. Coci, PhD Thesis (2018).
  * G. D. Moore and D. Teaney (2005) G. D. Moore and D. Teaney, Phys. Rev. C 71, 064904 (2005).
  * D. Banerjee, S. Datta, R. Gavai, and P. Majumdar (2012) D. Banerjee, S. Datta, R. Gavai, and P. Majumdar, Phys. Rev. D 85, 014510 (2012).
  * H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner (2012) H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner, Phys. Rev. D 86, 014509 (2012).
  * O. Kaczmarek (2014) O. Kaczmarek, Nucl. Phys. A 931, 633 (2014).
  * N. Brambilla, V. Leino, P. Petreczky, and A. Vairo (2020) N. Brambilla, V. Leino, P. Petreczky, and A. Vairo, Phys. Rev. D 102, 074503 (2020).
  * L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur, G. D.Moore, and H. T. Shu (2020) L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur, G. D.Moore, and H. T. Shu, Phys. Rev. D 103, 014511 (2020).
  * L. Tolos and J. M. Torres-Rincon (2013) L. Tolos and J. M. Torres-Rincon, Phys. Rev. D 88, 074019 (2013).
  * Vertesi (2019) R. Vertesi (ALICE), in  _18th Conference on Elastic and Diffractive Scattering_ (2019) arXiv:1910.01981 [nucl-ex] .
  * CMS Collaboration (2018b) CMS Collaboration, Phys. Rev. Lett. 120, 202301 (2018b).
  * Cao  _et al._ (2020) S. Cao, K.-J. Sun, S.-Q. Li, S. Y. F. Liu, W.-J. Xing, G.-Y. Qin, and C. M. Ko, Phys. Lett. B 807, 135561 (2020), arXiv:1911.00456 [nucl-th] .
  * H. A. Weldon (1983) H. A. Weldon, Phys. Rev. D 28, 2007 (1983).
  * E. Braaten and M. H. Thoma (1991a) E. Braaten and M. H. Thoma, Phys. Rev. D 44, 1298 (1991a).
  * E. Braaten and M. H. Thoma (1991b) E. Braaten and M. H. Thoma, Phys. Rev. D 44, R2625(R) (1991b).
  * O. Kaczmarek and F. Zantow (2005) O. Kaczmarek and F. Zantow, Phys. Rev. D 71, 114510 (2005).
  * J. D. Bjorken (1982) J. D. Bjorken, Fermilab Report No. PUB-82/59-THY (1982).
  * A. Beraduo, A. De Pace, W. M. Alberico, and A. Molinari (2009) A. Beraduo, A. De Pace, W. M. Alberico, and A. Molinari, Nucl. Phys. A 831, 59 (2009).
  * B. L. Combridge (1979) B. L. Combridge, Nucl. Phys. B 151, 429 (1979).