We start by outlining the main features of the interactions between open-heavy flavor mesons, i.e., mesons with one charm or bottom quark H ∋ { D , D s , D ∗ , D s ∗ , B ̄ , B ̄ s , B ̄ ∗ , B ̄ s ∗ } , and the light Goldstone bosons Φ ∋ { π , K , K ̄ , η } within the framework of HMChPT, an effective field theory that has been widely used in the last years to describe the interactions between open-heavy flavor mesons and light mesons [6, 43–46], with chiral symmetry and heavy-quark spin-flavor symmetry (HQSFS) as its guiding principles. When combined with a unitarization technique such as the solution of the Bethe-Salpeter equation in the on-shell factorization scheme [47, 48], the potentials of the HMChPT lead to the dynamical generation of quasi-bound states from the s-wave scattering of heavy flavored mesons off Goldstone bosons. In particular, it provides a description of the lightest scalar and axial vector open-charm states (i.e., D 0 ∗ ( 2300 ) , D s 0 ∗ ( 2317 ) , D ₁(2430) and D _s1(2460)) as hadronic molecules, as well as predictions for their counterparts in the open-bottom sector.

The Lagrangian of HMChPT expanded at next-to-leading order (NLO) in the chiral expansion and at leading order (LO) in the inverse of the mass of the heavy meson m _H reads [21, 49–52], L = L LO + L NLO , (1) with the subscripts LO and NLO referring to the chiral power counting, and L LO = L LO ChPT + 〈 ∇ μ H ∇ μ H † 〉 − m H 2 〈 H H † 〉 − 〈 ∇ μ H ∗ ν ∇ μ H ν ∗ † 〉 + m H 2 〈 H ∗ ν H ν ∗ † 〉 + i g 〈 H ∗ μ u μ H † − H u μ H μ ∗ † 〉 + g 2 m D 〈 V μ ∗ u α ∇ β H ν ∗ † − ∇ β V μ ∗ u α H ν ∗ † 〉 ϵ μ ν α β , (2a) L NLO = L NLO ChPT − h 0 〈 H H † 〉 〈 χ + 〉 + h 1 〈 H χ + H † 〉 + h 2 〈 H H † 〉 〈 u μ u μ 〉 + h 3 〈 H u μ u μ H † 〉 + h 4 〈 ∇ μ H ∇ ν H † 〉 〈 u μ u ν 〉 + h 5 〈 ∇ μ H u μ , u ν ∇ ν H † 〉 + h ̃ 0 〈 H ∗ μ H μ ∗ † 〉 〈 χ + 〉 − h ̃ 1 〈 H ∗ μ χ + H μ ∗ † 〉 − h ̃ 2 〈 H ∗ μ H μ ∗ † 〉 〈 u ν u ν 〉 − h ̃ 3 〈 H ∗ μ u ν u ν H μ ∗ † 〉 − h ̃ 4 〈 ∇ μ H ∗ α ∇ ν H α ∗ † 〉 〈 u μ u ν 〉 − h ̃ 5 〈 ∇ μ H ∗ α u μ , u ν ∇ ν H α ∗ † 〉 , (2b) where L LO ChPT and L NLO ChPT encode the chiral Lagrangians of the pure light-meson sector. In the charm sector, H and H μ ∗ denote the antitriplets of pseudoscalar D-mesons, D 0 D + D s + , and vector D ^∗-mesons, D μ ∗ 0 D μ ∗ + D s , μ ∗ + , respectively, while in the bottom sector they correspond to the pseudoscalar B ̄ -mesons, B ̄ − B ̄ 0 B ̄ s 0 , and vector B ̄ ∗ -mesons, B ̄ μ ∗ − B ̄ μ ∗ 0 B ̄ s , μ ∗ 0 . The octet of Goldstone bosons are contained in the unitary matrix u = exp ( i Φ / 2 f π ) in the building blocks u _μ = i(u ^† ∂ _μ u − u∂ _μ u ^†) and χ ₊ = u ^† χu ^† + uχu, with the quark mass matrix χ = diag ( m π 2 , m π 2 , 2 m K 2 − m π 2 ) . For our calculations, we rely on the values from the Fit-2B in Ref. [53]. We employ the relation { h i } = { h ̃ i } that is applicable at LO in the heavy-quark mass expansion. For the specific values, please refer to our previous works [37, 38]. In the bottom sector we take advantage of the heavy-quark mass scaling of the low energy constants (LECs), { h i B } M ̂ B − 1 = { h i D } M ̂ D − 1 , for h i H ∈ { h 0 H , h 2 H , h 3 H , h 4 H M ̂ H 2 , h 5 H M ̂ H 2 } .

The tree-level potential for the process H ^(∗)Φ → H ^(∗)Φ reads V i j s , t = 1 f π 2 C LO i j 4 s − u − 4 C 0 i j h 0 + 2 C 1 i j h 1 − 2 C 24 i j 2 h 2 p 2 ⋅ p 4 + h 4 p 1 ⋅ p 2 p 3 ⋅ p 4 + p 1 ⋅ p 4 p 2 ⋅ p 3 + 2 C 35 i j h 3 p 2 ⋅ p 4 + h 5 p 1 ⋅ p 2 p 3 ⋅ p 4 + p 1 ⋅ p 4 p 2 ⋅ p 3 , (3) where s = ( p 1 + p 2 ) 2 , t = ( p 1 − p 3 ) 2 , and u = ( p 1 − p 4 ) 2 are the standard Mandelstam variables, and the superindices i, j denote the incoming and outgoing channels from the coupled-channel basis. For instance, in the sector with strangeness S = 0 and isospin I = 1/2, which are the quantum numbers of the D 0 ∗ ( 2300 ) , we have { D π , D η , D s K ̄ } , and for S = 1 and I = 0, as for the D s 0 ∗ ( 2317 ) , we have {DK, D _s η}. We refer the reader to Ref. [38] for the values of the coefficients C k i j in the isospin basis.

The partial-wave projection with angular momentum ℓ is then obtained through the relation V ℓ i j s = 1 2 ∫ − 1 + 1 d cos ⁡ θ P ℓ cos ⁡ θ V i j s , t s , cos ⁡ θ , (4) where θ is the scattering angle between the initial and final particles in the center of mass, and P _ℓ (cos θ) the Legendre polynomial of order ℓ.

A hadron gas forms once the temperature is turned on. The relative abundance of each hadron species in thermal equilibrium is determined by the corresponding thermal distribution functions, i.e., the Bose-Einstein distribution for mesons and the Fermi-Dirac distribution for baryons. At temperatures T ∼ 100 − 150 MeV light mesons become the primary components of the medium and heavy mesons behave as Brownian particles scattering of the light mesons. We exploit this scenario and employ the HMChPT potential presented above to describe the dynamics of the heavy mesons with the light mesons in the bath.

In order to incorporate the effects of the hot medium, it is necessary to follow the techniques of thermal field theory. There exist two complementary formulations of thermal field theory that can be used to describe a system in thermal equilibrium, the “imaginary-time” and the “real-time” formalisms. In the imaginary-time formalism (ITF), also called Matsubara formalism owing to the pioneering work by Matsubara [54], time is treated as a purely imaginary quantity and then one performs an analytical continuation from Euclidean to Minkowski spacetime at the end of the calculation. In the real-time formalism, in contrast, the calculation is done in Minkowski spacetime, considering explicitly the evolution in real time [55, 56]. While the latter is capable of describing systems even outside thermal equilibrium, and thus appropriate to address out-of-equilibrium properties and transport coefficients, as we will discuss in the next section, the ITF has the advantage of resembling in a more intuitive way the zero temperature field theory. For instance, the main difference in the form of the propagators and the diagrammatic structure of the perturbative expansion, in the ITF, is the acquisition of thermal weights in the phase-space integrals compared with those at T = 0, as we will show below. An approach based on the ITF has been developed in the recent years to study the properties of open heavy-flavor mesons in hot hadronic matter at vanishing baryonic density [37, 38, 57].

To compute the thermal corrections to a given quantity, such as the two-meson propagator or the self-energy, the ITF provides some simple rules that basically consist in replacing the zeroth component of the four-momenta of the particles by discrete Matsubara frequencies iω _n , with ω _n = 2nπ/β for bosons and β = 1/T, and transforming the integration over internal energies into a summation over Matsubara frequencies. Then, by using some established computational techniques based on contour integrals and analytic continuation, the calculations can be done similarly as in the vacuum field theory. For details, we encourage the reader to consult the classical Refs. [58–62].

Using the rules above, the thermal two-meson propagator takes the form G E , p ; T = ∫ d 3 q 2 π 3 ∫ d ω ∫ d ω ′ S H ω , q ; T S Φ ω ′ , p − q ; T E − ω − ω ′ + i ε 1 + f ω , T + f ω ′ , T , (5) with p ^μ = (E, p ) the momentum in the center of mass of the two-meson system. In the above equation, in addition to the medium corrections arising from the ITF, i.e., the additional weighting factors containing appropriate combinations of Bose-Einstein distribution functions f ( ω , T ) = ( e β ω − 1 ) − 1 , the meson masses are dressed by the spectral functions S H and S Φ . Note that this is a compact expression with the integrals over energy extending from −∞ to +∞.

In the case of zero temperature, it is customary to regularize the vacuum contribution to the two-meson propagator, for example, by introducing a hard cutoff in the three-momentum integration. We use Λ = 800 MeV, a value that corresponds to the scale of the degrees of freedom that were integrated out when constructing the meson-meson interaction amplitude from the effective Lagrangian, and that is consistent with the regularization scheme used in Ref. [53], from where we adopted the values for the LECs of the NLO potential. This regularization scheme is straightforward to extend to finite temperature. The tree-level potential in Eq. 3 does not change at finite temperature because in the ITF thermal corrections enter in loop diagrams [60, 61].

The thermal effects on the unitarized scattering amplitude T i j ( s ) from an incoming channel i to an outgoing channel j are then obtained by solving the on-shell Bethe-Salpeter equation in a full coupled-channel basis, with the s-wave interaction kernel of Eq. 3 and the thermal two-meson propagator of Eq. 5 (see Figure 1A): T i j E , p ; T = V i j s + V i k s G k E , p ; T T k j E , p ; T . (6)

The spectral functions dressing the meson propagators in Eq. 5 take into account the modifications due to the presence of interactions with the medium. At finite temperature the heavy meson retarded propagator is defined by D H ω , q ; T = 1 ω 2 − q 2 − m H 2 − Π H ω , q ; T , (7) where m _H is the mass of the heavy meson in the vacuum, renormalized by the vacuum contribution of the retarded self-energy Π _H (see Figure 1B). For the purpose of our calculations, using the vacuum propagator for the light meson and thus a δ-type spectral function is a good approximation, as we discussed in our previous works [37, 38].

The light-meson contribution to the self-energy of the heavy meson can be obtained by closing the light-meson line in the corresponding T -matrix element (see Figure 1C), i.e., by integrating over the light-meson four-momenta q′ ^μ = (E′, q ′). In the ITF it is defined as Π H i ω n , q ; T = − 1 β ∫ d 3 q ′ 2 π 3 ∑ m D Φ i ω m − i ω n , q ′ T H Φ i ω m , p . (8) Using the spectral Lehmann representation for the light meson propagator and the T matrix, and performing the summation over the Matsubara frequencies ω _m of the internal HΦ, it reads Π H i ω n , q ; T = 1 π ∫ d 3 q ′ 2 π 3 ∫ d E ∫ d ω S Φ ω , p − q f E , T − f ω , T E − i ω n − ω I m T H Φ E , p ; T . (9)

We note that the self-energy entering in Eq. 7 can only contain thermal corrections after mass renormalization. However, the self-energy computed with Eq. 9 contains both vacuum and thermal corrections. We regularize it by dropping the vacuum contribution, which is identified with the expression obtained when taking the limit T → 0 of Eq. 9. See [38] for details.

Finally, the spectral function necessary to dress the heavy meson in the two-meson propagator is computed from the imaginary part of the retarded meson propagator, S H ω , q ; T = − 1 π I m D H ω , q ; T . (10)

Eqs 5–10 are interrelated to each other. As a result, solving this set of coupled equations requires an iterative approach until self-consistency is achieved. This process is outlined in Figures 1A–C, where the T -matrix amplitude is depicted as a hatched blob, the perturbative amplitude V ( s ) is represented by a small dot, and the propagator of the heavy meson dressed by the medium is shown with double red lines.

Since heavy mesons have large masses compared to the surrounding light mesons their time evolution is typically described by a Fokker-Planck (or Langevin) or a Boltzmann approach. The essential components of these approaches are transport coefficients, which can be calculated from the scattering amplitudes of heavy-flavor mesons with light mesons in the hadronic gas. These transport coefficients are typically derived assuming that the light scattering partners are in thermal equilibrium, and they are often calculated as functions of temperature and momentum.

In Ref. [39] we calculated the transport coefficients of D mesons in the hadronic phase incorporating medium corrections to the scattering amplitudes. To do so, we extended the kinetic theory of D mesons using the more general Kadanoff-Baym equations, so as to account for thermal and off-shell effects. This off-shell kinetic theory also applies to describe the propagation of B ̄ mesons, as it is valid for any heavy species that can be treated as Brownian particles propagating in a medium of light mesons. In fact, the separation of scales between the heavy-meson mass and the other scales in the system that is exploited to convert the off-shell kinetic equation into a Fokker-Planck equation is larger for the B ̄ meson than for the D meson. In addition, the quasiparticle approximation that we showed to be sound for the D mesons is even better for B ̄ , since their thermal width is of the same order as that of the D mesons, but their mass is considerably larger, as we will see in Section 3.

Let us summarize our main results. For a detailed derivation, we recommend to consult our previous work [39] and references therein. Starting from the Kadanoff-Baym equations and performing a Wigner transform along with a gradient expansion [55], we arrive to the following form of the off-shell transport equation for the time ordered Green’s function of the heavy meson G H < ( X , k ) , k μ − 1 2 ∂ Re  Π R X , k ∂ k μ ∂ ∂ X μ i G H < X , k = 1 2 i Π < X , k i G H > X , k − 1 2 i Π > X , k i G H < X , k . (11) The lesser and greater Green’s functions and the self-energies in Eq. 11 are functions of the center-of-mass coordinate X = (t, X ) and the four-momentum k = (k ⁰, k ) of the external heavy meson. Note that k ⁰ and k are independent variables, although related through the non-equilibrium spectral function S H ( X , k ) . Hence the reason we denote this kinetic equation to be “off-shell”, as the heavy meson is not on its mass shell. The self-energy Π H R ( X , k ) is the extension of the retarded self-energy of Eq. 8 to the non-equilibrium case, and the lesser Π^<(X, k) and greater Π^>(X, k) self-energies can be written in terms of the (retarded) T matrix of Eq. 6 in the so-called T -matrix approximation [55, 63, 64]. Inserting appropriate definitions of these quantities, Eq. 11 can be written in the following form: k μ − 1 2 ∂ R e Π R X , k ∂ k μ ∂ ∂ X μ i G H < X , k = 1 2 ∫ d k 1 0 2 π d 3 q 2 π 3 W k 0 , k + q , k 1 0 , q i G H < X , k 0 , k + q − W k 0 , k , k 1 0 , q i G H < X , k 0 , k . (12) The off-shell collision rate of a heavy meson with energy k ⁰ and momentum k to a final heavy meson with energy k 1 0 and momentum k − q , with momentum loss q ≡ k − k ₁, is defined as W k 0 , k , k 1 0 , q ≡ ∫ d 4 k 3 2 π 4 d 4 k 2 2 π 4 2 π 4 δ k 1 0 + k 2 0 − k 3 0 − k 0 δ 3 k 2 − k 3 − q × T k 1 0 + k 2 0 + i ε , k − q + k 2 2 i G Φ > X , k 2 i G Φ < X , k 3 i G H > X , k 1 0 , k − q . (13) The labels of the momenta correspond to the choice of a generic scattering process H(k) + Φ(k ₃) → H(k ₁) + Φ(k ₂). Each of the two terms of the right-hand side of Eq. 12 can be identified as the collision gain and loss terms respectively. It is important to note that in Eq. 13 there is an implicit sum over the different species Φ and H that can interact with the external off-shell heavy meson.

Next, one may exploit the separation of scales between the meson masses to arrive to an off-shell Fokker-Planck equation for i G H < ( t , k ) . While the derivation of a transport equation for heavy mesons required the use of real-time techniques, the actual calculation of the heavy-flavor transport coefficients can be addressed in a near-equilibrium regime, in which the temperature is at least locally well defined. While thermal local equilibrium can be safely considered for the light mesons, it is also reasonable to assume that the heavy mesons are not far from equilibrium. Then Eq. 12 can be written as follows: ∂ ∂ t i G H < t , k = ∂ ∂ k i A ̂ k ; T k i i G H < t , k + ∂ ∂ k j B ̂ 0 k ; T Δ i j + B ̂ 1 k ; T k i k j k 2 i G H < t , k , (14) with Δ ^ij = δ ^ij − k ⁱ k ^j / k ², and the transport coefficients A ̂ ( k 0 , k ; T ) , B ̂ 0 ( k 0 , k ; T ) , and B ̂ 1 ( k 0 , k ; T ) defined off shell and at temperature T. The drag force coefficient is given by A ̂ k 0 , k ; T ≡ ⟨ 1 − k ⋅ k 1 k 2 ⟩ , (15) and the transverse and longitudinal momentum diffusion coefficients read. B ̂ 0 k 0 , k ; T ≡ 1 4 〈 k 1 2 − k ⋅ k 1 2 k 2 〉 , (16a) B ̂ 1 k 0 , k ; T ≡ 1 2 〈 k ⋅ k − k 1 2 k 2 〉 . (16b)

The angle brackets ⟨ F ( k , k 1 ) ⟩ denote the average of the generic quantity F ( k , k 1 ) , which is defined as F k , k 1 = 1 2 k 0 ∑ λ , λ ′ = ± λ λ ′ ∫ − ∞ ∞ d k 1 0 ∫ ∏ i = 1 3 d 3 k i 2 π 3 1 2 E 2 2 E 3 S H k 1 0 , k 1 ; T × 2 π 4 δ 3 k + k 3 − k 1 − k 2 δ k 0 + λ ′ E 3 − λ E 2 − k 1 0 T k 0 + λ ′ E 3 , k + k 3 ; T 2 × f Φ λ ′ E 3 ; T f ̃ Φ λ E 2 ; T f ̃ H k 1 0 ; T F k , k 1 , (17) where f (E _i ; T) is the equilibrium occupation number, i.e., the Bose-Einstein distribution function, and f ̃ ( E i ; T ) = 1 + f ( E i ; T ) = − f ( − E i ; T ) is the Bose enhancement factor. Eq. 17 incorporates the equilibrium quantities presented in the previous section, i.e., the equilibrium thermal scattering amplitudes and spectral functions. In particular, it is a sum of four terms (λ, λ′ = ±) that evaluate the T matrix in different energy regions. The relevance of the contribution of each of these terms to the transport coefficients is discussed in the next section.

We start with the discussion of the self-consistent results of the self-energy of the ground-state heavy mesons in a pionic medium at finite temperature, displayed as a function of the energy in Figure 2, for zero three-momentum and scaled by the mass of the heavy meson in vacuum. We show the results for three different values of the temperature of the medium, T = 80, 120 and 150 MeV, in different line styles.

The real part of the self-energy is related to the thermal correction to the mass. This is evident from the expression of the heavy-meson retarded propagator at finite temperature in Eq. 7. In the quasiparticle approximation, which we will see is well-grounded for both D and B ̄ mesons, and if the thermal propagator’s pole is close to the vacuum pole, the mass shift is roughly given by Δm _H ≈ Re Π _H (m _H , 0; T)/(2m _H ). Therefore, as shown in the first row panels of Figure 2 for charmed mesons and in the third row for bottomed mesons, the negative character of the real part of the self-energy indicates that the masses of the heavy mesons will decrease as temperature rises. The fact that the real part of the self-energy is more negative for the non-strange mesons than for the strange mesons is explained by the large attractive interaction in the D ^(∗) π and B ̄ ( ∗ ) π channels after unitarization. Furthermore, one can see that the values of the real part of self-energy over the heavy meson mass at a particular temperature are similar for D _(s) and D ( s ) ∗ , as well as for B ̄ ( s ) and B ̄ ( s ) ∗ . Although it appears to be less negative for bottom than for charm, it is of comparable size in both sectors. These findings are closely connected to the HQSFS intrinsic of the interaction. It is also important to note that the quantitative comparison of the results in the two flavor sectors may be impacted by the details of the numerical calculations, e.g., by the choice of the infinite integration limits, or the limitations of the effective theory at high energies. The authors of Ref. [57] neglected the shift of the in-medium mass of the heavy mesons by setting to zero the real part of the respective self-energies. Although small compared to the vacuum mass, |Δm _H |/m _H ∼ 1 − 2%, we consider that it is important and keep the full self-energy for the calculation of the spectral function.

The imaginary part of the self-energy relates to the thermal width acquired by the heavy meson due to interactions with pions within the medium. The panels in the second and fourth rows of Figure 2 show the imaginary part of the self-energy of the charmed and bottomed mesons, respectively, over their respective mass in vacuum. The insets provide a zoom in the region E ≈ [m _H − 2m _π , m _H + 2m _π ]. Similar features are observed in all the panels. The imaginary part of the self-energy is essentially zero for energies below m _H − 2m _π , above which it starts decreasing mildly. This initial drop at around m _H is exclusively caused by the presence of the thermal medium, which allows for the absorption of two thermal pions. These absorption processes make it possible for the scattering amplitude to be non-zero even below the two-meson threshold due to the so-called Landau cut [58, 67] of the two-meson propagator (see also our discussion in Refs. [38–40]). Our calculations show that also the self-energy of the heavy meson can reveal the effects of the Landau cut, thanks to the self-consistency of our approach. This effect becomes more relevant at higher temperatures, i.e., T ≳ 100 MeV, where the pion density is larger. A substantially larger drop takes place at m _H + 2m _π , which is the energy where the heavy meson at rest can emit two pions. This later growth of the magnitude of the imaginary part of the self-energy takes place at similar rates for all temperatures, since the emission of two pions is also possible in vacuum for a large enough energy of an off-shell heavy meson, and it is related to the unitary cut of the propagator. As a result of the combination of the Landau and unitary cut effects, and by virtue of the relation between Im Π _H and the thermal decay width, we expect the ground-state spectral functions to broaden as the temperature of the thermal medium increases. Similarly as it happened for the real parts, the magnitude of Im Π _H /m _H at a given temperature is similar in size when comparing results for pseudoscalar and vector mesons, and it is somewhat larger in the bottom sector than in the charm sector. As a reminder, the real and imaginary parts of the self-energy are connected by analyticity constraints, which means that the oscillations seen in the real part are fully determined by the structure of the imaginary part that we have described in detail here.

The spectral function for the ground-state heavy mesons follows the standard definition in terms of the retarded propagator—see Eqs 7, 10—, in which the self-energy is responsible for the thermal corrections with respect to the vacuum propagator. In Figure 3 we show the energy dependence of the spectral function of the charmed mesons (top panels) and the bottomed mesons (bottom panels) at rest, at the same temperatures as for the self-energy described above. The vertical solid lines depict the corresponding value of the mass in vacuum. From these plots, the drop of the mass and the increase of the width anticipated from analyzing the self-energies become manifest. This is evident as the maximum of the spectral function shifts towards lower energies and it becomes wider with increasing temperature.

In the quasiparticle approximation, the spectral function admits a Lorentzian shape peaked at the quasiparticle energy E _k (T) (with M(T) ≡ E _k at rest) and a spectral width γ _k (T) ≪ E _k (T). For the spectral functions in Figure 3, which are narrow, the quasiparticle approximation is indeed justified. Figure 4 shows the values of the mass (left panels) and the decay width (right panels) as a function of the temperature determined by analyzing the position and the width of the peak of the spectral functions. For the charmed mesons, we find a reduction in mass ≈ 45  MeV and ≈ 25  MeV for the non-strange and the strange states, respectively, at the highest temperature T = 150 MeV, and corresponding thermal widths of ≈ 70  MeV and ≈ 20  MeV. For the respective bottomed mesons, the reduction in mass is ≈ 30  MeV and ≈ 20  MeV, and the acquired width is ≈ 90  MeV and ≈ 30  MeV. The thermal masses of the D ^(∗) and D s ( ∗ ) mesons were calculated using lattice QCD simulations in Ref. [68]. We show their results in the top left panel of Figure 4, although for the comparison with our calculations one has to keep in mind that a systematic shift is to be expected due to the use of heavier than physical pions in the lattice. Indeed, smaller thermal modifications at a given temperature are consistent with the lower abundances of heavier pions.

We remind that in the self-consistent calculations of the self-energies and the spectral functions, only the impact of the thermal pions is taken into account, arguing the contribution of other light mesons to be presumably suppressed. The contribution of each of the light mesons in the medium to the thermal width of the heavy meson can be easy analyzed from its definition in terms of the retarded self-energy in the quasiparticle approximation, Γ k = − z k E k I m Π E k , k ; T , (18) with z _k ≈ 1. Since Im Π(E _k , k ; T) is given by the integration over the imaginary part of the unitarized scattering amplitude, I m T , and we have access to all the matrix elements T i j , we can readily assess the effect of the four elastic channels for the interactions of a heavy meson with the light pseudoscalars ( π , K , K ̄ , η ) . We showed in Ref. [39] that the effect of the inelastic channels is negligible. In Figure 5 we show, in logarithmic scale, the contribution to the width of the D meson (left) and the B ̄ meson (right) coming from the different light mesons, averaged over momenta [39]. At low temperatures, the kaons and the η mesons have a negligible contribution because of their small abundances, as expected. The only relevant contribution is that of the pions. Close to T = 150 MeV, the more massive mesons contribute several MeV to the heavy-meson decay width, but are still subdominant compared to the pion.

The process of unitarization of the scattering amplitude of Eq. 6 leads to the emergence of two poles in the sectors with strangeness S = 0 and isospin I = 1/2 that correspond to the two-pole structure of the D 0 ∗ ( 2300 ) , in the case of total angular momentum J = 0, and of the D ₁(2430), in the case of J = 1. The same applies for the counterparts in the bottom sector. In the sectors with (S, I) = (1, 0), the poles of the D s 0 ∗ ( 2317 ) and the D _s1(2460) emerge for J = 0 and J = 1, respectively, as well of their bottomed analogues. The characterization of these states requires the analytical continuation of T to complex energies. The pole position in the complex-energy plane provides the mass and the half-width of these states. While this is a well-established procedure at T = 0, the poles search at finite temperature is a complex task for two reasons. Firstly, one has to deal with the analytic continuation of imaginary frequencies, and secondly, a numerical search of a singularity in the complex plane within self-consistency is computationally challenging. Alternatively, to determine the spectral properties of the dynamically generated states at finite temperature, we use the imaginary part of T on the real axis, shown in Figure 6, as a proxy for their spectral shape. From the several coupled channels in each sector, we choose to plot the diagonal element T i i for the channel i to which the state couples more strongly in vacuum. The numerical values of the vacuum properties and effective couplings are given in our previous works [38, 40].

In the cases with zero strangeness, the proximity of the position of the resonances to channel thresholds gives rise to complicated structures. However, one can still clearly recognize a gradual evolution of the peaks and widths as T increases. The strange sectors are more straightforward. The T = 0 delta-type spectral function (i.e., bound state) acquires a non-zero width, and the shift and widening of the peak is comparable to that of the ground states. Nevertheless, an increase in strength is visible on the right-hand side of the distributions. This asymmetry can be explained by the fact that the channel threshold is not sharp anymore due to the widening of the D ^(∗) or B ̄ ( ∗ ) meson, and it is lowered in energy due to the decrease of the heavy-meson mass with temperature. Both of these effects open the phase space for decay into this channel at lower energies.

Now, we discuss the results of the transport coefficients for a heavy meson propagating through the hadronic medium. We will specifically focus on the comparison between D and B ̄ mesons. We start with two of the transport coefficients defined in momentum space, the drag force in Eq. 15 and the transverse diffusion coefficient of Eq. 16a. To reduce the number of variables we will present the results in the so-called static limit k → 0, i.e., for a low-momentum heavy meson. In the off-shell version of the Fokker-Planck equation (Eq. 14), these transport coefficients also depend on the value of k ⁰. Since the quasiparticle approximation is excellent for the temperatures considered, we will set k ⁰ to the quasiparticle mass, that is k ⁰ = E _k→0 = M(T). Then, the only remaining dependence is on the temperature. In the left panel of Figure 7 we present the comparison of the drag force coefficient for the two heavy flavors as functions of temperature, while the right panel displays the transverse diffusion coefficients. In the static limit, we have checked that the longitudinal diffusion coefficient B ₁(k ⁰, k → 0; T) is degenerate with B ₀(k ⁰, k → 0; T). As explained in Ref. [39], these coefficients can be computed with different degrees of approximation, but in this work we only present the complete off-shell computation. This calculation incorporates: 1) the exact thermal spectral function of the heavy meson as required in the average of Eq. 17; 2) the full thermal T -matrix appearing in the same equation; 3) quantum effects encoded in the Bose enhancement factors; and 4) all kinematic processes allowed by energy-momentum conservation, including number-conserving (2 ↔ 2) and number-violating processes (1 ↔ 3).

The 2 ↔ 2 scatterings are described by the λ = λ′ terms in Eq. 17, while the 1 ↔ 3 processes take λ = −λ′ This fact can easily be grasped by looking at the signs of the energy conservation delta. In Ref. [39] we reported that the number-violating processes contribute very little to the transport coefficients, while the 2 ↔ 2 collisions make the leading contributions. Among the latter, the case λ = λ′ = + corresponds to the standard term in which the binary collision is taking place at an energy corresponding to the sum of the incoming energies k ⁰ + E ₃ (cf. Eq. 17). However the case λ = λ′ = −, corresponds to a binary collision in which the scattering matrix is evaluated at the energy difference k ⁰ − E ₃. For typical energies around the quasiparticle masses, this difference probes the kinematic region below the two-particle threshold. For interactions computed in vacuum, this region has a vanishing T -matrix amplitude, and this entire process can be safely neglected. However, for interactions self-consistently calculated at T ≠ 0 the T -matrix has a non-vanishing support in this region, due to the Landau cuts (see Section 3.1). We have proven in Ref. [39] that the contribution of these processes cannot be neglected and it becomes comparable to the contribution stemming from the unitary cut. Eventually, we have obtained transport coefficients that are up to a factor of three larger with respect to previous results at the highest temperatures. We nonetheless agree at low temperatures, where the Landau cut disappears.

To compare the results between D and B ̄ mesons in Figure 7 we recall that a simple nonrelativistic approximation for the momentum-space diffusion coefficient B ₀ has no leading-order dependence on the heavy mass [69], but a dependence proportional to the total cross section, B ₀ ∝ σ. From Figure 6 we have learnt that the cross sections (proportional to Im  T ) of B ̄ mesons are a factor 2-3 larger than the those of D mesons. This explains why the B ₀ coefficient is 2-3 times larger for bottom than for charm. On the other hand, the drag force coefficient does have a leading dependence on the heavy mass M _H . Again from Ref. [69] the nonrelativistic expression goes like A ∝ σ/M _H . Therefore, going from charm to bottom, the gain factor from the cross section is approximately compensated by the reducing factor due to the increasing mass, as M D / M B ̄ ≃ 1 / 3 . Therefore we expect that A B ̄ m e s o n ≲ A D m e s o n , which is what we observe in the left panel of Figure 7.

Finally we plot in Figure 8 the so-called spatial diffusion coefficient, D s T = lim k → 0 T 2 B 0 E k , k ; T = lim k → 0 T A E k , k ; T M H , (19) where we express it with calligraphic font not to confuse it with D _s mesons.

In the left panel of this figure we show our result for this coefficient at low temperatures for D mesons (solid red line) and for B ̄ mesons (dotted-dashed blue line). As expected, D s is lower for B ̄ mesons than for D mesons, while the general trend is a monotonically decreasing function of temperature. As for the high temperature side, we plot two extractions from relativistic heavy-ion collisions at high energies using Bayesian analyses to estimate the temperature-behavior of this coefficient for the two flavors separately [70]. While a clear ordering cannot be settled in this extraction, a likely continuous matching can be observed around T _c ≃ 155 MeV. According to these results, the absolute minimum of the spatial diffusion coefficient might happen at the transition temperature T _c . On the right panel of Figure 8 we plot our results at low temperatures together with results of a quasiparticle model for the quark-gluon plasma at high temperature, which can distinguish charm and bottom quarks [71]. An approximate matching is also seen in T _c and the flavor-mass ordering of D s is consistent in both sides of the transition. Close to T _c a more refined model including a mixed phase with hadronization processes should be able to fill the gap.