Thermal Energy of a Charm-meson Molecule in a Pion Gas

The thermal corrections to the propagator of a loosely bound charm-meson molecule in a pion gas are calculated to next-to-leading order in the heavy-meson expansion using a zero-range effective field theory. Ultraviolet divergences in the charm-meson-pair self energy are canceled by corrections to the charm-meson-pair contact vertex. Terms that are singular at the charm-meson-pair threshold can be absorbed into thermal corrections to the rest energies and kinetic masses of the charm-meson constituents. The remaining terms reduce to a thermal correction to the binding momentum that is proportional to the pion number density and suppressed by the pion/charm-meson mass ratio. The correction gives a tiny decrease in the binding energy of the charm-meson molecule relative to the charm-meson-pair threshold in the pion gas and a change in its thermal width that is small compared to the thermal widths of the charm-meson constituents. These results are encouraging for the prospects of observing $X(3872)$ and $T_{cc}^+(3875)$ in the expanding hadron gas produced by heavy-ion collisions.


I. INTRODUCTION
Nature has provided us with at least two exquisite examples of loosely bound hadronic molecules: X(3872) (or, more concisely, X), also known as χ c1 (3872), which was discovered by the Belle collaboration in 2003 [1], and T + cc (3875) (or, more concisely, T + cc ), which was discovered by the LHCb collaboration in 2021 [2].Since their discovery, there have been numerous studies of these and other exotic heavy hadrons.For recent reviews, see Refs.[3,4].Their binding energies relative to the appropriate charm-meson-pair thresholds and their decay widths are both over an order of magnitude smaller than the energy scale associated with pion exchange between the charm mesons: m 2 π /M ≈ 10 MeV, where m π is the pion mass and M is the charm-meson mass.In addition, the charm-meson constituents have decay widths that are over two orders of magnitude smaller than m 2 π /M .The X and T + cc are much more loosely bound than the deuteron, whose binding energy is about 40% of the energy scale 5 MeV associated with pion exchange between the nucleons.We refer Ref. [5] for a review of the universality of a loosely bound molecule and Ref. [6] for a general introduction to hadronic molecules.
The tiny binding energies of X and T + cc can be exploited to develop quantitative treatments of their properties, including their interactions with other hadrons.The simplest effective field theory (EFT) that can be applied to a loosely bound charm-meson molecule is a zero-range effective field theory (ZREFT) for nonrelativistic charm mesons, in which they interact only through contact interactions [7].This ZREFT is the analog for charm mesons of the pionless EFT that has been applied to the few-nucleon problem [8].An alternative EFT with a greater range of validity is XEFT, which describes nonrelativistic charm mesons and pions [9].Calculations in XEFT can be simplified by using a Galilean-invariant formulation of XEFT that exploits the fact that the D * -D mass difference is approximately equal to the pion mass [10,11].
A challenging problem to which EFT's for charm mesons may be able to provide some insight is the production of X and T + cc in relativistic heavy-ion collisions.A central relativistic heavy-ion collision is believed to produce a region of quark-gluon plasma that remains in local thermal equilibrium as it expands and cools.At a temperature of 156 MeV, there is a hadronization transition from the quark-gluon plasma to a hadron resonance gas, which continues to expand and cool in local thermal equilibrium.At a temperature around 115 MeV, the hadron gas reaches kinetic freeze-out and goes out of thermal equilibrium.The system then continues to expand by free-streaming of the hadrons.
The CMS collaboration has observed the production of X(3872) in Pb-Pb collisions at the LHC [12].The ratio of the production rates of X and ψ(2S) seems to be about an order of magnitude larger than that from hadronic production mechanisms, such as prompt production in pp collisions and exclusive production in B-meson decays.It is particularly surprising that a system with a binding energy less than 1 MeV can survive in a hadronic environment whose temperature is greater than 100 MeV.
The behavior of a loosely bound charm-meson molecule in the hadron resonance gas near the hadronization temperature is a very challenging problem.One complication is the many strongly interacting light-hadron resonances that must be taken into account.Another complication is the restoration of chiral symmetry near the hadronization temperature, which requires scalar and pseudoscalar charm mesons to become degenerate and requires axialvector and vector charm mesons to become degenerate.
A much simpler problem is the behavior of a loosely bound charm-meson molecule in the hadron gas near the kinetic freeze-out temperature.The most abundant hadrons in the hadron gas are pions.Near kinetic freeze-out, the abundance of kaons is smaller by about a factor of 5 and other hadrons are even less abundant.The hadron gas can therefore be approximated by a pion gas.The temperature of the pion gas is high enough that the pions must be treated as relativistic particles.However the temperature is low enough that it may be possible to describe pion interactions using a chiral effective field theory (χEFT).Since the kinetic freezeout temperature is orders of magnitude larger than the binding energy of a loosely bound molecule, one might expect that ZREFT or XEFT are simply not applicable at such a high temperature.We show in this paper that ZREFT can be applied to the loosely bound charm-meson molecule in the pion gas by first integrating out thermal pions in favor of temperature-dependent modifications of the parameters of ZREFT.The thermal properties of the molecule can then be calculated in terms of those T -dependent parameters.We calculate thermal corrections to the parameters of ZREFT using the heavymeson expansion, which is an expansion in powers of charm-meson kinetic energies divided by the pion mass.At leading order (LO) in the heavy-meson expansion, the only T -dependent changes in the ZREFT parameters are thermal corrections to the complex rest energies of the charm mesons.At next-to-leading order (NLO) in the heavy-meson expansion, the changes in the ZREFT parameters are much more complicated.However, we find that the only correction to the energy of the loosely bound molecule relative to the T -dependent charmmeson-pair threshold comes from a small T -dependent correction to the complex binding momentum of the molecule.
The rest of this paper is organized as follows.In Section II, we describe the amplitude whose pole corresponds to a loosely bound molecule.In Section III, we introduce notation for the properties of charm mesons and pions and we describe the pion gas that can be produced by a relativistic heavy-ion collision.In Section IV, we calculate the self energies of charm mesons in the pion gas to NLO and we determine thermal corrections to the parameters of ZREFT.In Section V, we calculate the self-energy of a charm-meson pair in the pion gas to NLO.In Section VI, we determine the thermal mass shift and the thermal width of the loosely bound charm-meson molecule.We summarize our results in Section VII.In Appendices A and B, we give results for integrals over the pion momentum distribution and for integrals over the relative momentum of a charm-meson pair.In Appendix C, we give the Feynman rules used in our calculations.

II. LOOSELY BOUND MOLECULE
In this Section, we present the simplest approximation to the amplitude whose pole corresponds to a loosely bound molecule.We also calculate the effects on the amplitude from corrections to the propagators of the constituents of the molecule.

A. Pair Propagator
We consider a loosely bound molecule X whose two constituents are a particle D with kinetic mass M and rest energy ε and a particle D * with kinetic mass M * and rest energy ε * .The nonrelativistic propagator for D with energy E (relative to M ) and momentum p is i/[E − ε − p 2 /(2M ) + iϵ].The propagator for D * has an analogous expression.The amplitude for the propagation of D and D * between contact interactions can be represented by the bubble diagram in Fig. 1.The amplitude with nonrelativistic propagators is linearly ultraviolet divergent.The ultraviolet cutoff on the loop momentum k is most conveniently implemented by subtracting and adding to the integrand a term proportional to 1/k 2 and then imposing the cutoff on the additional term.The amplitude for total energy E (relative to M + M * ) and total momentum P is where µ = M * M/(M * + M ) is the D * D reduced mass.We have imposed a momentum cutoff |k| < (π/2)Λ on the ultraviolet-divergent integral.The square-root function S 0 in Eq. ( 1) is where E cm is the center-of-mass energy, and M X = M * + M is the kinetic mass of X.
The D * and D can scatter through a contact interaction with vertex i C 0 , which is represented by the first diagram in Fig. 2. They can also scatter through bubble diagrams, as in the second and third diagrams in Fig. 2. The bubble diagrams form a geometric series.For the molecule to be loosely bound, there must be a fine tuning of C 0 to nearly cancel the term proportional to Λ in Eq. ( 1): where γ X is small compared to Λ.The sum of the geometric series of diagrams in Fig. 2 is This amplitude has a square-root branch point in E at E cm = ε * + ε.If γ X > 0, it also has a nearby pole on the real axis at E cm = E X , where is the energy of the loosely bound molecule at zero 3-momentum.Its binding momentum is γ X , and its binding energy relative to the D * D scattering threshold is γ 2 X /(2µ).The amplitude in Eq. ( 5) is the D * D scattering amplitude in ZREFT at leading order.It can also be interpreted as the propagator for a local operator that annihilates D and D * .This amplitude can be obtained diagrammatically by omitting the first diagram in Fig. 2, amputating the D and D * legs on the remaining diagrams, and then taking the limit Λ → ∞.The propagator differs from the scattering amplitude by the absence of the term i C 0 from the first diagram in Fig. 2, but that term goes to 0 in the limit Λ → ∞.We will refer to the amplitude on the right side of Eq. ( 5) as the D * D propagator.The corrections to the D * D propagator, which is a 2-point Green function, are much simpler than the corrections to the D * D scattering amplitude, which is a 4-point Green function.
The complete D * D propagator includes all the corrections to the sum of the bubble diagrams in Fig. 2. The corrections to the bubble amplitude for the propagation of D * D between successive contact interactions give an additive correction Σ(E cm , P ) to the denominator of the D * D propagator in Eq. ( 5).We will refer to Σ, which has dimensions of momentum, as the D * D self energy.Unless there is exact Galilean invariance, Σ can have additional dependence on P beyond its dependence through E cm .There are D * D-irreducible corrections to each of the contact vertices in Fig. 2. The momentum-independent corrections can be taken into account by replacing C 0 by a contact vertex C 1 , which also must be fine tuned to ensure that the pole remains near the threshold.
If X is a loosely bound molecule, the complete D * D propagator must have a square-root branch point near E cm = ε * + ε, a pole near E cm = ε * + ε when P is small, and no other nearby singularities.The branch point comes from a function S 1 (E cm , P ) that can be chosen such that the only dependence of S 2  1 on E is an additive term −2µE.The behavior of the complete D * D propagator near the branch point must then have the form where Z X and δγ X are constants independent of E cm and P .The additional corrections to the denominator in Eq. ( 7) are an expansion in powers of S 1 beginning at order S 2 1 .The factor Z X in the numerator is determined by the condition that the coefficient of S 1 in the denominator is 1.The constant δγ X in the denominator is a correction to the binding momentum of the molecule.

B. Constituent Propagator Corrections
There are corrections to the propagators of D and D * from their interactions.The corrections to the D propagator can be expressed in terms of a self energy Π(E, p).The effects of the self energy include a shift δε in the rest energy of D and a correction δZ to the residue of the pole in its propagator.If Galilean invariance is not an exact symmetry, the inverse kinetic mass 1/M can also be modified by a multiplicative factor 1 + ζ.The behavior of the complete D propagator near its pole is where the additional corrections to the denominator are an expansion in powers of p 2 beginning at order p 4 .The factor 1 + δZ in the numerator is determined by requiring the coefficient of E in the denominator to be 1.The corrections to the D * propagator can be expressed in terms of a self energy Π * (E, p).The behavior of the complete D * propagator near its pole can be expressed in a form analogous to Eq. ( 8) with constants δε * , δZ * , and ζ * .
The form of the complete D propagator in Eq. ( 8) implies a change in the kinetic energy of D by the factor 1 + ζ.If ζ and ζ * are small compared to 1, the changes in the kinetic energies are small at small momentum p, but they become increasingly large as p increases.They therefore can have a large effect on the ultraviolet behavior of Green functions, and this must be compensated by changes in the parameters of ZREFT.These parameters include the D * D contact vertex C 0 in Eq. ( 4), which must be fine tuned to a new value C 1 to compensate not only for short-distance corrections to the D and D * propagators but also for short-distance corrections to the D * D contact interaction.
In the geometric series of bubble diagrams in Fig. 2, the D propagator corrections can be taken into account by replacing the D propagator by the right side of Eq. ( 8) with the denominator truncated after the p 2 term.The D * propagator corrections can be taken into account by making the corresponding change in the D * propagator.Changes in the contact vertex can be taken into account by replacing C 0 by C 1 .The sum of the geometric series analogous to Eq. ( 5) is where the square-root function S 1 is The modified reduced mass µ 1 is obtained from µ by making the substitutions The constant ζ X is The D * D contact vertex C 1 required to compensate for constituent propagator corrections can be determined by identifying the ultraviolet-divergent terms −Λ in the denominators of Eqs. ( 5) and ( 9).The first term in the denominator of Eq. ( 9) must then be identical to the corresponding constant (2π/µ)/C 0 in Eq. ( 5).The required contact vertex can be expressed as 2π/µ We have replaced C 0 by C 0 + δC to allow for additional corrections to the contact vertex required to compensate for the effects of interactions between constituents.The D * D self energy Σ(E cm , P ) in Eq. ( 7) must include terms that change S 0 (E cm ) in Eq. ( 2) into S 1 (E cm , P ) in Eq. (10).It must therefore include terms that are singular at the branch point of S 0 (E cm ).The expansion of S 1 (E cm , P ) to first order in δε * , δε, and ζ X is At successively higher orders, there are increasingly singular terms in the form of increasing powers of 1/S 0 (E cm ).These singular terms must be resummed to all orders to change the function S 0 (E cm ) into S 1 (E cm , P ).

III. CHARM MESONS IN A PION GAS
In this Section, we introduce notation for some of the properties of charm mesons and pions.We describe the pion gas that can be produced by a heavy-ion collision, and we specify the momentum distribution for pions in the pion gas.We also define thermal averages that appear in the self energies of the charm mesons and of the charm-meson pair in the pion gas.

A. Charm Mesons and Pions
We consider a loosely bound charm-meson molecule that is a bound state of a vector charm meson and a pseudoscalar charm meson.We denote the molecule by X and the charm mesons by D * a and D b , where a and b are charm-meson flavor indices.(These charm mesons may contain a charm antiquark instead of a charm quark.)In the case where X is T + cc (3785), the constituents are D * + D 0 , which are the (a, b) = (1, 2) members of a charm-meson isospin doublet.In the case where X is X(3872), the constituents are (D * 0 D0 + D * 0 D 0 )/ √ 2, which is a superposition of the (a, b) = (2, 1) and (1, 2) members of charm-meson isospin doublets.
We treat pions as relativistic particles.We denote the mass of the pion with flavor i, which can be +, 0, or −, by m πi .The flavor-averaged pion mass is m π = 138.0MeV.We represent the propagator for π by a dashed line.The self interactions of pions in chiral effective field theory (χEFT) at leading order are determined by the pion decay constant f π = 131.7 MeV.
We treat charm mesons as nonrelativistic particles.We denote the masses of D represent the propagator for D by a solid line.Since the mass of D * is close to the sum of the masses of Dπ, we represent the propagator for D * by a double line (solid+dashed).We take the vertices for interactions between charm mesons and pions to be the nonrelativistic form of those in heavy-hadron χEFT at leading order.The interaction parameters are f π and a dimensionless coupling constant g π = 0.520, whose value is determined from the decay rate for D * + → D 0 π + .The Feynman rules used in our calculations are given in Appendix C.
We treat the loosely bound charm-meson molecule using a ZREFT for the nonrelativistic charm mesons analogous to that described in Sec.II A. The energy of X relative to the D * a D b threshold can be expressed as ε X = −γ 2 X /(2µ), where γ X is the binding momentum and µ = M M * /(M * + M ) is the reduced kinetic mass.The assumption that X is loosely bound is equivalent to The most accurate determinations of the mass and width of the X(3872) resonance have been made by the LHCb collaboration [13,14].A fit using a Breit-Wigner line shape gives E BW = −0.07± 0.12 MeV for the energy relative to the D * 0 D0 threshold and Γ BW = 1.19 ± 0.19 MeV for the width.Using a Flatté parameterization that takes into account the nearby D * 0 D0 threshold and the D * 0 width, the best fit to the X line shape in the J/ψ π + π − decay channel gives (+25 − 140 i) keV for the pole energy relative to the threshold [13].The complex binding momentum γ X can be determined by identifying this energy with −γ 2 X /(2µ) − i Γ * 0 /2, where µ is the D * 0 D0 reduced mass and Γ * 0 = 55.4 keV is the predicted D * 0 decay width.The solution with a positive real part, corresponding to a bound state, is γ X = (9.3+ 11.6 i) MeV.The positive imaginary part takes into account short-distance decay modes of X, such as J/ψ π + π − and J/ψ π + π − π 0 .These decays account for most of the width Γ X = 280 keV determined by the pole energy of X in Ref. [13].
Accurate determinations of the mass and width of the T + cc (3875) resonance have been made by the LHCb collaboration [15].A fit to the T + cc line shape in the D 0 D 0 π + decay channel using a Breit-Wigner line shape gives E BW = −273 ± 61 keV for the energy relative to the D * + D 0 threshold and Γ BW = 410 ± 165 keV for the width.Using a unitarized model that takes into account the nearby D * + D 0 threshold, the central value of the pole energy relative to the threshold is (−360 − 24 i) keV [15].The complex binding momentum γ T can be determined by identifying the pole energy with −γ 2 T /(2µ) − i Γ * + /2, where µ is the D * + D 0 reduced mass and Γ * + = 83.4keV is the measured D * + decay width.The solution with a positive real part, corresponding to a bound state, is γ T = (26.4− 0.6 i) MeV.The small negative imaginary part takes into account decays from the small D * 0 D + component of the bound state.The D * 0 D + scattering threshold is higher than the D * + D 0 threshold by only 1.4 MeV.The D * 0 D + component can decay into D + D 0 π 0 through a decay of the D * 0 constituent, and this interferes destructively with the decay D * + → D + π 0 from the dominant D * + D 0 component.Estimates of the decay width of T + cc that take into account the coupled D * + D 0 and D * 0 D + channels range from 36 to 78 keV [16][17][18][19][20][21][22][23], all of which are smaller than the width Γ * + of D * + .The width Γ T = 48 keV determined by the pole energy of T + cc in Ref. [15] is a little more than half of Γ * + .

B. Pion Gas
A central relativistic heavy-ion collision can produce a region of quark-gluon plasma that expands and cools while in thermal equilibrium.When it reaches the hadronization temperature of 156 MeV, it undergoes a transition to a hadron resonance gas, which continues to expand and cool while in thermal equilibrium until kinetic freeze-out.For Pb-Pb collisions at the LHC, the kinetic freeze-out temperature T kf is estimated to be 115 MeV [24].At temperatures T near T kf , the hadron resonance gas can be approximated by a pion gas.The momentum distribution of the pions is a Bose-Einstein distribution 1/(e βωq − 1), where ω q = m 2 π + q 2 and β = 1/T .If the pions are in thermal equilibrium at temperature T , the number density for each of the three pions is At T kf = 115 MeV, the equilibrium pion number density is 1/(3.9fm) 3 .
After kinetic freeze-out, the hadron gas is no longer in thermal equilibrium, but it continues to expand by the free-streaming of hadrons.It can be approximated by a pion gas with a decreasing number density n π and a fixed temperature T kf .(Such a system could also be produced by the isothermal expansion of a pion gas that was in thermal equilibrium at temperature T kf .)The pion momentum distribution in both the expanding and cooling pion gas before kinetic freeze-out and the expanding pion gas after kinetic freeze-out can be described by Before kinetic freeze-out, T decreases with time and n π = n (eq) π is determined by T using Eq. ( 15).After kinetic freeze-out, n π decreases with time and the temperature remains fixed at T = T kf .A dimensionless number that characterizes the size of the effects of thermal pions is n π /(f 2 π m π ).At T kf = 115 MeV, this number is small: n (eq) π /(f 2 π m π ) = 0.052.There have been many previous discussions of the production of exotic heavy hadrons in heavy ion collisions [24][25][26][27][28][29][30], including some that focused on X(3872) and T + cc (3875) [31][32][33][34][35][36][37][38].Most treatments of hadronic molecules have taken into account their size, but are otherwise uninformed about the physics of loosely bound molecules.There have been previous calculations of the thermal mass shift and the thermal width of X(3872) in a hadron gas [39,40].The primary goal of this paper is the calculation of these properties for a loosely bound charm-meson molecule in a pion gas.
We proceed to summarize the important energy and momentum scales of the system.There is a large energy scale set by the temperature T , which is comparable to m π .The typical momentum q and energy ω q of a pion are order m π .There is a much larger energy scale set by the charm-meson masses M * and M .Isospin splittings and the D * -to-Dπ mass differences ∆ cd − m πi are all smaller than m 2 π /M ∼ 10 MeV.The strong inequality m 2 π /M ≪ m π allows a charm-meson propagator with momentum of order m π to be expanded in powers of its kinetic energy, which is order m 2 π /M , and in powers of isospin splittings, which we take to also be order m 2 π /M .We refer to this expansion as the heavy-meson expansion.In an amplitude for the propagation of a charm-meson pair, the heavy-meson expansion produces integrals over the relative momentum of a charm-meson pair that are ultraviolet divergent.These integrals can be defined by imposing the ultraviolet momentum cutoff Λ of ZREFT, which we take to be order m π .The binding momentum γ X of the charm-meson molecule is assumed to be much smaller than m π , which implies that the binding energy γ 2 X /(2µ) is much much smaller than m 2 π /M .The hierarchy of energy scales can be succinctly summarized as

C. Thermal Averages
We use angular brackets to denote the average over the Bose-Einstein momentum distribution of a pion.The thermal average of a function F (q) of the pion momentum is The thermal average depends on the temperature T but not on the pion number density n π .It may be sensitive to the flavor i of the pion, in which case the pion energy in Eq. ( 17) should be replaced by ω iq = m 2 πi + q 2 .The propagators of charm mesons in the pion gas involve thermal averages over the pion momentum distribution.Some of the thermal averages have the simple form ⟨(q 2 ) n /ω m q ⟩, but other thermal averages are more sensitive to isospin splittings.The thermal average that appears in the charm-meson self energies at leading order (LO) in the heavy-meson expansion is where ∆ cd is the and m πcd is an alternative notation for the mass m πi of the pion produced by the transition D * c → D d π i .The thermal averages that appear in the charm-meson self energies at next-to-leading order (NLO) in the heavy-meson expansion include where n is 1 or 2, and The thermal averages in Eqs. ( 18) and ( 19) can be defined by the limit ϵ → 0 + .In Appendix A 2, their real parts are expressed as principle-value integrals over a momentum variable and their imaginary parts are evaluated analytically.The thermal averages can be expanded in powers of isospin splittings divided by m π .The leading terms in the expansions of their real and imaginary parts are given in Appendix A 3. The thermal average H cd in Eq. ( 20) has terms that diverge as 1/ϵ as ϵ → 0. The divergence can be regularized by taking into account the widths of the charm mesons as described in Appendix A 4.

IV. CHARM-MESON SELF ENERGIES
In this Section, we calculate the self energies of pseudoscalar and vector charm mesons in the pion gas to NLO in the heavy-meson expansion.The self energies are used to identify thermal corrections to parameters of ZREFT.

A. Corrections from Pion Forward Scattering
In a system of pions and nonrelativistic charm mesons with temperature T , the charm meson has a typical kinetic energy of order T and a typical momentum of order √ M T .If T is order m π , the typical charm-meson kinetic energy is order m π .However, for simplicity we will calculate the self energies of charm mesons with a parametrically smaller kinetic energy E of order m 2 π /M .The self energies are calculated at leading order in the pion interactions and to NLO in the heavy-meson expansion.
The effect of the pion gas on the propagation of charm mesons can be taken into account using thermal field theory.In the imaginary-time formalism, the self energy of a D meson at leading order is the sum of the two one-loop diagrams in Fig. 3.The second diagram, which has an internal pion line and an internal charm-meson line, includes a zero-temperature part, a thermal part from the coherent forward scattering of an on-shell pion, and a thermal part from the coherent forward scattering of an on-shell charm meson.We ignore the zerotemperature part, because it is taken into account in the parameters of the ZREFT.We ignore the thermal charm-meson part, because the thermal distribution of the charm meson has an exponential suppression factor exp(−M/T ).The thermal pion part can alternatively be represented by the forward-scattering diagram in Fig. 4.An on-shell pion emerges from one open circle with momentum q and flavor i and is scattered back into the other open circle with the same momentum and flavor.The amplitudes must be added coherently by multiplying them by f π (ω q )/(2ω q ), where f π (ω q ) is the pion momentum distribution, integrating over q with measure d 3 q/(2π) 3 , and summing over the three flavors i.

B. Pseudoscalar charm meson
The complete propagator for a pseudoscalar charm meson D b with energy E (relative to its kinetic mass M ) and momentum p has the form on the left side of Eq. ( 8), with rest energy ε b and self energy Π b (E, p).The expansion of the denominator near the pole in E has the form   E, p) is the sum of the diagrams with D legs and pion legs amputated, weighted by f π (ω q )/(2ω q ), integrated over q, and summed over the pion flavor i.
in Fig. 4 is where the sum is over the flavor of the intermediate and m πcb is the mass of the pion produced by the transition D * c → D b π.The angular brackets, which are defined in Eq. ( 17), indicate the average over the Bose-Einstein distribution for the momentum q of the pion with flavor cb.
The terms ω cbq and ∆ cb in the denominator of the self energy in Eq. ( 22) are order m π .The heavy-meson expansion is obtained by expanding Eq. ( 22) in powers of the small energies E − ε b and (p ± q) 2 /(2M * ), which we consider to be order m 2 π /M .Terms with odd powers of q vanish after averaging over the directions of the pion momentum.The LO term is a constant independent of E and p: where F cb is the thermal average defined in Eq. (18).The NLO term is a linear function of E and p 2 : where G n,cb is the thermal average defined in Eq. (19).We can read off the corrections δε b , δZ b , and ζ b defined by the expansion in Eq. ( 21): The D self-energy in Eq. ( 22) comes from P-wave pion interactions through the second diagram of Fig. 4. The contribution from S-wave pion interactions through the first diagram in Fig. 4 is 0. This zero implies a cancellation between the contributions from πD scattering in the channels with total isospin 1  2 and 3 2 .The S-wave contribution to the D self energy at leading order in the chiral expansion can be expressed as where the two canceling terms in the last factor come from total isospin 1 2 and 3 2 .At higher orders in the chiral expansion, the canceling factor (1−1) is replaced by terms suppressed by powers of T /(4πf π ) or m π /(4πf π ).For a pion gas at T kf = 115 MeV, the canceling S-wave contributions to the D thermal mass shift at leading order in Eq. ( 26) are ±7.2MeV.These canceling S-wave contributions are larger than the P-wave contribution in Eq. (25a).This suggests that S-wave contributions suppressed by powers of m π /(4πf π ) or T /(4πf π ) could be significant.
The exact contribution to the D self-energy from the coherent forward scattering of lowenergy pions in the S-wave channel can be expressed as where µ = M m π /(M + m π ) is the πD reduced mass and a (1/2) πD and a (3/2) πD are the S-wave scattering lengths in the channels with total isospin 1 2 and 3 2 .At leading order in the chiral expansion, the scattering lengths are a ).These leading-order scattering lengths are a has been calculated using lattice QCD [41,42].The extrapolation to the physical pion mass gives a (3/2) πD = −0.16(4)fm [41].The scattering length a (1/2) πD has been calculated using lattice QCD [43,44], but it has not been extrapolated to the physical pion mass.
There have been several calculations of the S-wave πD scattering lengths beyond leading order in the chiral expansion [42,[45][46][47].The S-wave scattering lengths for all the pseudoscalar light mesons and pseudoscalar charm mesons through NLO in the chiral expansion are given in Ref. [45].The canceling factor (1−1) in the S-wave self energy at LO in Eq. ( 26) is replaced at NLO by (h 1 − 3h 3 − 6h ′ 5 )m π /M , where h 1 , h 3 , and h ′ 5 are dimensionless lowenergy constants.The constant h 1 is determined by charm-meson mass splittings, but h 3 and h ′ 5 have not been determined definitively.The calculations beyond leading order have used various prescriptions to fix the unknown low-energy constants.The resulting factors replacing (1 − 1) in Eq. ( 26) range from −0.36 to +0.67 [42,[45][46][47].The largest values were obtained using unitarized chiral perturbation theory constrained to reproduce the mass of the strange charm meson D * s0 (2317) [42,45].All the results for the S-wave contribution to the D thermal mass shift are smaller than the I = 1 2 term at LO, which is n π /f 2 π .
FIG. 5. Diagrams for the D * self energy from pion forward scattering.The last two diagrams are summed over the two directions for routing the pion momentum.The contribution to −i Π mn * (E, p) is the sum of the diagrams with D * legs and pion legs amputated, weighted by f π (ω q )/(2ω q ), integrated over q, and summed over the pion flavor i.

C. Vector charm meson
The complete propagator for a vector charm meson D * a with energy E (relative to its kinetic mass M * ) and momentum p has the form on the left side of Eq. ( 8), with rest energy ε * a , kinetic mass M * , and self energy Π * b (E, p).The behavior of the denominator near the pole in E has the form where δε * a is the thermal rest energy, δZ * a is the change in the residue of the pole, and ζ * a determines the change in the kinetic mass.The constants δε * a , δZ * a , and ζ * a can be identified as corrections to parameters of ZREFT.At leading order in the pion interactions, the thermal contribution to the D * self energy comes from pion forward scattering through the diagrams in Fig. 5.The contributions from the first diagram cancel upon summing over the pion flavor.The self energy tensor from the last two diagrams is The terms ω adq , ∆ ad , and ω acq in the denominator of the self-energy tensor in Eq. ( 29) are order m π .The heavy-meson expansion is obtained by expanding in powers of the small energies E − ε * a , E − ε * c , (p ± q) 2 /(2M ), and (p ± q) 2 /(2M * ).After averaging over the directions of the pion momentum, terms with odd powers of q vanish and the self-energy tensor reduces to the diagonal tensor Π * a (E, p) δ mn .The LO term is a constant independent of E and p: where F * ad is the complex conjugate of the thermal average defined in Eq. ( 18).The NLO term is a linear function of E and p 2 : where G * n,ad is the complex conjugate of the thermal average defined in Eq. (19).We can read off the corrections δε * a , δZ * a , and ζ * a defined by the expansion in Eq. ( 28):

D. Charm-meson Rest Energies
The complex thermal rest energies of the charm mesons in Eqs.(25a) and (32a) can be expressed as double expansions in isospin splittings divided by m π and in pion/charm-meson mass ratios using results given in Appendix A. The leading terms in the thermal mass shifts for the D and D * do not depend on their flavors:  The temperature dependence of the real parts of the charm-meson rest energies is illustrated in Fig. 6.The range of the temperature T is only up to 130 MeV, which may already be beyond the range of validity of our results.A solid curve is the sum of the rest energy at T = 0 and the thermal mass shift.The thermal mass shifts are smaller than the isospin splittings in this temperature range.The thermal mass shifts for D + and D 0 increase monotonically as T increases.The thermal mass shifts for D * + and D * 0 have local minima near 130 and 120 MeV, respectively.A dashed curve in Fig. 6 uses the approximation for the thermal mass shift in Eqs.(33).The approximations for D and D * in Eqs.(33a) and (33b) match well with the more accurate results in Eqs.(25a) and (32a) at sufficiently low temperatures.The approximation for D has the same qualitative behavior as the more accurate result at higher temperatures.The approximation for D * has the same qualitative behavior as the more accurate result only at temperatures lower than about 100 MeV.The difference at higher temperatures comes from the increasingly significant contributions from the thermal averages G * 2,ad and q 4 /ω 3 acq in Eq. (32a) at higher temperatures.The imaginary parts of the charm-meson thermal rest energies in Eqs.(25a) and (32a) can be evaluated analytically.The leading terms in their expansions in isospin splittings divided by m π and in pion/charm-meson mass ratios can be expressed in terms of the partial decay rates for D * → Dπ: The leading terms in the thermal widths of the charm mesons are The widths for the charm mesons as functions of the temperature T are shown in Fig. 7.For a pseudoscalar charm meson D b , the solid curve is the thermal width δΓ b .For a vector charm meson D * a , the solid curve is the sum Γ * a + δΓ * a of the decay width and the thermal width.The thermal widths all increase monotonically with T .A dashed curve in Fig. 7 uses the approximation for the thermal width in Eqs.(35).The approximations for D and D * in Eqs.(35a) and (35b) match well with the more accurate results in Eqs.(25a) and (32a) at sufficiently low temperature.The approximations have the same qualitative behavior as the more accurate results at higher temperatures.

E. Comparison with Previous Work
There have been several previous calculations of the thermal energies of charm mesons in a pion gas [48][49][50][51][52].The temperatures ranged from 0 to beyond the hadronization temperature, which is near 150 MeV.Our results are not applicable at temperatures as high as 150 MeV.We will therefore only compare our results with previous calculations at 100 MeV.The comparisons are presented in Table I.
The complex thermal energies of the charm mesons D and D * in a pion gas have been calculated by Fuchs et al. using a resonance approach [48].Their πD ( * ) elastic scattering amplitudes are the sums of relativistic Breit-Wigner amplitudes corresponding to the known S-wave, P-wave, and D-wave πD ( * ) resonances.These amplitudes are not consistent with the low-energy constraints from chiral symmetry.The thermal self energies from pion forward scattering were obtained by integrating the πD ( * ) scattering amplitudes over the pion momentum distribution.The results of Fuchs et al. for the D and D * thermal mass shifts and the D and D * thermal widths at T = 100 MeV are given in Table I.Similar results for the D thermal width were obtained by He et al. using a resonance approach [49].
The thermal widths of D and D * in a pion gas have been calculated by Cleven et al. using a self-consistent unitarized approach [50].Their tree-level πD ( * ) amplitudes are vertices from an effective lagrangian for pseudoscalar and vector mesons with an approximate SU (4) symmetry.They did not include diagrams with intermediate D or D * propagators, like those in Figs. 4 and 5. Their unitarized πD ( * ) amplitudes are the analytic solutions of the corresponding regularized Lippmann-Schwinger equations.A thermal D ( * ) self energy can be obtained numerically as a function of the 3-momentum by integrating the πD ( * ) amplitude over the loop momentum and over the pion momentum distribution.A Lippmann-Schwinger integral equation for the thermal πD ( * ) amplitude can be obtained by inserting the thermal D ( * ) self energy into the D ( * ) propagator in the loop.Self-consistent solutions for the thermal D ( * ) self energy and the thermal πD ( * ) amplitude were obtained by iterating the equations several times numerically.Their results for the thermal mass shifts of D and D * were little more than numerical noise.The results of Clevens et al. for the D and D * thermal widths at T = 100 MeV are given in Table I.
The complex thermal energies of D in a pion gas and in a gas of π, K and η mesons have been calculated by Montaña et al. using a different self-consistent unitarized approach [51].Their low-energy πD scattering amplitude is the sum of the contact vertex in Fig. 4 and the next-to-leading order vertices in heavy-hadron χEFT.They did not include the diagram with the intermediate D * propagator in Fig. 4. Their unitarized amplitude is the analytic solution of the corresponding regularized Lippmann-Schwinger equation.The unitarization reduces the accuracy to leading order in the chiral expansion.Self-consistent solutions to the integral for the thermal D self energy and the Lippmann-Schwinger integral equation for the thermal πD amplitude were obtained numerically.The results were extended to other charm mesons including D * in Ref. [52].The results of Montana et al. for the D and D * thermal mass shifts and the D and D * thermal widths at T = 100 MeV are given in Table I.
The thermal mass shifts and thermal widths for D and D * at T = 100 MeV calculated in Refs.[48][49][50][51][52] differ at most by about a factor of 2. Our results in Table I for the flavor averages at T = 100 MeV of the thermal mass shifts for D from Eq. (25a) and for D * from Eq. (32a) are +1.07MeV and −0.15 MeV, which have opposite signs.The thermal mass shifts for D and D * in Refs.[51,52]  for the flavor averages at T = 100 MeV of the thermal widths for D from Eq. (25a) and for D * from Eq. (32a) are 64 keV and 15 keV.The thermal widths for D and D * in Refs.[51,52] are approximately equal and larger than our results by about 3 orders of magnitude.
The orders-of-magnitude discrepancies between our results and those in Refs.[51,52] are a serious issue.We cannot exclude the possibility that higher-order terms in the chiral expansions of the thermal rest energies for charm mesons are much larger than the leadingorder terms, even at temperatures as low as 100 MeV.In this case, our thermal masses and thermal widths for charm mesons could be much too small.One possible source of large corrections is from S-wave pion interactions, since the contributions from total isospin 1 2 and 3 2 cancel at leading order.All previous calculations of S-wave πD scattering lengths at higher order in the chiral expansion, including those obtained using a chiral unitarization prescription, imply that the S-wave contribution to the D thermal mass shift is smaller than the I = 1 2 contribution at leading order, which can be obtained from Eq. ( 26) by replacing the factor (1-1) by 1.For a pion gas at T kf = 100 MeV, the I = 1 2 contribution is 4.3 MeV, which is about 4 times larger than the flavor-averaged thermal mass shift from P-wave interactions.The D mass shift in Ref. [51] is larger by another factor of 3, and it has the opposite sign.So this possible source of large corrections seems to be insufficient.Another possible explanation for the large discrepancies between our results and those in Refs.[51,52] is that they arise from the iteration of the integral equation used to obtain self-consistent solutions.This integral equation was obtained by generalizing the integral equation for a chiral unitarization prescription to nonzero temperature.The chiral unitarization prescription is just a model, and the proposed generalization to finite temperature may give thermal mass shifts and thermal widths for charm mesons that are much too large at low temperatures.In Refs.[51,52], the S-wave contributions to πD scattering were taken into account, but the P-wave contributions were ignored.This is appropriate only at temperatures much smaller than m π .The effects of the P-wave contributions could be significant at T = 100 MeV.
The first lattice QCD calculations of charm meson masses as functions of the temperature were presented in Ref. [53].The lattice QCD simulations were carried out with 3 flavors of dynamical light quarks and a pion mass of 239 MeV.The results have not been extrapolated to the physical pion mass.The masses of D and D * were calculated at 6 temperatures ranging from 47 to 169 MeV.At the 4 lowest temperatures, which range from 47 MeV to 127 MeV, the thermal shifts in the masses of D and D * are consistent with zero to within the errors, which are about 6 MeV.The lowest temperature where the thermal mass shifts are significant is 152 MeV, where the shift is −20 ± 7 MeV for D and −43 ± 10 MeV for D * .The D thermal mass shift at 150 MeV in Refs.[51,52] is larger by a factor of 2, which is about 3 standard deviations, while the D * thermal mass shift is compatible with Ref. [53].At a temperature of 109 MeV, the thermal shifts in Ref. [53] are 0±6 MeV for D and +4±6 MeV FIG. 8. One-loop diagrams for the D * D self energy from pion forward scattering.The diagrams are summed over the two directions for the routing of the pion momentum.The contribution to i(µ/2π)Σ is the sum of the diagrams with pion legs amputated, weighted by f π (ω q )/(2ω q ), integrated over q, and summed over the pion flavor i.The diagrams with a πD contact vertex or a πD * contact vertex are not shown, because each diagram cancels upon summing over the pion flavor.
for D * .The D and D * thermal shifts at 100 MeV in Refs.[51,52] are outside the error bar by 2.2 and 2.7 standard deviations, respectively.Thus lattice QCD is beginning to provide useful results for the thermal shifts in charm-meson masses.

V. CHARM-MESON PAIR SELF ENERGY
In this Section, we calculate the D * D self energy in the pion gas to NLO in the heavymeson expansion.For simplicity, we calculate it only for charm-meson pairs with kinetic energies of order m 2 π /M .The self-energy is complicated by terms that diverge at the charmmeson-pair threshold and by terms that correspond to ultraviolet divergences in ZREFT.

A. Leading Order
At leading order in the pion interactions, the thermal contributions to the D * D selfenergy come from pion forward scattering.Four of the corresponding diagrams for the D * D self energy are shown in Fig. 8.There are two additional diagrams with πD or πD * contact vertices.They are not shown, because they cancel upon summing over the pion flavor.Each pion-forward-scattering diagram is the sum of two terms: an on-shell pion emerges from one of the two open circles with momentum q and flavor i and it scatters into the other open circle with the same momentum and flavor.The amplitudes must be added coherently by multiplying them by f π (ω q )/(2ω q ), integrating over q with measure d 3 q/(2π) 3 , and summing over the flavors i.
The self energy diagrams in Fig. 8 are functions of the total energy E and total momentum P of the D * D pair.They are integrated over the pion momentum q, the loop momentum k, and the loop energy ω.The integral over ω can be evaluated by contours.In the first three diagrams in Fig. 8, we use the pole of the single propagator connecting the D * D pair vertices.The large pion momentum of order m π must then flow through the single propagator connecting the pion vertices.The other two propagators have identical denominators with small energies of order m 2 π /M .In the fourth diagram in Fig. 8, we integrate over ω using the poles of the two adjacent propagators on the lower side of the loop.The large pion momentum of order m π must flow through the two propagators on either the right side or the left side of the loop.The remaining propagator on the upper left or upper right of the loop has a denominator with a small energy of order m 2 π /M .In a charm-meson propagator carrying the large momentum q, the large energies in the denominator are in terms of the form ω cdq and ∆ cd , which are order m π .The heavy-meson expansion is obtained by expanding in powers of small energies of order m 2 π /M , which include the energy difference E − ε * a − ε b , the kinetic energies of the charm mesons, and isospin splittings.Only the first three diagrams in Fig. 8 have contributions at LO in that expansion.All three diagrams factor into an integral over q and the same convergent integral over k, which is given in Eq. (B1) of Appendix B. The LO contributions from the third diagram cancel upon summing over the two directions for the routing of the pion 4-momentum, so it is actually NLO.In the first two diagrams, the integrand of the integral over q can be simplified by summing over the two directions for the routing of the pion 4-momentum.The contributions to the D * D self energy at LO from the first two diagrams in Fig. 8 are where F cd is the thermal average defined in Eq. ( 18) and S 0 (E cm ) is the function of the center-of-mass energy defined in Eq. ( 2).
The LO contribution to the D * D self energy is the sum of Eqs.(36).It can be expressed as µ(δε * a + δε b )/S 0 (E cm ), where δε b and δε * a are the thermal rest energies of the charm mesons at LO, which are given by the F cb term in Eq. (25a) and by the F * ad term in Eq. (32a).This agrees with the LO term in the expansion of S 1 (E cm , p) around S 0 (E cm ) in Eq. ( 14).Thus the D * D self energy at LO can be completely absorbed into the complex thermal rest energies of the charm-meson constituents.

B. Next-to-Leading Order from One-loop Diagrams
The D * D self energy has NLO contributions from all four diagrams in Fig. 8.The NLO contributions from the first three diagrams are obtained by expanding the propagator carrying the large momentum q to first order in the small energies.The resulting integrals over the loop momentum can be reduced to the forms found in Appendix B: Eq. (B1), which gives a factor of 1/S 0 (E cm ), and Eqs.(B2b) and (B2c), which are ultraviolet divergent.The NLO contribution from the fourth diagram in Fig. 8 can be obtained by setting the small energies in the two propagators carrying the large momentum equal to 0. The resulting integral over the loop momentum reduces to Eq. (B2a), which is ultraviolet divergent.The NLO terms in the four self-energy diagrams are where Λ is the ultraviolet cutoff, ∆ = M * − M , and G n,cd and H ab are the thermal averages defined in Eqs. ( 19) and (20).There are terms in the D * D self energy at NLO that diverge at the branch point of S 0 (E cm ).We compare them with the terms expected from the expansion of the square-root function S 1 (E cm , P ) in Eq. ( 14).The sum of the singular terms proportional to 1/S 0 (E cm ) in Eq. ( 37) matches the contribution to the term µ(δε * + δε)/S 0 (E cm ) in Eq. ( 14) from the NLO thermal D b and D * a rest energies δε b and δε * b in Eqs.(25a) and (32a).The sum of the singular terms proportional to P 2 /S 0 (E cm ) in Eq. ( 37) can be expressed in the form of the last term in Eq. ( 14) with the coefficient An expression for ζ X that takes into account propagator corrections for the constituents was deduced in Eq. (12).Its expansion to first order in ζ and ζ * gives Upon inserting the NLO expressions for ζ b and ζ * a in Eqs.(25c) and (32c), we reproduce the expression for ζ X in Eq. (38).Thus all the terms in the D * D self energy at NLO that diverge at the branch point of S 0 (E cm ) can be absorbed into the square-root function S 1 (E cm , P ) defined in Eq. (10).

C. Next-to-Leading Order from a Two-loop Diagram
The only two-loop diagram for the D * D self energy at NLO is shown in Fig. 9.The pion momentum can be routed through two D * propagators and through the contact vertex that attaches the two loops.The Feynman rule for the D * D contact vertex in the ab channel is i C 0 δ mn , where m and n are the vector indices of the attached D * lines.The NLO contribution from this diagram is obtained by setting the small energies to 0 in the two charm-meson propagators carrying the pion momentum.The integral over the loop momentum can be reduced to Eq. (B2a), which is ultraviolet divergent.The resulting contribution to the D * D self energy is If we use Eq. ( 4) to replace C 0 /(2π/µ) by 1/(Λ − γ X ) and then take the large Λ limit, the 2-loop contribution to the D * D self energy at NLO reduces to

D. Contact vertex corrections
The ultraviolet divergent part of the D * D self energy at NLO is the sum of the terms proportional to Λ in Eqs.(37) and Eq. ( 41): These ultraviolet divergences must be canceled by ultraviolet divergences from corrections to the contact vertex.The contact vertex C 1 required to compensate for D and D * propagator corrections at short distances as well as short-distance vertex corrections is given in Eq. ( 13).The expansion of (2π/µ)/C 1 to NLO is  25c) and (32c).The δC term allows for additional corrections to the contact vertex that might be needed to compensate for short-distance vertex corrections.At NLO in the heavy-meson expansion, there are D * D vertex corrections from pion forward scattering from the diagrams in Fig. 10.The diagrams have two charm-meson propagators that carry a large energy of order m π that flows through the D * D contact vertex.The NLO contribution from those diagrams is obtained by setting the small energies of order m 2 π /M in the denominators of those propagators equal to 0. The changes in the D * D contact vertex required to cancel the vertex corrections from the first two diagrams in Fig. 10 are The change in the D * D contact vertex required to cancel the vertex correction from the third diagram in Fig. 10 is Upon using Eq. ( 4) to replace the factor of (2π/µ)/C 0 in Eq. ( 43) by Λ − γ X , we obtain terms linear in the ultraviolet cutoff Λ: The terms proportional to Λ cancel the ultraviolet divergences in Eqs.(42).
The denominator of the complete D * D propagator has the form where C 1 is the complete contact-interaction vertex and Σ(E cm , P ) is the D * D self energy.The D * D self energy at NLO is the sum of Eqs. ( 36), (37), and (41).The ultraviolet divergences in the D * D self energy are canceled by the corrections from the contact interaction vertex C 1 in Eq. ( 46).The net effect is essentially to replace Λ in the D * D self energy by γ X .The terms in the D * D self energy proportional to 1/S 0 (E cm ), which diverge at the branch point of S 0 (E cm ), can be absorbed by replacing S 0 (E cm ) in Eq. ( 47) by the square-root function S 1 (E cm , P ) in Eq. ( 10).An appropriate resummation of higher order corrections would replace the remaining S 0 (E cm ) terms in the D * D self energy by S 1 (E cm , P ).The resulting expression for the denominator of the D * D propagator at NLO is

VI. THERMAL ENERGY OF THE MOLECULE
In this Section, we use the D * D self-energy at NLO in the heavy-meson expansion to determine the thermal correction to the binding momentum of a loosely bound charm-meson molecule in the pion gas.We then calculate the thermal mass shifts and thermal widths of X(3872) and T + cc (3875) in the pion gas.

A. Correction to the binding momentum
The denominator of the D * D propagator through NLO in Eq. ( 48) can be expressed in the form of the denominator in Eq. ( 7) by factoring out the coefficient of S 1 into a multiplicative factor Z −1 X given by After factoring out Z −1 X and re-expanding to NLO, the denominator of the D * D propagator reduces to the much simpler expression This has the same form as the denominator of the amplitude for a loosely bound molecule in Eq. ( 5) but with S 0 (E cm ) replaced by the thermally modified square-root function S 1 (E cm , P ) in Eq. ( 10) and with a correction δγ X to the binding momentum.The NLO thermal correction to the binding momentum is The residue factor Z −1 X in Eq. ( 49) can be absorbed into the normalization of the local composite operator whose 2-point Green function is the D * D propagator.The NLO corrections in Z −1 X therefore have no physical effects.The remaining NLO corrections to the D * D self energy have been completely absorbed into thermal corrections to parameters of ZREFT.The specific parameters are those that determine the propagators of the individual charm mesons at NLO, which appear in the square-root function S 1 (E cm , P ) in Eq. ( 10), and the complex binding momentum, whose thermal correction is given in Eq. (51).
The thermal averages in Eq. ( 51) can be expressed as double expansions in isospin splittings divided by m π and in pion/charm-meson mass ratios.The leading term in the real part of the fractional change in the binding momentum is Re [For a pion gas at T kf = 115 MeV, its value is Re[δγ X /γ X ] = −3.4×10−4 .The more accurate results from Eq. ( 51) differ by −11% for X and by −15% for T + cc .The leading term in the imaginary part of the fractional change in the binding momentum is where q cd = ∆ 2 cd − m 2 πcd is the momentum of the pion in the decay D * c → D d π.This expression is different for X(3872) with constituents D * 0 D0 /D 0 D * 0 and for T + cc (3875) with constituents D * + D 0 .For X, the linear combination of pion momenta in Eq. ( 53) is (4q 00 + 6q +0 )/3.For a pion gas at T kf = 115 MeV, the imaginary part of the fractional change in the binding momentum is Im[δγ X /γ X ] = 2.8 × 10 −4 .For T + cc , the linear combination of pion momenta in Eq. ( 53) is (3q 00 +q ++ +8q +0 )/3.For a pion gas at T kf = 115 MeV, the imaginary part of the fractional change in the binding momentum is Im[δγ T /γ T ] = 3.4 × 10 −4 .The more accurate results for Im[δγ X /γ X ] from Eq. ( 51) differ from Eq. ( 53) by about −8% for both X and T + cc .The only thermal average in the inverse propagator in Eq. ( 48) that does not have a well defined limit as ϵ → 0 + is H ab , which has a term that diverges as 1/ϵ in this limit.A more careful treatment of this thermal average would regularize the divergence by replacing ϵ by appropriate decay widths of charm mesons, as in Eq. (A9) of Appendix A. The thermal average H ab is therefore sensitive to these decay widths.Since H ab appears in the residue factor Z −1 X in Eq. ( 49) but not in the denominator of the D * D propagator in Eq. ( 50), it does not affect the energy-momentum relation for X.We therefore do not give any analytic approximation for H ab .

B. Thermal mass shift and thermal width
The zero of the inverse propagator in Eq. ( 50) determines the energy-momentum relation for X.The pole energy of X with zero 3-momentum at NLO can be expressed as where δγ X is the NLO correction to the binding momentum in Eq. ( 51).We have taken into account the decay width of D * a at T = 0 by replacing iϵ in Eq. ( 10) with iΓ * a /2.For X(3872), the pole energy at T = 0 relative to the real D * 0 D0 threshold determined in Ref. [13] is (0.025 − 0.140 i) MeV.In a pion gas at T kf = 115 MeV, the pole energy relative to the real D * 0 D0 threshold from Eq. ( 54) is (+1.64−0.21i) MeV.The change in the pole energy comes primarily from the complex thermal energy shifts of the D * 0 D0 /D 0 D * 0 constituents, whose sum is (+1.61 − 0.07 i) MeV.The contribution from the complex thermal correction δγ X to the binding momentum is (+0.04 + 0.08 i) keV.Its real and imaginary parts are smaller than those from the energy shifts of the constituents by orders of magnitude.
For T + cc (3875), the pole energy at T = 0 relative to the real D * + D 0 threshold determined in Ref. [15] is (−0.36−0.024i) MeV.In a pion gas at T kf = 115 MeV, the pole energy relative to the real D * + D 0 threshold from Eq. ( 54) is (+1.20 − 0.10 i) MeV.The change in the pole energy comes primarily from the complex thermal energy shifts of the D * + D 0 constituents, whose sum is (+1.56 − 0.08 i) MeV.The contribution from the complex thermal correction δγ T to the binding momentum is (+0.20 − 0.23 i) keV.Its real and imaginary parts are smaller than those from the energy shifts of the constituents by orders of magnitude.
The temperature dependence of the pole energies of X(3872) and T + cc (3875) is illustrated in Fig. 11.The range of the temperature T extends only up to 130 MeV, which may already be beyond the range of validity of our results.In the left panel of Fig. 11, a solid curve is the real part of E X , which is the sum of the rest energy of the molecule at T = 0 and its thermal mass shift.A dashed curve is the T -dependent charm-meson pair threshold.The thermal mass shifts of X and T + cc increase monotonically with T .Their binding energies relative to the charm-meson pair thresholds in the pion gas are the differences between the dashed and solid curves.The T -dependence of the binding energies, which comes from the thermal correction to the binding momentum, is too small to be visible in Fig. 11.In the right panel of Fig. 11, a solid curve is the imaginary part of E X multiplied by −2, which is the sum of the decay width of the molecule and its thermal width.The thermal widths increase monotonically with T .Their dependence on T comes almost entirely from the thermal rest energies of the charm-meson constituents.The remaining T -dependence from the complex thermal correction to the binding momentum is too small to be visible in Fig. 11., and their constituents in the pion gas as functions of the temperature T .In the left panel, the solid curves are the sums of the rest energy at T = 0 and the thermal mass shift for X (red lower curve) and T + cc (blue upper curve).The dashed curves and the horizontal lines are the corresponding charm-meson pair thresholds in the pion gas at temperature T and at T = 0, respectively.In the right panel, the solid curves are the sums of the decay width at T = 0 and the thermal width for X (red upper curve) and T + cc (blue lower curve).The horizontal lines are the decay widths of X (red upper line) and T + cc (blue lower line) at T = 0 MeV.The dashed curves are the sums of the widths of the constituents D * 0 and D 0 of X (red lower curve) and the constituents D * + and D 0 of T + cc (blue upper curve).

C. Comparison with Previous Work
There have been two previous calculations of the thermal energies of X(3872) in a pion gas [39,40].The temperatures ranged from 0 to 150 MeV.Our results are not applicable at temperatures as high as 150 MeV, which is near the hadronization temperature.We will therefore only compare our results with the previous calculations at 100 MeV.The comparisons are presented in Table II.
The effect of a hot pion bath on the X(3872) has been studied previously by Cleven, Magas, and Ramos (CMR) [39].The X was identified with a dynamically generated bound state with energy 2.5 MeV below the D * D threshold.This energy is the pole in a unitarized DD * amplitude given by the analytic solution of a regularized Lippmann-Schwinger equation.The DD * tree amplitude in that equation is the vertex from an effective Lagrangian for pseudoscalar and vector mesons with an approximate SU (4) symmetry, as in Ref. [50].The thermal D * D amplitude T DD * is obtained by the numerical solution of the Lippmann-Schwinger integral equation with charm meson propagators that include the self-consistent thermal D ( * ) self energies calculated in Ref. [50].The energy of X and its width were identified with the energy at the peak of |T DD * | and the width of the peak, respectively.As T increases from 0 to 100 MeV, the energy of X increases from −2.5 to +3 MeV.Although it was not stated in Ref. [39], the D * D threshold was actually held constant at its T = 0 value, according to Ref. [40].The thermal energy of X in Ref. [39] should therefore be interpreted as the difference between its thermal energy and an unknown T -dependent shift in the D * D threshold.As T increases from 0 to 100 MeV, the width of X increases from 0 to 30 MeV.The effect of a hot pion bath on the X(3872) has also been studied previously by Montaña, Ramos, Tolos, and Torres-Rincon (MRTT) [40].They also considered its effect on a possible P-wave charm-meson molecule called X(4140).The X was identified with a dynamically generated bound state with energy 4.1 MeV below the D * D threshold.This is the energy of a pole in a unitarized DD * amplitude given by the analytic solution of a regularized Lippmann-Schwinger equation.The DD * tree amplitude in that equation is obtained from the effective Lagrangian for pseudoscalar and vector mesons with an approximate SU (4) symmetry used in Ref. [52].The thermal D * D amplitude T DD * is obtained from the numerical solution of the Lippmann-Schwinger integral equation with charm meson propagators that include the self-consistent thermal D ( * ) self energies calculated in Ref. [52].The energy of X and its width were identified with the energy of the peak in −Im[T DD * ] and the width of the peak, respectively.As T increases from 0 to 100 MeV, the energy of X decreases from −4 to −30 MeV.Most of the change comes from the thermal mass shifts of the charm-meson constituents, whose sum decreases from 0 to −27 MeV.The energy of X relative to the Tdependent charm-meson-pair threshold increases from −4 to −3 MeV.As T increases from 0 to 100 MeV, the width of X increases from 0 to 30 MeV.Most of the change comes from the thermal widths of the charm-meson constituents, whose sum increases from 0 to 34 MeV.The thermal width agrees with that in Ref. [39] at 100 MeV.
We proceed to compare our results for X(3872) with those in MRTT for the pion gas at T = 100 MeV.Our result in Table II for the sum of the thermal mass shifts of the constituents D * 0 and D0 from Eqs. (25a) and (32a) is +1.00 MeV.This is an order of magnitude smaller than the sum of the thermal mass shifts of D * and D in Ref. [52] and it has the opposite sign.Our result in Table II for the real part of the pole energy of X relative to the complex charm-meson-pair threshold from Eq. ( 54) is +25 keV at T = 0 and it increases by only 0.03 keV at T = 100 MeV.The corresponding energy at T = 100 MeV in MRTT is about two orders of magnitude larger and it has the opposite sign.Our result in Table II for the sum of the thermal widths of the constituents D * 0 and D0 from Eqs. (25a) and (32a) is 0.11 MeV.This is smaller than the sum of the thermal widths of D * and D in Ref. [52] by a factor of about 300.Our result in Table II for the difference between the width of X from Eq. ( 54) and the sum of the widths of the constituents D * 0 and D0 is 225 keV at T = 0 and it decreases by only 0.11 keV at T = 100 MeV.The corresponding difference in MRTT is 0 at T = 0 and it decreases to −4 MeV at T = 100 MeV.
The orders-of-magnitude discrepancies between our results for X and its constituents and those in Ref. [40] are a serious issue.We cannot exclude the possibility that higher orders in the chiral expansion are much larger than the leading order, even at temperatures as low as 100 MeV.Another possibility is that the proposed generalizations to nonzero temperature of the integral equations in Ref. [40] may give thermal mass shifts and thermal widths for the charm-meson molecule that are much too large at low temperatures.

VII. SUMMARY
We have calculated the thermal energy of a loosely bound charm-meson molecule in a pion gas.The molecule consists of a bound pair of charm mesons in a channel we denoted by D * D. It is associated with a pole in the D * D propagator very close to the D * D threshold.At zero temperature, the molecule can be described by a ZREFT for nonrelativistic charm mesons with a nonperturbative contact interaction in the D * D channel.At leading order in the ZREFT, the D * D propagator has the simple form in Eq. ( 5).We assumed that the pole in the D * D propagator and the branch cuts associated with the D * D threshold remain the only nearby singularities even in the thermal environment provided by the pion gas.The temperature of the pion gas was assumed to be low enough that the interactions between charm mesons and pions can be described by heavy-hadron χEFT at leading order.
The D * D self-energy was calculated to NLO in the heavy-meson expansion.At LO in that expansion, the D * D self-energy is the sum of the contributions in Eqs.(36) from two one-loop diagrams.Both terms diverge at the D * D threshold, but these diverging terms can be absorbed into the parameters of ZREFT by adding the LO thermal corrections δε and δε * to the charm-meson rest energies ε and ε * in the square-root function S 0 (E cm ) in Eq. (2).Thus there are no thermal corrections to the binding energy of the mloecule at this order.
At NLO in the heavy-meson expansion, the D * D self-energy is much more complicated.It is the sum of the contributions in Eqs.(37) from one-loop diagrams and in Eq. ( 41) from a two-loop diagram.The D * D self-energy at NLO has linear ultraviolet divergences as well as terms that diverge at the D * D threshold.The terms that diverge at the D * D threshold can be absorbed into the parameters of ZREFT by replacing S 0 (E cm ) in the D * D propagator by the thermally corrected square-root function S 1 (E cm , P ) in Eq. (10), which takes into account thermal corrections to the charm-meson rest energies and their kinetic masses.The ultraviolet divergences in the D * D self-energy at NLO are canceled by the correction to the strength of the contact vertex in Eq. ( 46), which compensates for thermal corrections to the charm-meson propagators as well as thermal vertex corrections.
After canceling the ultraviolet divergences and absorbing the divergences at the D * D threshold into charm-meson propagator corrections, the denominator of the complete D * D propagator at NLO reduces to Eq. ( 48).This form is consistent with our assumption that the molecule in the pion gas can be described by ZREFT.After factoring out the residue factor Z −1 X in Eq. ( 49), the denominator of the complete D * D propagator at NLO reduces to Eq. (50).Aside from the thermal corrections to the charm-meson propagators that appear in the square-root function, the only other thermal correction is the NLO correction to the binding momentum in Eq. (51).It is proportional to the pion number density n π , and it is suppressed by a factor of ∆/M X ≈ m π /M X .The thermal mass shift and thermal width of the molecule is therefore dominated by the thermal mass shifts and thermal widths of the charm-meson constituents.
The thermal energy of the X(3872) in a pion gas has been calculated previously in Ref.
[39] (CMR) and Ref. [40] (MRTT).In MRTT, the thermal mass shift and the thermal width of X(3872) were calculated as a function of T up to 150 MeV.Our results are certainly not valid at such a high temperature, but we can compare results near the kinetic freeze-out temperature.At such a temperature, the thermal mass shift for X comes almost entirely from the thermal mass shifts of the D * D constituents.The sum of our thermal mass shifts for D * and D is positive while the sum in MRTT is negative.At T = 100 MeV, the sum in MRTT is more than an order of magnitude larger in absolute value than ours.At a temperature near kinetic freezeout, the thermal contribution to the width of X comes almost entirely from the thermal widths of the D * D constituents.At T = 100 MeV, the sum of the thermal widths for D * and D in MRTT is more than two orders of magnitude larger than ours.
Since the results for X in CMR and MRTT are completely numerical, it is difficult to identify the reasons for the orders-of-magnitude discrepancies from our results.We cannot exclude the possibility that thermal chiral EFT is not applicable even at temperatures as low as 100 MeV, in which case our thermal mass shift and thermal width for the loosely bound molecule could be much too small.Another possibility is that the iteration of the integral equations used to obtain self-consistent thermal D ( * ) self energies in Refs.[39,40] give thermal mass shift and thermal width for X that are much too large at low temperatures.
We have calculated the D * D self energy to NLO in the heavy-meson expansion.At nextto-next-to-leading order (NNLO) in the heavy-meson expansion, there will be additional terms that diverge at the D * D threshold and additional ultraviolet divergences.It would be interesting to know whether all the NNLO terms that diverge at the D * D threshold can be absorbed into D * and D propagator corrections and whether all the ultraviolet divergences at NNLO can be cancelled by renormalizations of parameters of ZREFT.The NLO correction to the binding momentum is suppressed by a factor of ∆/M .If the NNLO correction to the binding momentum has no such suppression factor, it could be comparable to the NLO correction.
We have assumed the interactions of charm mesons and pions in the hadron gas can be described by χEFT at leading order at temperatures at least as high as that of kinetic freeze-out.Some insight into the applicability of thermal χEFT could be provided by explicit next-to-leading order calculations.There have been very few calculations in thermal χEFT beyond leading order.The thermal energy of the pion has been calculated to next-to-leading order by Schenk [54].Thermal corrections to the pion decay constant and the pion mass have been calculated to next-to-leading order by Toublan [55].The corrections are small enough to justify optimism towards the applicability of χEFT at temperatures as high as 115 MeV.A next-to-leading order calculation of the thermal energy of a charm meson should reveal whether thermal χEFT is applicable at that temperature.It could also provide a test of the validity of the self-consistent iteration methods used in Refs.[39,40].
Our results suggest that loosely bound charm-meson molecules, such as X(3872) and T + cc (3875), can remain loosely bound and narrow in the thermal environment of a hadron gas at sufficiently low temperature.This is consistent with the surprisingly large rate for the production of X(3872) in Pb-Pb collisions at the LHC observed by the CMS collaboration [12].Our results are more encouraging for the study of loosely bound charm-meson molecules in heavy-ion collisions than the previous results for X(3872) in Refs.[39,40].Further studies of X(3872) along with studies of T + cc (3875) in heavy-ion collisions will provide essential insights into the behavior of loosely bound hadronic molecules.
where ε a = M a − M is the rest energy of D a .The propagator for a nonrelativistic vector charm meson D * with energy E (relative to the kinetic mass M * ) and 3-momentum p is iδ ab δ mn E − ε * a − p 2 /(2M * ) + iϵ , (C3) where ε * a = M * a − M * is the rest energy of D * a .The D * decay width can be taken into account by replacing iϵ with iΓ * a /2.

Vertices
In heavy-hadron χEFT at leading order, the vertices for the interactions between charm mesons and pions are determined by the pion decay constant f π and a dimensionless constant g π .The D ( * ) π → D ( * ) π contact vertices, with the incoming and outgoing pions having 3momenta q and q ′ , are D a π i (q) → D b π j (q ′ ): + i 4f 2 π [σ i , σ j ] ab (ω q + ω q ′ ), (C4a) where σ i is a Pauli matrix and ω q = m 2 π + q 2 .The vertices for the transitions D ( * ) → D ( * ) π, with the outgoing pion having 3-momentum q, are D a → D * b n π i (q) : +i g π √ 2f π σ i ab q n , (C5a) where ε mnk is the Levi-Civita symbol.The corresponding vertices for the transitions D ( * ) π → D ( * ) , with the incoming pion having 3-momentum q, are obtained by replacing q with −q.

FIG. 1 .
FIG. 1. Bubble diagram for the amplitude for the propagation of D and D * between contact interactions.The propagator for the D and D * are represented by a solid line and a double (solid+dashed) line, respectively.The Feynman rule for the open circles is 1.
* c and D d by M * c and M d .We decompose them into kinetic masses M * and M and rest energies ε * c and ε d : M * c = M * + ε * c and M d = M + ε d .Convenient choices for M and M * are the isospin averages of the charm-meson masses: M = (M + + M 0 )/2, M * = (M * + + M * 0 )/2.The rest energies ε * c and ε d are then comparable to isospin splittings.The charm-meson mass differences ∆ cd = M * c −M d differ from m π by amounts comparable to isospin splittings.The value of ∆ cd averaged over the four D * → D transitions is ∆ = M * − M = 141.3MeV.We

FIG. 3 .
FIG.3.One-loop diagrams for the D self energy in thermal field theory.
where δε b is the thermal rest energy, δZ b is the change in the residue of the pole, and ζ b determines the change in the kinetic mass.The constants δε b , δZ b , and ζ b can be identified as corrections to parameters of ZREFT.At leading order in the pion interactions, the thermal contribution to the D self energy comes from pion forward scattering through the diagrams in Fig.4.The contributions from the first diagram cancel upon summing over the pion flavor.In the contributions from the second diagram, the sum over the pion flavor i of the Clebsch-Gordan factors from the pion vertices can be evaluated using a simple identity for Pauli matrices: i σ i cd 2 = 2 − δ cd , where c and d are flavor indices for charm mesons.The self energy from the second diagram

FIG. 4 .
FIG.4.Diagrams for the D self energy from pion forward scattering.The second diagram is summed over the two directions for the routing of the pion momentum.The contribution to −i Π(E, p) is the sum of the diagrams with D legs and pion legs amputated, weighted by f π (ω q )/(2ω q ), integrated over q, and summed over the pion flavor i.

(1/ 2 )
πD = 0.233 fm and a (3/2) πD = −0.116fm.The scattering length a (3/2) πD where the sums are over the flavors of the intermediate D d and the intermediate D * c and ∆ ad is the D * a -D d mass difference.The angular brackets indicate the average over the Bose-Einstein distribution for the momentum q of the pion with flavor ad or ac.

)
For a pion gas at T kf = 115 MeV, these thermal mass shifts are Re[δε b ] = +1.26MeV and Re[δε * a ] = −0.42MeV.The more accurate results for D + and D 0 from Eq. (25a) differ by factors of 1.27 and 1.41, respectively.The more accurate results for D * + and D * 0 from Eq. (32a) differ by factors of 0.51 and 0.38, respectively.The spin-weighted average of the simple thermal mass shifts in Eq. (33) is zero.At T kf = 115 MeV, the spin-weighted average of the more accurate thermal mass shifts for the four charm mesons is +0.28 MeV.

FIG. 6 .
FIG.6.Real rest energies for D + and D 0 (left panel: blue upper and red lower curve) and for D * + and D * 0 (right panel: blue upper and red lower curve) in a pion gas as functions of the temperature T .The solid curve is the sum of the rest energy ε at T = 0 and the thermal mass shift Re[δε].The dashed curve uses the simple approximation for the thermal mass shift in Eqs.(33).The horizontal line is the rest energy at T = 0.

FIG. 7 .
FIG.7.Widths for D + and D 0 (left panel: blue lower and red upper curve) and for D * + and D * 0 (right panel: blue upper and red lower curve) in a pion gas as functions of the temperature T .A solid curve is the sum of the decay width Γ at T = 0 and the thermal width δΓ.A dashed curve uses the simple approximation for the thermal width in Eqs.(35).A horizontal line is the decay width at T = 0.

FIG. 9 .
FIG. 9. Two-loop diagram for the D * D self energy from pion forward scattering.The D * flavors are constrained by the vertices to be D * a , so the forward-scattered pion must be π 0 .

FIG. 10 .
FIG.10.Diagrams for vertex corrections to the D * D contact interaction from pion forward scattering.

FIG. 11 .
FIG.11.The real parts of the poles of the pair propagators (left panel) and the thermal widths (right panel) for X(3872), T + cc (3875), and their constituents in the pion gas as functions of the temperature T .In the left panel, the solid curves are the sums of the rest energy at T = 0 and the thermal mass shift for X (red lower curve) and T + cc (blue upper curve).The dashed curves and the horizontal lines are the corresponding charm-meson pair thresholds in the pion gas at temperature T and at T = 0, respectively.In the right panel, the solid curves are the sums of the decay width at T = 0 and the thermal width for X (red upper curve) and T + cc (blue lower curve).The horizontal lines are the decay widths of X (red upper line) and T + cc (blue lower line) at T = 0 MeV.The dashed curves are the sums of the widths of the constituents D * 0 and D 0 of X (red lower curve) and the constituents D * + and D 0 of T + cc (blue upper curve).
The real part of the pole of the pair propagator Re[E X (T )], the real part relative to the T -dependent D * 0 D 0 threshold δε 0 (T ) + δε * 0 (T ), the thermal width δΓ X (T ), and the thermal width relative to the sum δΓ 0 (T ) + δΓ * 0 (T ) of the thermal widths of the constituents for X(3872) at T = 0 and T = 100 MeV.The energies and widths are in MeV.
are both negative, approximately equal, and larger than our results by about one and two orders of magnitude, respectively.Our results in TableIRe[δε] Re[δε * ] δΓ δΓ *