Effective field theory approach to open heavy flavor production in heavy-ion collisions

We develop a version of Soft Collinear Effective Theory (SCET) which includes finite quark masses, as well as Glauber gluons that describe the interaction of collinear partons with QCD matter. In the framework of this new effective field theory, labeled SCET$_{\mathrm{M,G}}$, we derive the massive splitting functions in the vacuum and the QCD medium for the processes $Q\to Qg$, $Q\to gQ$ and $g\to Q\bar Q$. The numerical effects due to finite quark masses are sizable and our results are consistent with the traditional approach to parton energy loss in the soft gluon emission limit. In addition, we present a new framework for including the medium-induced full splitting functions consistent with next-to-leading order calculations in QCD for inclusive hadron production. Finally, we show numerical results for the suppression of $D$- and $B$-mesons in heavy ion collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV and 2.76 TeV and compare to available data from the LHC.


Introduction
Inclusive open heavy flavor production in proton-proton and heavy-ion collisions is considered to be one of the most important tests of our understanding of QCD and, in particular, of quark mass effects both in the vacuum and in a QCD medium. Measurements of heavy flavor meson cross sections have been performed at the Tevatron [1][2][3] and at the Relativistic Heavy Ion Collider (RHIC) [4][5][6]. More recently, the ATLAS, CMS and ALICE experimental collaborations at the Large Hadron Collider (LHC) have provided high precision data [7][8][9][10][11][12][13][14] for several center of mass (CM) energies, and they will continue to extend the currently existing suite of data sets with future measurements. For a recent review on heavy flavor production, see Ref. [15]. Theoretical and experimental advances in understanding the nuclear modification of light hadrons, heavy mesons, as well as jets and jet substructure in nucleus-nucleus reactions have been a highlight of the heavy ion programs at RHIC and the LHC [16]. Such highly energetic particles and jets are powerful and valuable probes of the quark-gluon plasma (QGP) produced in these collisions. In particular, open heavy flavor production plays a crucial role in elucidating the properties of QGP and has received growing attention from the experimental and theoretical communities in recent years. Earlier data from RHIC and preliminary measurements from the CMS collaboration [13] at √ s NN = 5.02 TeV show that the suppression rates for D 0 mesons are in fact the same as for light charged hadrons within the experimental uncertainty, contrary to the early expectation from traditional parton energy loss framework [17]. This has stimulated a series of work addressing the interaction of heavy quarks with the QCD medium in the literature, see e.g. [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] and references therein. Different than light hadrons, the heavy quark mass introduces an additional perturbative scale besides the large transverse momentum p T at which the heavy meson is produced. This feature makes the description of heavy quark dynamics both challenging and also particularly interesting, as it can reveal unique information about the QCD medium.
In the past few year, there has been a new development in the theoretical description of observables in heavy-ion collisions. In [34,35], the powerful techniques of Soft Collinear Effective Theory (SCET) [36][37][38] were first applied to describe the interactions of hard probes (highly energetic particles and jets) with the QCD medium. The underlying idea is to include a Glauber mode that describes the interaction of highly energetic partons with the QCD medium. The effective field theory describing the medium interactions is commonly referred to as SCET G . In [35,[39][40][41] the in-medium splitting functions were derived to first order in opacity. By taking the soft emission limit, the full in-medium splitting functions reduce to the results obtained within traditional approaches to parton energy loss [42,43]. Using SCET G allows to systematically go beyond these traditional approaches. In [44][45][46][47] several applications of the in-medium splitting functions were developed. In [44,45] a medium-modified DGLAP evolution approach was successfully applied to describe the suppression of light charged hadrons in the medium. In [46,47] the SCET G based splitting functions were used to describe both inclusive jet production and jet substructure observables in heavy ion collisions.
In this paper, we perform the next logical step in this line of work by including finite mass effects in the SCET G Lagrangian and, thus, enable the effective theory study the of interactions of heavy quarks with the QCD medium. The SCET Lagrangian in the vacuum with quark masses was first derived in [48]. The corresponding theory in the vacuum is commonly referred to as SCET M . Consequently, we label the new effective field theory presented in this work SCET M,G . With this new theory at hand, we extend the in-medium splitting functions to the massive case. The newly derived results can be used to describe the suppression of open heavy flavor production in heavy-ion collisions. We introduce a new way to implement the in-medium corrections consistently at next-toleading order (NLO) in perturbative QCD. This can be achieved by formally introducing medium-modified fragmentation functions based on the SCET M,G splitting functions. We present first numerical results in this work and compare to data taken at the LHC. As it turns out, the description of the underlying proton-proton baseline plays an important role. Several different approaches are available in the literature to deal with heavy quark masses in the fragmentation process [49][50][51][52][53][54][55][56]. The suppression rates in heavy-ion collisions crucially depend on whether the heavy meson is produced by a fragmenting heavy quark or a gluon. Gluons lose more energy than heavy quarks when they undergo multiple scatterings and splittings in the medium before they eventually fragment into the observed heavy meson. In this work we analyze the different production mechanisms and study the associated suppression rates in Pb+Pb collisions.
The remainder of this paper is organized as follows. In section 2, we derive the Lagrangian of the new effective field theory SCET M,G that includes both finite quark masses and Glauber gluons that describe the interaction with the QCD medium. We derive the massive vacuum and in-medium splitting functions and make the connection to previous results in the literature in the soft emission limit. In addition, we study the numerical impact of finite quark masses. In section 3, we start by introducing our new formalism that treats both vacuum and in-medium corrections consistently to NLO in QCD. We present numerical results for both proton-proton and heavy-ion collisions and compare to currently available data from the LHC. In section 4, we conclude and give an outlook.

Effective field theory for massive quarks in the medium
In this section, we first introduce the basic SCET ingredients relevant to our calculation. We then derive the Lagrangian for the effective theory SCET M,G , which describes the interactions of collinear heavy quarks with the QCD medium through Glauber gluon exchange. We use this new version of SCET to derive the massive splitting kernels in the vacuum and in nuclear matter for the splitting processes Q → Qg, Q → gQ and g → QQ, where Q represents a heavy quark. In order to establish the connection between our newly derived in-medium splitting functions and the results derived previously within the traditional parton energy loss approach, we further consider the limit of the splitting kernels where the emitted parton becomes soft. Finally, we present numerical results for the splitting functions and compare to the massless case [39], as well as the results from traditional parton energy loss calculations.

Basic SCET ingredients
SCET [36][37][38] is an effective field theory describing the dynamics of soft and collinear quarks and gluons in the presence of hard interactions. SCET has been applied successfully to hard-scattering processes at the LHC, in particular to the production of highly energetic hadrons and jets. We adopt the following convention for the light-cone notation. We define two reference vectors n µ = (1, 0, 0, 1),n µ = (1, 0, 0, −1) satisfying n 2 =n 2 = 0 and n·n = 2. Any four-vector p µ can be written as p µ = (p + , p − , p ⊥ ) = (n · p, n · p, p ⊥ ). In other words, The hierarchy between the hard, collinear and the soft scale is determined by the SCET power counting parameter λ. The SCET degrees of freedom have the following momentum scaling using light-cone coordinates p µ = (p + , p − , p ⊥ ) ∼ p + (1, λ 2 , λ) for a collinear mode and p µ ∼ p + (λ 2 , λ 2 , λ 2 ) for a soft mode. Note that in our notation, p + is a hard scale. The Glauber modes that we consider in the next section scale as p µ ∼ p + (λ 2 , λ 2 , λ). For a collinear quark, one separates the momentum as p =p + k, wherep =n(n · p)/2 + p ⊥ is the large label momentum and k is the residual momentum. The label momentum component is removed by defining a quark field ψ n,p (x) through where ψ(x) is the standard QCD quark field. The four component field ψ n,p (x) has two large components ξ n,p (x) and two small components ξn ,p (x). One defines the following two projections with ψ n,p (x) = ξ n,p (x)+ξn ,p (x). Operators in SCET are defined in terms of gauge invariant quark and gluon fields which are given by with iD µ n⊥ = P µ n⊥ + gA µ n⊥ and P µ is the label momentum operator. Furthermore, W n is a Wilson line of collinear gluons, withP =n · P.

SCET with quark masses and Glauber gluons
When an energetic parton traverses a dense and/or hot QCD medium, as produced in heavy-ion collisions, the formation of an in-medium parton shower can be described using perturbative methods. Following [57], the medium can be thought of as color-screened quasi-particles that generate a Coulomb-like potential that effectively leads to a background field for the partons traveling through the QCD medium. The energetic parton that eventually produces a jet of particles undergoes multiple elastic interactions with the quasi-particles. In the vacuum, the parton shower forms by standard soft and collinear splittings. These processes are described by the original SCET Lagrangian. In the medium, one needs to consider additional medium-induced splitting processes. The interaction of collinear massless quarks and gluons with t-channel off-shell gluons is described by the SCET G Lagrangian as derived in [34,35]. So far, SCET G only contains the interaction of collinear quarks and gluons with the medium through the so-called Glauber gluon exchange, which has the momentum scaling ∼ p + (λ 2 , λ 2 , λ). The corresponding soft sector has not been derived yet, nor the interactions between the Glauber gluons and soft gluons and quarks [58]. In general, Glauber modes play an essential role in various aspects of QCD, see for example [58][59][60] for more details. We leave the complete formulation of SCET G with soft modes in the medium for future work. The purpose of this section is to derive an extension of the effective field theory SCET G , as presented in [35] for massless quarks and gluons, by including the effects of heavy quark masses. The resulting new version of the effective field theory, labeled SCET M,G , will enable us to study the interactions of energetic heavy quarks with a hot QCD medium as a first example. We start by reviewing the previous work [35], in which several source fields and gauges were considered and the scaling of the in-medium background field was derived. The structure of the SCET G Lagrangian is where L SCET (ξ n , A n ) is the vacuum SCET Lagrangian, with ξ n and A n the collinear quark and gluon fields, respectively. The second term L G (ξ n , A n , A G ) is given by which contains the interactions between the Glauber gluons A G and the collinear quarks and gluons traversing the QCD medium. This result corresponds to the static source as described in more detail in [35]. Moreover, the result presented here for L SCET G corresponds to the so-called hybrid gauge, where a different gauge is chosen for the collinear gluons (light-cone gauge) than for the Glauber gluons (covariant R ξ gauge). This is a valid gauge choice as it was shown in [35] since both sectors of the effective theory SCET G are separately gauge invariant. This gauge choice simplifies the calculations of the in-medium splitting functions considerably as presented in the next section. In the hybrid gauge both the collinear Wilson line as well as the transverse gauge link which appear in the effective field theory reduce to unity. The corresponding Feynman rules for the interaction of collinear massless quarks and gluons with the Glauber modes can be obtained directly from Eq. (2.7). We now present the extension of the above effective field theory SCET G by including finite quark masses. The vacuum SCET Lagrangian involving finite quark masses was first derived in [48]. As mentioned above, the corresponding effective field theory is typically denoted by SCET M . In this section, we recall its derivation and, in addition, we include Glauber modes that describe the interaction with the QCD medium. The resulting new effective field theory that involves both massive quarks and Glauber modes describing the interaction with the medium is denoted by SCET M,G . We only focus on the collinear quark sector of the Lagrangian density for which we now take into account finite quark masses. We start from the standard QCD Lagrangian where the covariant derivative is given by iD µ = ∂ µ +gA µ . Here the gauge field A µ consists of three contributions where A µ c,s,G are the collinear, the soft, and the Glauber gluon gauge fields, respectively. The collinear and soft fields scale like the corresponding collinear and soft momenta as described above. The Glauber gluon has the momentum scaling as p µ ∼ p + (λ 2 , λ 2 , λ). However, the scaling of the corresponding Glauber gluon field A G is different and needs to be derived by expressing the Glauber gluon field in terms of the QCD current of the source and the gluon propagator. By working out the scaling of every term in the resulting expression, the scaling for A G can be obtained which eventually also depends on the gauge choice, see [34,35] for more details. In the hybrid gauge we use the covariant R ξ gauge for the medium Glauber gluons. In this case the corresponding Glauber gluon field scales as A G ∼ p + (λ 2 , λ 2 , λ 3 ) as it was derived in [35]. We start by substituting Eqs. (2.2) and (2.3) into the QCD Lagrangian in Eq. (2.8), and we find We can now integrate out the small component of the collinear quark field ξn ,p by making use of the equation of motion With this relation, we obtain from Eq. (2.10) the following result for the leading-order (in λ) SCET M,G Lagrangian density where we introduced the label momentum operator also in the exponential for notational convenience and we have The mass-dependent terms L m of the leading-order Lagrangian L 0 in Eq. (2.12) are given by (2.14) Whether L m contributes at leading-order in λ still depends on the scaling of m as discussed further in the next section. Note that the factor in · D in Eq. (2.12) contains the collinear, the soft, and the Glauber gluon field The last term here gives the interaction vertex of a collinear quark with the medium off-shell Glauber gluons. It is exactly the same as in the massless case, see Eq. (2.7) above. As it turns out, there is no modification for this vertex due to the finite quark mass when we work in the hybrid gauge. As mentioned above, the Glauber gluon field scales as A G ∼ p + (λ 2 , λ 2 , λ 3 ). Following the usual power counting arguments, the only situation where we obtain such an interaction term is from Eq. (2.15). All other possible contributions are power corrections in λ to the leading-order Lagrangian. In particular, there is no Glauber gluon term in L m to leading-order in λ. In other words, the massive part of the vacuum SCET Lagrangian is not modified at leading-order when including medium interactions.
To summarize, we find that the Feynman rules that describe the interaction with the medium remain the same as in the massless case [35]. This statement holds in the hybrid gauge for a static source, which is the relevant case for all practical purposes considered in this work. All Feynman rules derived from the massive SCET M Lagrangian in the vacuum also remain the same. In the next sections, we are going to make use of them when calculating the massive splitting functions in vacuum as well as in the QCD medium. Since SCET M,G is a direct combination of SCET M in the vacuum and SCET G for the medium interactions, we do not repeat the relevant Feynman rules here but instead refer the reader to [35,48].

Massive splitting functions in the vacuum derived from SCET M
In this section, we derive the massive splitting functions in vacuum for the channels Q → Qg, Q → gQ, and g → QQ, using SCET M . The same results obtained in the standard perturbative QCD can be found in [61], where the authors derived their results using the so-called "quasi-collinear" limit. Such a limit is an extension of the "collinear limit" used for deriving massless splitting functions, where the on-shell quark masses are kept of the same order as the invariant mass of the two final state partons. See [61] for more details. Using SCET M , this limit has already been taken into account at the level of the SCET M Lagrangian when the mass is taken to scale as m ∼ p + λ. Therefore, the massive splitting functions can be obtained directly without having to make any further approximations.  Feynman diagrams for the two splitting processes involving massive quarks Q → Qg (left) and g → QQ (right). The off-diagonal splitting process Q → gQ can be obtained from the result for Q → Qg via crossing. J and J ν,b represents the remaining amplitude that produces the incident highly energetic parent quark or gluon respectively.
We adopt the notation introduced in [35,39]. We first consider the splitting process Q(p 0 ) → Q(p) + g(k), where we labeled the four momenta of the involved partons. The splitting process is illustrated in Fig. 1 (left). Throughout this paper, we adopt the convention to label a heavy quark of mass m by Q whereas a massless quark is denoted by q. As mentioned above, we work in the hybrid gauge for the in-medium splitting functions. Therefore, we adopt the light-cone gauge also for the vacuum calculation. Using the onshell conditions p 2 = m 2 and k 2 = 0 for the two outgoing partons as well as momentum conservation, we can parametrize the involved momenta as follows (2.16) Here, we choose to work in the frame in which the parent parton has no transverse momentum, x = k + /p + 0 is the momentum fraction taken away by the emitted gluon, and k ⊥ is the gluon momentum transverse to the parent parton momentum.
The amplitude A vac Q→Qg for the splitting process can be written as where J stands for the remaining amplitude producing the massive parent quark Q that undergoes the splitting process, see Fig. 1 (left). The factor R µ (p, k, m) is associated with the splitting vertex. Using the Feynman rules for SCET M , we find (2.18) The gluon polarization vector ε µ (k) in Eq. (2.17) can be written as [43,62] ε µ (k) = 0, satisfying both k · ε(k) = 0 andn · ε(k) = 0. Using the parametrization of the momenta as defined in (2.16), we can express the factor resulting from the massive quark propagator asn · (p + k) (2.20) Furthermore, we can write the product R(p, k, m) · ε(k) as Here, we separated the resulting three terms into two pieces where Γ i m contains only the third term proportional to the mass m in the first line of Eq. (2.21), since Γ i m has an odd number of γ ⊥ matrices. In summary, we can write the amplitude for the splitting process as To proceed we square the amplitude A vac Q→Qg and average over spin and color configurations of the initial quark: where we used Note that the second relation is the same for massless and massive collinear quarks. We can now write the expression in Eq. (2.23) involving the Γ i matrices inside the trace as 25) which is a scalar in both Dirac and color space. Here we used and (Σ 3 ) † = Σ 3 and (Σ 3 ) 2 = I. Hence, we may now write the splitting amplitude squared in a factorized form where the first factor is the leading-order amplitude squared, i.e. without any emission. In order to obtain the correct normalization for the splitting function, we need to factor out |A vac Q | 2 . Thus the splitting function for the process Q → Qg will be given by the second factor on the right-hand side of Eq. (2.27).
We still need to take into account the phase space for the quark (momentum p) and the gluon (momentum k) in the final state. Following for example [63], we write the corresponding 2-particle phase space as Here, the first factor is the leading-order phase space (momentum p 0 ) which needs to be factored out together with the leading-order amplitude squared |A vac Q | 2 in Eq. (2.27), in order to obtain the correct normalization for the splitting function. Eventually, the Q → Qg massive splitting function in the vacuum is given by which reduces to the massless splitting function derived analogously in [35,39] when taking m → 0. Different to the massless case, the dependence on x and k ⊥ does not factorize anymore. Likewise, we can derive the splitting kernel for the process Q → gQ. Being a crossed process of Q → Qg, the splitting function for Q → gQ can be obtained from Eq. (2.29) by substituting x → 1 − x. We now turn to the process g → QQ.

g → QQ
Next, we consider the splitting process where a gluon splits into a massive quark anti-quark pair: g(p 0 ) → Q(p) +Q(k), as shown in Fig. 1 (right). Analogous to Eq. (2.16), we may now write the momenta of the involved partons as The amplitude for the splitting process g → QQ can be written as The function J ν,b represents the remaining amplitude that produces the parent gluon, with Lorentz and color indices ν and b respectively. We assume again the physical polarization for the parent gluon, and thus J ν,b satisfiesn · J = (p + k) · J = 0, cf. Eq. (2.19). In addition, we have Note that the factor R µ (p, k, m) associated with the splitting vertex depends on the direction of the momentum flow. Analogous to Eq. (2.21), we can rewrite the combination Figure 2. Single-and double-Born Feynman diagrams that contribute to the massive in-medium splitting functions to first order in opacity. The topology is the same for all three splitting processes Q → Qg, Q → gQ and g → QQ. Figure adapted from [39].
In the second line, we separated again the m-dependent part from the rest as it contains only one γ ⊥ matrix. Squaring the amplitude A vac g→QQ in Eq. (2.31) and averaging over gluon polarizations and colors, we obtain the following result = δ ii δ bb which is the appropriate relation for the calculation of the spin averaged splitting function, see e.g. [64]. We further evaluate the expressions containing the Γ i 's and find We continue by evaluating the remaining trace part in Eq. (2.35) and by including the appropriate phase space factors as before. After taking into account the normalization to the leading-order amplitude squared, we obtain the following result for the g → QQ splitting function in the vacuum which reduces to the massless splitting function derived analogously in [39]. Again, the dependence on x and k ⊥ does not factorize as in the massless case.

Massive splitting functions in the medium derived from SCET M,G
The calculation of the massive in-medium splitting functions follows roughly the same steps as in the massless case [35]. For completeness, we outline the basic steps of the calculation. For every Glauber interaction of the energetic collinear parton with the i'th scattering center, we need to integrate over the Glauber gluon momentum. We denote the corresponding integral by dΦ i . Following [35], we introduce the following notation where q i is the momentum exchanged between the incident parton and the QCD medium (through the Glauber gluon exchange) and q i⊥ is its transverse component. Moreover, we have δx i = x i − x 0 , where x 0 is the space-time position where the initial energetic parton was created and x i is the position of the i'th interaction with the quasi-particles of the medium. The transverse component of δx i is denoted by δx i . The functions v(q i ) and v(q i⊥ ) are related to the elastic scattering cross section. We have v(q i ) = 2πδ(q 0 i )ṽ(q i⊥ ) and to lowest order for a Yukawa-screened potential where C 2 (R) and C 2 (T ) denote the quadratic Casimir invariants in the representation of the incident highly energetic parton and the target (source) respectively. See [35] for more details. The delta function in v(q i ) makes the q + integral trivial. We can now write dΦ i as where δz i is the distance along the z-axis between the scattering center i and the point where initial energetic parton was created. Remaining integrals over q ⊥ will eventually be performed numerically at the end using a realistic model for the medium. However, the q − i integrals still need to be performed analytically. The corresponding longitudinal integrals can be evaluated in terms of contour integrals in the complex plane. The results including non-vanishing mass terms can be obtained analogously to the techniques outlined in [35] for the massless case. Formally, the amplitudes for the three in-medium splitting processes, can be written as where we use the same labelling of the involved parton momenta as for the vacuum case in the previous section. Here, Q → ab corresponds to either the diagonal splitting Q → Qg or the off-diagonal case Q → gQ. The gauge invariant quark and gluon fieldsχ n and B µ n⊥ were defined in (2.4) and L SCET M,G was derived in section 2.2 above. In order to obtain the in-medium splitting functions to first order in opacity, we need to take into account both single-and double-Born diagrams for every scattering center i, see [43]. We denote the corresponding single-and double Born amplitudes by A med SB and A med DB respectively. The double-Born diagrams are evaluated in the "contact limit" where δz 1 = δz 2 . After squaring the sum of all amplitudes, we have to calculate schematically where A vac is the vacuum splitting amplitude without any interaction with the QCD medium. In Fig. 2, all the relevant diagrams are shown that correspond to Eq. (2.43). The topology is the same for all three massive splitting processes. On the left hand side of Fig. 2, the square of the three single-Born diagrams is shown. On the right hand side, all double-Born diagrams are shown that give a non-zero result in the contact limit. For all double-Born diagrams, the interference with the vacuum leading-order splitting diagram is shown. After adding and squaring the relevant amplitudes, we still need to perform the sum over all scattering centers i = 1, . . . N . Following [35], this sum can be turned into a continuous integral that gives a delta function, which in turn can be used to perform one of the remaining transverse momentum integrals. We would like to stress that to first order in opacity, we take into account the single-and double-Born diagrams shown in Fig. 2 for every scattering center i. See [35] for more calculational details.
Since the intermediate steps of the calculations are very similar to the massless case [35], we will skip them and present only the final results. On the other hand, for comparison and for later convenience, we also list the results for the massless in-medium splitting functions as calculated in [35,39]. We define the transverse momentum vectors where again x and k ⊥ are the longitudinal momentum fraction and the transverse momentum of the emitted parton relative to the parent parton respectively. Furthermore, q ⊥ is the transverse momentum introduced by the Glauber gluon exchange. In addition, we have the following phases , , .

(2.45)
We reproduce the light parton in-medium splitting functions for reference and subsequent comparison. The result for q → qg is given by (2.46) Here λ g,q (z) is the gluon (quark) mean free path in the medium. (1/σ el )dσ med el /d 2 q ⊥ is the normalized in-medium elastic scattering cross section, see (2.39) above. We are left with a three dimensional integral over q ⊥ and ∆z that need to be evaluated using a realistic model for the QCD medium as it is produced in heavy-ion collisions. The latter is an integral over the interactions with the medium quasi-particles 0 < ∆z < L, where L is the size of the medium. Note that the in-medium splitting function for q → gq can be obtained via crossing x ↔ 1 − x. The q ⊥ integral will be evaluated numerically and is subject to phase space constraints. The results for the in-medium splitting processes g → gg and g → qq are given by Note that for all four massless splitting functions, the x-dependent vacuum splitting function factors out. Next, we present our final results for the massive in-medium splitting functions Q → Qg, Q → gQ and g → QQ. All transverse momentum vectors defined in Eq. (2.44) remain the same, so are the phase factors Ω 2 − Ω 3 and Ω 5 in Eq. (2.45). However, the remaining three phases in Eq. (2.45) are modified as follows , , where the variable ν is given by for the three different massive splitting processes, respectively. The full in-medium splitting function for Q → Qg is given by where the ellipses denote analogous terms as in the first square bracket following the pattern as indicated. The variable ν for the process Q → Qg was defined in Eq. (2.48), i.e. ν = x m. The expressions in both square brackets have the same structure as the full massless inmedium splitting functions. The result for Q → gQ can be obtained via crossing. The splitting function g → QQ is given by Here, the ellipses denote again analogous terms as in the first square bracket following the pattern as indicated, and ν = m as given in Eq. (2.50). Note that in all three cases, the massive vacuum splitting functions given in Eqs. (2.29) and (2.37) do not factor out as it was the case for the massless in-medium splitting functions.

Soft gluon approximation
In this section, we consider the soft gluon approximation (SGA) of the full massive splitting functions in Eqs. (2.51) and (2.52) by taking the limit x → 0. We then make the connection with the results obtained in the traditional picture of parton energy loss [18]. For comparison, we first present the massless results in the SGA as derived in [35,39] and earlier in [43] using the traditional approach of parton energy loss: (2.53) Note that we keep the O(x) expressions for the off-diagonal splitting functions even though they actually vanish for strictly x → 0. Next, we consider the massive splitting functions in the SGA. We would like to point out that there is some ambiguity for defining the massive small-x result. This ambiguity arises for the diagonal splitting process Q → Qg, where the mass term is proportional to x 2 m 2 . This means that any mass dependence vanishes for strictly x → 0. Therefore, in order to keep a finite mass correction even for x → 0, one conventionally chooses to keep the first order correction of the mass in the denominator. However, this convention leaves some ambiguity at what stage of the derivation one should keep the first order correction in the denominator. When deriving the SGA, one ends up with the following expression where we kept the mass terms ∼ x 2 m 2 in the denominators. This structure is eventually multiplied by a factor involving a cosine similar to the massless case in Eq. (2.53). Keeping the masses at this stage would be consistent with the convention in [18], where the massive small-x result was derived within the conventional approach to parton energy loss. In [18], the final result is cast in the following form (2.55) Note that this result involves three factors in the denominator which is different than in the massless case. However, this expression can also be written as where now the first term has a similar structure as the massless result in Eq. (2.53). The second term is proportional ∼ x 2 m 2 and vanishes for strictly x → 0. Therefore, we choose to keep only the first term in Eq. (2.56) in the SGA. This version of the massive SGA for Q → Qg is more similar to the massless result and more importantly, it is consistent with the results for the off-diagonal splitting functions in the SGA where there are no ambiguities.
Note that for g → QQ, the first mass correction has no x-dependence, see Eq. (2.52) and also Eq. (2.59) below. For Q → gQ, the mass correction is proportional ∼ (1−x) 2 m 2 which also becomes ∼ m 2 for x → 0. Since there is no ambiguities for the off-diagonal small-x results, we choose to define the massive diagonal SGA result for Q → Qg in analogy to them. The complete result for Q → Qg in the SGA is (2.57) Note that we also keep a finite mass correction in the phase of the cosine. We would like to stress again that this convention is different than the one chosen in [18] where the structure in Eq. (2.54) was used instead. We now continue presenting the results for the off-diagonal splitting processes. In the small-x limit, we have the following result for Q → gQ and for g → QQ Note that for the two off-diagonal splitting functions, the mass correction is directly ∼ m 2 without any dependence on x as discussed above.

Numerical results
In this section, we present numerical results for the massive in-medium splitting kernels. We compare the full splitting kernels with the soft gluon approximated results and study the finite mass effects. We perform the q ⊥ and ∆z integrations in Eqs. For comparison, we also plot the massless results q → qg for both the full splitting function (dashed black) and the small-x limit (green). The q ⊥ and ∆z integrals are evaluated numerically with a realistic model for the medium and physical phase space cuts, see text and [45]. As an example, we choose the incident parent parton energy as E 0 = p + 0 /2 = 20 GeV (left) and E 0 = 100 GeV (right).  In Figs. 3 and 4, we show the results for the intensity spectra x(dN/dx), which are obtained by integrating over k ⊥ up to k ⊥,max = 2E 0 x(1 − x). This choice is just for illustrational purposes. In section 3 below, we consider a different upper integration limit for k ⊥ that is required by how the in-medium splitting functions have to be treated when they appear in an actual cross section. In Fig. 3, we consider two initial parton energies E 0 = p + 0 /2 = 20 GeV (left) and E 0 = 100 GeV (right). We show the full massive splitting function for Q → Qg in blue, whereas the corresponding soft gluon approximated result is shown in red. For comparison, we also plot the massless results q → qg for both the full splitting function (dashed black) and the small-x limit (green). Note that we take into account a finite Debye mass m D = gT for the medium, which appear for both the massless and the massive cases. It can be seen clearly that the mass effects are a large-x effect since all four curves for both choices of E 0 in Fig. 3 are very close together at small-x. This is to be expected as the mass corrections are of the form ∼ x 2 m 2 for Q → Qg. By comparing the two results for the full splitting functions for massive quarks (blue) and massless quarks (dashed black), one finds a significant difference in the large-x region for x > 0.4. Interestingly, the rise of the massless result at large-x completely disappears when considering a finite bottom quark mass. As expected, the finite mass results are more relevant for E 0 = 20 GeV (left), where the differences are clearly larger.
In Fig. 4, we present analogous numerical results for the off-diagonal splitting functions Q → gQ (left) and g → QQ (right) for E 0 = 100 GeV. The finite mass effects are even more pronounced here than for the diagonal splitting Q → Qg, and can be relevant for both the large and the small-x region. The enhanced effect in the small-x region is consistent with the fact that the mass corrections for the processes Q → gQ and g → QQ are proportional ∼ (1 − x) 2 m 2 and ∼ m 2 respectively, and thus remain finite when taking x → 0. Although the mass corrections can be large in the small-x region, it is instructive to keep in mind that both off-diagonal splitting functions vanish when x → 0, as can be seen clearly from Eqs. (2.58) and (2.59). Therefore, the overall numerical impact of the finite mass effects at the level of x(dN/dx) from these regions is not directly translated to the cross section in heavy-ion collisions.
In summary, we find that finite mass effects are indeed very significant at the level of the intensity spectra x(dN/dx). Eventually this can have a sizable numerical impact for the suppression of heavy mesons in heavy ion collisions as discussed in the next section.

Application to PbPb → HX at NLO
In this section, we first introduce a new framework for including in-medium effects consistent with next-to-leading order calculations in QCD for inclusive hadron production in heavy ion collisions. This can be achieved by making use of the in-medium massive slitting functions derived in last section and by effectively introducing in-medium fragmentation functions. We then consider the cross section for open heavy flavor production in protonproton collisions, and provide numerical results in the so-called zero mass variable flavor number scheme (ZM-VFNS).

In-medium fragmentation functions
In this section we derive a framework to include in-medium interactions for PbPb → HX which is consistent with NLO calculations in the vacuum for pp → HX. Initial state Cold Nuclear Matter (CNM) effects will be included only for numerical evaluations at the very end, with the actual implementation explained in more details in [45]. An improved treatment of CNM energy loss using SCET G -based initial-state splitting functions [41] will be left for future work. The framework that we develop in this section is related to jet calculations in [46,47] and to some extend it corresponds to a first order expansion of the DGLAP formalism developed in [44,45]. For a discussion of a medium-modified DGLAP in semi-inclusive DIS, see [65].
Interactions with the hot and dense QCD medium affect partons after the hardscattering event but before they eventually fragment into hadrons. To NLO in the strong coupling constant, we have to consider one-loop real and virtual corrections for the outgoing final state parton. As an example, we consider the corrections to the leading-order hard process qq → qq as shown in Fig. 5. In the vacuum, a splitting process such as that shown in Fig. 5 (left) needs to be taken into account. Such a contribution will eventually lead to the DGLAP evolution of the vacuum fragmentation function. On the other hand, in the QCD medium, besides the vacuum splitting process, an medium-induced splitting process as shown in Fig. 5 (right) will also happen, which leads to additional contributions to the cross section for the hadron production in heavy ion collisions. Of course, when squaring the amplitude corresponding to the medium-induced diagram on the right hand side of Fig. 5, we actually need take into account all relevant single-and double-Born diagrams to first order in opacity as discussed in section 2. The gray ellipses denote the standard vacuum fragmentation function for both situations. We start by rederiving the vacuum case and we then continue by describing how this calculation can be extended to the medium case. At one-loop order, the relevant part that describes the splitting of the final state parton can be schematically written as Here,σ i is the leading-order hard-scattering cross section to produce a parton i, P ji is the leading-order Altarelli-Parisi splitting function for i → j and D H j is the parton-to-hadron fragmentation function. This generic structure in Eq. (3.1) can be found in standard textbooks such as [66]. The symbols ⊗ denote convolution products. This structure can be obtained by calculating a parton-to-parton fragmentation function to one-loop order or by considering the relevant splitting process to be part of the hard-scattering function, see [67,68]. Note that these two possibilities are fully equivalent to first order. We choose to present the calculation for a parton-to-parton fragmentation function. Conceptually, we want to treat the two splitting processes shown in Fig. 5 (vacuum and medium case) to be part of the first order correction to the leading-order process qq → qq.
We start by calculating the massless partonic quark and gluon fragmentation functions in the vacuum. Massive in-medium splitting functions can be implemented in a straightforward way as well. We will comment on the extension to massive quarks below. For q → q and g → g, we need to take into account both real and virtual corrections. Using the method of [69], one can express the contributions of the virtual graphs in terms of splitting functions derived from real emission graphs. This is consistent with the so-called flavor and momentum sum rules [70]. From here on, we switch to the more traditional convention, where for any given splitting process, the radiated parton carries a momentum fraction 1 − z instead of z as in the previous section. For q → q, we have where we use the notation k ⊥ = |k ⊥ |, and k ⊥ is integrated between a lower scale Q 0 and an upper cutoff µ that is usually identified as the relevant hard scale of the process in consideration. If one takes the derivative of Eq. (3.2) with respect to the upper integration limit µ, one will derive the DGLAP evolution equations. Note that the integral over x in the first line is divergent by itself but it is cancelled between real and virtual contributions, and we are left with a regularized plus distribution in the second line. For the two off-diagonal fragmentation functions at one-loop order, we have The process g → g is slightly more involved and also needs special attention in the medium case For notational convenience and easy generalization to the medium case, let us introduce the functions P i→jk (z, µ) for every splitting process i → jk, where k corresponds to the emitted parton carrying away the momentum fraction 1 − z. The functions P i→jk (z, µ) are related to the splitting functions (dN/dz/d 2 k ⊥ ) i→jk defined in section 2 as We use this identification both for the vacuum and the medium case. To be specific, we always include a superscript "vac" or "med" below. Note that we included the k ⊥ integral in the definition of P i→jk (z, µ) as it is always the same. For example, for the splitting process q → qg in the vacuum, we have (3.7) Using this notation, we can now write the partonic vacuum fragmentation functions derived above as which contain both real and virtual contributions. With such notations, it is straightforward to extend these results to the medium case, where one has to consider both the vacuum splitting function as well as a medium induced part. Therefore, we can directly make the following substitutions in Eq. (3.8) above Here, the P vac i→jk (z, µ) are the vacuum splitting functions as discussed above and P med i→jk (z, µ) are the in-medium splitting functions as derived in section 2 and integrated over k ⊥ as defined in Eq. (3.6). Eventually, we need to "match" onto the standard vacuum fragmentation functions. We now continue by evaluating this matching procedure and calculate the convolution with the standard fragmentation functions. Schematically, we have the following structure for the medium case Effectively, the functions D H,med i (z, µ) can be considered as medium-modified fragmentation functions [71,72]. Even though we consider the medium induced splittings as a correction to the vacuum hard-scattering function, it is notationally more convenient to think of it as a medium-modified fragmentation function. The medium-modified FF will then be convolved with the leading-order hard-scattering cross section. As it is formally a one-loop correction, we are then going to add it to the NLO calculation in the vacuum.
Following the definition of the medium-modified quark and gluon fragmentation functions D H,med i in Eq. (3.10), we find At this point one can make a direct connection between the new treatment of the medium effects as proposed here and the approach considered in [44,45]. In these two papers, the authors derived medium-modified DGLAP equations. The DGLAP equations including medium effects can be obtained by taking the derivative of the above Eq. (3.11) with respect to the scale µ. Using medium-modified DGLAP equations essentially leads to an exponentiation of the in-medium branchings. In [44,45], a close connection was established between the medium-modified DGLAP equations and traditional parton energy loss calculations. In our case, we only consider the first order correction in α s evaluated in the opacity series expansion. The numerical results of the two approaches turn out to be very similar but the approach proposed in this paper is easier to implement. In addition, the new approach has a close connection to NLO calculations in the vacuum which we use as proton-proton baseline as discussed in the next section. Differences between the two approaches are expected only when the observed hadrons have a relatively small transverse momentum p T . Although schematically correct and overall finite, Eq. (3.11) cannot be used as they are since the individual terms can become numerically divergent at the phase space boundaries z → 0, 1 due to the behavior of the splitting functions P i→jk (z, µ), see Eqs. (2.46), (2.47), and (2.51). To rewrite the expressions into separate "numerically stable" pieces, we introduce plus distributions similar to [45]. For example, the quark in-medium fragmentation function can be directly written as We can now write the sum of both real-and virtual gluon contributions to the g → g splitting function as where we used h med g→gg (1 − x, µ) = h med g→gg (x, µ). We choose a slightly different convention for expressing the plus distribution as in [44,45]. All conventions are equivalent as long as the divergence in the x integral is cancelled. The version here is the minimally required one, where only h med g→gg (z, µ)/(1 − z) are written in terms of a plus distribution. Including the off-diagonal contribution, we can now write the in-medium fragmentation function D H,med where dσ H,NLO pp is the NLO cross section in the vacuum, and dσ H,med PbPb is the one-loop medium correction, and schematically we have (3.17)

Numerical results for pp → HX at NLO within the ZM-VFNS
In this section we present numerical calculations for open heavy meson production in proton-proton collisions, pp → HX. We choose to work within the ZM-VFNS [73,74], in which one neglects all the heavy quark mass corrections in the partonic hard-scattering functions [75,76]. Thus, this approximation is justified for µ m c,b , where µ is the characteristic scale of the process and m c,b are the charm and bottom quark masses respectively. In the vacuum, the only characteristic scale of the process is the large transverse momentum of the produced hadron µ = p T . Therefore, the ZM-VFNS is expected to be applicable in the high-p T range: p T m c,b . The exact range of validity of the ZM-VFNS needs to be checked by comparing theory and experimental data. For this reason, we perform several exemplary numerical calculations for pp → HX below.
On the other hand, the medium contribution is also sensitive to the properties of the QCD medium which introduces much lower energy scales than p T . The characteristic scale for in-medium interactions is given by the typical momentum exchange between the incident parton and the QCD medium, which can be even smaller than the heavy quark mass. Therefore, while one can set the heavy quark masses to zero in the hard-scattering functions because of p T m c,b , we do take into account the masses of both charm and bottom quarks in the medium-modified FFs for studying the medium contribution as discussed in the previous section 2. Such a set-up for the medium contribution is similar to [77].
We use the pp → HX NLO framework developed in [67,68] and typically applied for the production of light hadrons [75,[81][82][83][84][85][86][87][88]. The double differential cross section can be written in the following way where a,b,c stands for a sum over all the parton flavors including light and heavy quarks and gluons, and s, p T and η correspond to the center of mass energy, the hadron transverse momentum and hadron rapidity, respectively. Moreover, f a,b (x a,b , µ) are the parton distribution functions for the two incoming protons. The hard functions dσ c ab (ŝ,p T ,η, µ) are functions of the corresponding variables at the parton level: the partonic CM energŷ s = x a x b s, the partonic transverse momentump T = p T /z c and the partonic rapidity The integration limits in (3.18) are customarily written in terms of the hadronic variables V, Z, 21) and are given by On the other hand, D H c (z c , µ) are the heavy meson fragmentation functions. For charmed mesons, we use the fragmentation functions of [53,[78][79][80], whereas for B-mesons we use the ones from [54,89,90]. 10 100 √ s = 7 TeV, pb/GeV of these terms is small for sufficiently large p T . See for example [49-52, 55, 56, 91-95] for other sets and possible approaches to heavy meson fragmentation functions.
We start by presenting a comparison of NLO results in the ZM-VFNS with D-meson data taken by the ATLAS and CMS collaborations in Fig. 6. In Fig. 6 (left), we show the ZM-VFNS results for D * ± mesons at a CM energy of √ s = 7 TeV as a function of the transverse momentum. The rapidity is integrated over an interval of |η| < 2.1. Note that here D * ± does not correspond to the average but to the sum D * ± = D * + + D * − . Following [53,[78][79][80], we choose µ R,F = m T = p 2 T + m 2 c as the central values for the renormalization and factorization scales. The band is obtained by varying µ R,F independently around their central values by a factor of 2 and by taking the envelope. Along with the theoretical calculation, we show the ATLAS data of [11]. Throughout this section, we always show the combined statistical and systematic errors added in quadrature. Analogously, in Fig. 6 (right), we show the results for D 0 production at √ s = 5.02 TeV with |η| < 1 comparing to preliminary data from CMS [13]. Keeping in mind that the D-meson fragmentation functions of [53,[78][79][80] are fitted to e + e − data only, we find that the agreement between theory and data is indeed remarkably good. In addition, we would like to emphasize that the agreement between the NLO calculation within the ZM-VFNS and the data is still good even at relatively low values of p T of a few GeV. Similarly, in Fig. 7, we show analogous comparisons for B-mesons. We choose to only show two exemplary comparisons of the NLO calculation in the ZM-VFNS and inclusive pp → B + X data from ATLAS [10] (left) and CMS [7] (right). For both data sets, B + mesons are identified via the exclusive decay channel B + → J/ψK + → µ + µ − K + . The corresponding multiplicative branching fractions are taken into account in our calculations.
The CMS data is integrated over |η| < 2.4, whereas the ATLAS data is presented for four different rapidity intervals in the range of η = 0 − 2.25. Both data sets were taken at a CM energy of √ s = 7 TeV. The default scale for the NLO calculation is now µ R,F = m T = p 2 T + m 2 b following [54,89,90]. Again, the band is obtained by varying µ R,F independently around their default choice by a factor of 2 and by taking the envelope. Similar to the inclusive D-meson production, the agreement between theory and data is remarkably good even down to relatively low p T .
We would like to point out an important difference to several earlier calculations in the literature. Often the heavy meson production is calculated differently and only the modification of the heavy-quark-to-heavy-meson fragmentation process is taken into account in Pb+Pb collisions. However, in the ZM-VFNS, there is a large gluon-to-heavy-meson contribution, even though the gluon-to-heavy-meson fragmentation function itself is much smaller than the corresponding heavy-quark-to-heavy-meson FF. The smallness of the gluon FF itself is compensated by the large gluon production cross section in proton-proton collisions at high CM energies. In fact, as illustrated in Fig. 8, the gluon-to-heavy-meson contribution (shown in blue) is of the order of 50% for both D-(left) and B-meson (right) production in pp collisions at √ s = 7 TeV. The percentage contribution of the heavy-quark-to-heavymeson fragmentation process is shown in red for both D-and B-meson production. We note that the light-quark-to-heavy-meson fragmentation contribution turns out to have a marginal effect only. We would like to stress again that heavy meson FFs are generally extracted from e + e − data only. In this case, the gluon-to-heavy-meson FF only enters at the one-loop level and through evolution effects. This leads to the fact that the gluon fragmentation function is relatively poorly constrained from e + e − alone. In the future, it will be very helpful to obtain heavy meson fragmentation functions within a global analysis including in particular pp → HX data as it is customarily done for light hadrons [82,83,85].
In addition, including in-jet fragmentation data pp → (jetH)X [96], an observable for which a new theory framework was recently developed [97][98][99][100][101][102][103][104][105], is expected to lead to great improvements of the corresponding global fits. Whether a gluon-to-heavy-meson fragmentation function is included in the calculation or not is especially relevant for the in-medium calculation. The energy loss of heavy quarks and gluons in the medium is very different. Therefore, it is absolutely crucial to understand how heavy mesons are formed in order to obtain a reliable quantitative understanding of their suppression in heavy-ion collisions. In fact, not only the relative percentage of gluon and heavy quark fragmentation as shown in Fig. 8 is important but also the exact shape of the fragmentation functions is relevant. The inclusive hadron spectra pp → HX are relatively well described by currently available sets of fragmentation functions as shown above. However, for example, the disagreement between theory and data for heavy mesons measured inside jets clearly shows that the currently available fragmentation functions are still not well enough understood so far, see [101]. Performing new global fits is beyond the scope of this work but we are planning to addressed this issue in future publications.

Suppression of Dand B-mesons in Pb+Pb collisions at the LHC
In this section, we present the results for the nuclear modification factor R AA defined as Here, N bin is the average number of binary nucleon-nucleon collisions for a given centrality and dσ H pp /dηdp T , dσ H PbPb /dηdp T are the double differential cross sections for inclusive heavy meson H production in proton-proton and Pb+Pb collisions respectively. The cross section for proton-proton collision was given in Eq. (3.18) and its modification for Pb+Pb collisions was discussed in the previous section, cf. Eq. (3.16). Our calculations depend on one parameter and the result of initial-state effects. Firstly, there is the coupling constant g that describes how strongly the hard partons couple to the QCD medium. As in several earlier publications [44,45], we choose this parameter around g ≈ 2. When presenting numerical results for the nuclear modification factor R AA , we typically vary this parameter around its central value by ±0.1 and plot the obtained band. Eventually, the coupling strength g will have to be constrained by comparing to data. The second set of effects are Cold Nuclear Matter (CNM) effects, which happen before the formation of the QGP. These include isospin effects, coherent power corrections for heavy quarks [22,106], the Cronin effect [107] and cold nuclear matter energy loss [41,108,109]. By implementing isospin effects we take into account that the Pb nucleus is made up of both protons and neutrons. As discussed above and illustrated in Fig. 8, we find that for heavy meson production only heavy quark and gluon fragmentation functions turn out to be relevant. These processes are isospin symmetric. Therefore, for all our results presented here, the isospin effects are very small. The effect of power corrections is limited to small p T and does not affect the ZM-VFNS region of applicability. The Cronin effect and CNM energy loss effects can partly be constrained by p+Pb collisions for example. It is clear that there is always a non-trivial interplay between the value of the coupling strength g and possible CNM effects. To some extent a variation of, say, g can be absorbed by changing the strength of CNM effects. That being said, we would like to point out that the two effects are not completely interchangeable and it is indeed possible to constrain both effects from precision data in Pb+Pb collisions when observables are carefully selected, see also [45].
We start by presenting results for the R AA of D 0 -mesons using the new framework of SCET M,G (hatched red band) in Fig. 9. In addition, we present results obtained within the traditional approach to parton energy loss (green band). We choose a CM energy of √ s NN = 5.02 TeV and integrate over the rapidity interval |η| < 1. The charm mass is chosen as m c = 1.3 GeV. The results are presented for central collisions with centrality 0 − 10%. The bands are obtained by varying the coupling strength around its central value g = 1.9 ± 0.1. On the left (right) hand side of Fig. 9, we show the results without (with) CNM effects. It can be seen that the CNM effects can affect the R AA both at low and highp T . As it turns out, they lead to a rise at relatively low-p T , whereas a suppression in the high-p T region is observed. We find that both the SCET M,G results and the results based on traditional parton energy loss are quite similar in the large-p T region. However, at low-p T , the two results differ significantly. We would like to point out that the difference between the two results is not entirely due to the different theoretical approaches to the in-medium interaction SCET M,G vs. traditional parton energy loss. For both D-and B-meson (see Fig. 10) production within SCET M,G , we use modern fits of fragmentation functions that include both heavy quark and gluon fragmentation. Instead, the results presented here using the traditional picture of parton energy loss are calculated using only heavy quark fragmentation functions based on a model calculation [94], as it is conventionally done in the literature, see e.g. [110]. The different choice of fragmentation functions can lead to a very different result for the R AA in particular at low p T . We discuss this point in more detail below. In Fig. 10, we present analogous results for the nuclear modification factor for B + meson production. The bottom mass is chosen as m b = 4.5 GeV. By comparing Figs. 9 and 10 one notices that there is indeed a difference of the R AA between D 0 -and B + -meson suppression independent of the approach (medium-induced splitting only). Firstly, this is in part due to the different fragmentation functions. Secondly, the different masses for charm and bottom quarks can affect the R AA even at relatively large p T . As can be seen from Fig. 10, the difference between SCET M,G based results and the results from traditional parton energy loss differ more significantly at low-p T than it is the case for D 0 -mesons. The large difference between the two approaches at low p T is mainly due to the fact that for the traditional parton energy loss calculation we only take into account heavy quark fragmentation functions.
As already pointed out in the previous section, it is of great relevance to understand the relative contributions of heavy quark and gluon fragmentation to the heavy meson production cross sections [24,29]. Fig. 11 illustrates the implications for the obtained nuclear modification factor R AA . In black, we show the combined calculation for the suppression of D 0 (left) and B + mesons (right) in heavy-ion collisions as before in Figs. 9 and 10. In blue, the suppression is shown for the "heavy quark only" case. These results are  Figure 11. Nuclear modification factor R AA for D 0 mesons (left) and B-mesons (right). As an example, we choose √ s NN = 5.02 TeV, |η| < 1 and 0 − 10% centrality. In black the standard R AA is shown as in Figs. 9 and 10. In blue, we show the suppression for the "heavy quark only" case. This result is obtained by calculating both the proton-proton baseline as well as the medium effects with only charm (left) and bottom quark fragmentation functions (right). In red, the analogous result is shown for the "gluon only" case.
obtained by calculating both the proton-proton baseline as well as the in-medium effects with only charm quark (left) and bottom quark (right) fragmentation functions. In red, we show the analogous result for the "gluon only" case. While the suppression for heavy quarks and gluons is similar in the high p T region, it turns out that their suppression is very different in the low p T region. The difference between fragmenting gluons and heavy quarks is more pronounced for the heavier B + -mesons. The very different suppression rates for heavy quarks and gluons can lead to a significantly different picture of how the QCD medium affects open heavy flavor. Besides these important differences, there are two main sources of uncertainties at low p T . Firstly, the gluon-to-heavy-meson fragmentation function is still relatively poorly constrained. This concerns both the exact functional form as well as the total contribution to the cross section of gluons as discussed above and illustrated in Fig. 8. Therefore, we would like to stress that a more reliable picture of the in-medium interactions requires further improvements already at the level of the protonproton baseline calculation. Secondly, in the low p T region higher order terms in the opacity series expansion are expected to play a more important role. In the future we plan to address these issues in order to systematically improve the current framework allowing an extension to lower values of p T . In addition, other effects like collisional energy loss and dissociation are expected to play a role at low p T as well [20,21]. However, given the currently remaining uncertainties at low p T , it is clear that further improvements are needed before making any definitive statements about where exactly other effects are relevant or even dominate. We note that our conclusions here are different than in [32], where the authors concluded that the energy loss for gluons that fragment into heavy mesons is small because gluons quickly split into heavy quark anti-quark pairs which then fragment into the observed heavy mesons. Instead, here we are motivated by QCD factorization saying that the actual hadronization is a long-distance effect which happens at much later time scales than the hard-scattering event. At least within the ZM-VFNS, gluon fragmentation has been put on an equal footing as the heavy-quark fragmentation function within the QCD factorization formalism, as can be seen clearly in Eqs. (3.16) and (3.18). In this sense, our approach is in direct analogy to light charged hadron production. We would like to add that a unique opportunity to test and improve the current theoretical framework would be to measure and calculate the in-jet fragmentation of heavy mesons both in proton-proton and heavy-ion collisions [96,101]. Firstly, the poorly known gluon-to-heavy-meson fragmentation functions can be studied at a more differential level. Secondly, the in-jet fragmentation will allow to disentangle better the modification of the two main fragmentation contributions to heavy meson production.
Despite the remaining uncertainties at low p T , we expect to have a reliable description of the in-medium effects as long as the transverse momentum of the observed heavy meson is sufficiently large p T 10 GeV. We proceed by comparing our new SCET M,G based calculations with currently available experimental data from CMS and ALICE as an example. In Fig. 12, we present a comparison to the preliminary CMS data of [13] for the nuclear modification factor R AA for D 0 -mesons. The data was taken for √ s NN = 5.02 TeV, |η| < 1 and 0 − 10% centrality. Note that for all the data presented in this section, we show the statistical errors as standard error bars and the systematic ones as yellow boxes. On the left hand side, we compare the data to our calculation without CNM effects and we choose the coupling strength as g = 2.0 ± 0.1. Instead, on the right hand side, we include  Figure 13. Nuclear modification factor R AA for D-mesons (D 0 , D + and D * + average) without (left) and with (right) CNM effects in comparison to preliminary ALICE data of [9]. We have √ s NN = 2.76 TeV, |η| < 0.5 and 0 − 7.5% centrality. When CNM effects are (not) included, we choose the coupling strength as g = 1.9 ± 0.1 (g = 2.0 ± 0.1).
CNM effects as discussed above. As illustrated in Fig. 12, it turns out that CNM effects lead to a larger suppression at high-p T . Therefore, we choose a lower coupling strength g = 1.9 ± 0.1 for the comparison to data. Again, we would like to point out that one effect can not be compensate entirely by the other one. Both calculations based on the newly derived SCET M,G agree very well with the data in the expected p T region. While the result without CNM effects seems to agree slightly better with the data, one can not make any definitive statement given the experimental uncertainties. Next, we compare to ALICE data [9] at √ s NN = 2.76 TeV in Fig. 13. The nuclear modification factor R AA is shown for D-mesons (D 0 , D + and D * + average) with |η| < 0.5 and 0 − 7.5% centrality. Again we present our numerical results including CNM effects for g = 2.0 ± 0.1 (left) and the results without CNM effects for g = 1.9 ± 0.1 (right). We find that the data is well described by our calculations for p T 10 GeV. Our result without CNM effects seems to agree slightly better with the data.
Finally, in Fig. 14 we compare to the R AA data from CMS [8] for non-promt J/ψ production which originate from the decay of B-mesons. The data was taken at √ s NN = 2.76 TeV for |η| < 2.4 and is available only for minimum bias collisions (0−100% centrality). For comparison, we show our results for 0-10% centrality and for mid-peripheral collisions with 30-50% centrality. Again, we present our theoretical B-meson results without CNM effects (left, g = 2.0 ± 0.1) and with CNM effects (right, g = 1.9 ± 0.1). We find that our central results (0 − 10% centrality) agrees very well with the data. The minimum bias results are dominated by central collisions as they are weighted with the number of collisions. It will be instructive if the experiments can provide these data sets for fixed centrality bins, and thus to further test our theoretical framework. We look forward to  Figure 14. Nuclear modification factor R AA for non-prompt J/ψ production which originate from B-meson decays. The CMS data [8] was taken at √ s NN = 2.76 TeV for |η| < 2.4 and 0 − 100% centrality. For comparison, we present our B-meson results for central (0 − 10% centrality, hatched red band) and mid-peripheral (30 − 50% centrality, green band) collisions. Again, we present our theoretical results without CNM effects (left, g = 2.0 ± 0.1) and with CNM effects (right, g = 1.9 ± 0.1).
more experimental data at the LHC in the near future.

Conclusions and outlook
We have derived a version of Soft Collinear Effective Theory that includes both the interactions with the medium that are mediated by Glauber gluon exchange and heavy quark masses. Using the new effective field theory, we obtained vacuum and in-medium massive splitting functions for the Q → Qg, Q → gQ and g → QQ processes. Despite some ambiguities, we found agreement in the soft emission limit with earlier results in the literature where traditional approaches to parton energy loss in the QCD medium were used. In addition, we proposed a new formalism to include in-medium effects consistently at next-to-leading order in QCD. We presented numerical open heavy flavor results for proton-proton collisions in the ZM-VFNS. Comparing with currently available data, we found good agreement even for relatively low p T . Our numerical results for the suppression of open heavy flavor production in Pb+Pb collisions are applicable for p T 10 GeV. We observed good agreement between theory and currently existing data sets for both D-and B-meson production at √ s NN = 5.02 TeV and 2.76 TeV. As it turns out, the low-p T suppression rates for open heavy flavor production are very sensitive to the relative contributions of heavy quark and gluon fragmentation. The currently available sets of fragmentation functions may not be sufficiently well constrained to make quantitative predictions in the very low p T region. In the future, there are several possible ways to improve the current framework. Firstly, our study clearly motivates global fits of heavy meson fragmentation functions which are currently only constrained from e + e − data alone. Including both pp → HX and hadron-in-jet pp → (jetH)X data will greatly improve the sensitivity in particular to the gluon-to-heavy-meson fragmentation functions. Secondly, it would be interesting to calculate the full in-medium splitting functions to second order in opacity. This way, it may be possible to extend the current framework down to lower p T , where correlated multiple interactions with the medium become more relevant. This way, the full range of applicability of the current framework can be assessed and additional effects in the medium can also be taken into account.