Non-relativistic bound states at finite temperature (II): the muonic hydrogen

We illustrate how to apply modern effective field theory techniques and dimensional regularization to factorise the various scales which appear in QED bound states at finite temperature. We focus here on the muonic hydrogen atom. Vacuum polarization effects make the physics of this atom at finite temperature very close to that of heavy quarkonium states. We comment on the implications of our results for these states in the quark gluon plasma. In particular, we estimate the effects of a finite charm quark mass in the dissociation temperature of bottomonium.


I. INTRODUCTION
In a previous paper [1] we showed how to apply modern effective field theory (EFT) techniques to the hydrogen atom at finite temperature. They provide a systematic way to separate the physics occurring at the various dynamical scales involved in that system, which makes calculations simple and transparent. The main motivation of that work was to pave the way to a QCD based quantitative study of heavy quarkonium states in the quark gluon plasma (several works in this direction have recently appeared in the literature [2][3][4][5][6][7]), which share with the hydrogen atom a number of important features. The main qualitative difference, as far as the bound state dynamics is concerned, between heavy quarkonium states and the hydrogen atom is that vacuum polarization effects in the latter are very suppressed. For muonic hydrogen, however, the vacuum polarization effects provide the leading corrections to the Coulomb potential, as it is the case for heavy quarkonium states. This is our main motivation to study muonic hydrogen in detail here.
Muonic hydrogen is under current research at the Paul Scherrer Institute in the so called Muonic Hydrogen Lamb Shift experiment [8]. It allows for precision tests of QED which, among other things, probe the electromagnetic structure of the proton [9] or the size of the proton [10]. Theoretical calculations, on one hand, have achieved an impressive precision [11][12][13][14][15][16], and a number of experimental results are available [17,18]. However thermal effects on this atom due to black body radiation, or to electron-positron plasmas, have not been considered to our knowledge, either theoretically or experimentally. Current experimental facilities may now produce electron-positron plasmas [19], which are also the target of intense theoretical studies [20,21] (see [22] for a recent review). Making muonic hydrogen atoms slowly travel through an electron-positron plasma would be an ideal experiment to probe how well we understand thermal effects in non-relativistic bound states at relatively high temperatures. Recall that an analogous experiment at lower temperature with blackbody radiation on Rydberg atoms [23] first detected thermal level shifts, in the early eighties. This would mean a further step in taking advantage of the similarities of the electron-positron plasma with the quark-gluon plasma (see [24,25] for reviews) in order to learn about the non-trivial properties of the latter, as it has already been advocated by some authors [26].
In the center of mass frame, the proton of a muonic hydrogen is essentially at rest and the muon moves at small velocities v α ≪ 1. Hence, the relevant scales at zero temperature are those of a non-relativistic system [27]: the muon mass m µ (hard), the typical momentum p which is of the order m µ α/n (soft) and the energy E ∼ m µ α 2 /n 2 (ultrasoft), n being the principal quantum number. Unlike the hydrogen atom case, vacuum polarization effects introduce a new scale in the muonic hydrogen atom, the electron mass m e , which is of the order of the soft scale for the lower lying states (n = 1, 2), but larger for the remaining ones. At finite temperature further scales are introduced, not only T , the temperature, but also eT ∼ m D , the Debye mass, and others that will be discussed later. In order to efficiently deal with the physics at each of these scales we will use the effective theories of non-relativistic QED (NRQED) [27], suitable for energies much smaller than the hard scale, potential NRQED (pNRQED) [28], suitable for energies much smaller than the soft scale, and Hard Thermal Loop Effective Theory (HTL) [29], suitable for energies much smaller than the temperature, in a way analogous to Ref. [1]. Recall that pNRQED facilitates enormously the iteration of the Coulomb potential in ultrasoft contribution (at the scale E), and the HTL action does the same for soft thermal photon resummations at the scale eT . When the scales eT and E coincide the combination of pNRQED and HTL is crucial in order to obtain consistently both the iteration of the Coulomb potential and the resummation of soft thermal photons.
We will use the real-time formalism [30], which is mandatory for the study of the propagation of a (non-thermalised) non-relativistic system in a thermal bath. We shall restrict ourselves to temperatures much smaller than the muon mass, and hence the thermal bath does not affect the free muon propagator which remains the same as at zero temperature. The same holds true for the free proton propagator, which will be further approximated by that of a static source. Thermal propagators will in general be necessary for the photons, electrons and positrons. Recall that in the real-time formalism a doubling of degrees of freedom is required to properly account for the thermal propagation. External propagators can only correspond to type "1" fields, vertices contain either type "1" fields or type "2" fields. Propagators can be "11", "12", "21" or "22". When drawing Feynman diagrams we will understand that all possible types of vertices and propagators compatible with a given diagram are added up, and will not display each type explicitly (for muons and protons only the "11" propagator must be considered). The techniques and results we shall use have been reviewed in ref. [31]. We reproduced the basic ones in the Appendix B.
We distribute the paper as follows. In Section II we study the ideal case in which the electron mass is set to zero (m e = 0). This not only makes calculations simpler, but also makes the system closer to the heavy quarkonium case. In section III we focus on the actual case m e = 0. These two sections are divided in subsections in which the cases T ≪ p, T ∼ p and T ≫ p are studied, p being the typical relative momentum in the bound state (or the inverse Bohr radius). Section IV is devoted to discussion and conclusions.
II. me = 0 CASE Let us first consider an ideal case in which the electron mass m e is taken to be zero. This case is in fact closer to the one in heavy quarkonium states, particularly in charmonium, than the actual case with m e = 0, which we will study in the next section. It has already been discussed in the past in order to clarify subtle issues on the renormalization group structure of non-relativistic effective theories [32].
For temperatures much smaller than the soft scale (m µ α in this case) we can study the thermal effects in the atom starting from the pNRQED Lagrangian at zero temperature, up to exponentially suppressed contributions ∼ e −p/T . This means that the potentials are the same as the ones at zero temperature. The only difference with respect to the hydrogen atom case is that these potentials contain now O(α) corrections due to vacuum polarization effects produced by virtual electron-positron pairs. Like in the hydrogen atom, the ultrasoft photons, electrons and positrons are responsible for the finite temperature effects. Our starting point in this section is then eq. (6) from [33].
where S(t, r, R) is the muon wave-function field, r being its distance to the proton and R the position of the proton; e(t, x) the electron Dirac field. α = e 2 /4π is the electromagnetic coupling constant, and c 0 , c s and d 2 are matchings coefficients which can be found at one loop in [34]. Let us separate the cases T E and T ≫ E, which are analysed the two following sections: In this case, the leading temperature-dependent contributions are given by the diagram in fig. 1, in an analogous way to the hydrogen atom case. Virtual ultrasoft electron-positron pairs give rise to O(α) corrections and no soft thermal photon resummation is necessary at the scale E. Hence, there is no qualitative difference with respect to the hydrogen atom, and we will not further discuss it. We refer to [1] for the relevant formulas for the spectrum and decay widths 1 .

T ≫ E
In this case the scale T can be integrated out before calculating the spectrum and decay widths, we call the resulting effective theory pN RQED >T (for further explanations about the notation see Appendix A). In the photon and electron-positron sector this gives rise to the HTL action [29]. In the atom sector, the pNRQED Lagrangian gets additional temperature-dependent potentials. At leading order (LO) in α , they arise from the diagram in fig. 1 upon expanding E − H in the integrals, and have been calculated in [1] δV H ∼ E ∼ α/r. Possible O(αr 4 ET 4 ) terms arising from higher orders in the multipole expansion cancel out. Note that the dominant contribution above is a constant mass shift, and the r-dependent part is (m µ α 2 /n 2 T ) 2 suppressed. Hence, vacuum polarization corrections to the photon propagator may compete with the r-dependent part of the LO potential displayed above and must be calculated. From the diagrams in fig. 2 we obtain the next-to leading order (NLO) in α, with m 2 D = (eT ) 2 /3 and c a numerical constant. The computation have been done in dimensional regularization with ǫ = (4 − d)/2 → 0. The first line of this result also appears in an analogous calculation that has already been carried out in the static limit of the QCD case [6]. The second line is subleading. We have displayed it to match the precision of (2) when eT ∼ E. In order to eventually check the correct cancellation of the 1/ǫ poles we have calculated analytically the leading IR behavior in the appendix C. The constant c remains unknown. We see that indeed δV , except for the global energy shift given by the first term in (2). In fact, it provides the dominant term in the potential for the energy of the photon transitions between two states belonging to this case. The IR divergences arising above are canceled by UV divergences divergences arising from contributions at smaller scales (E, eT , . . .). These contributions are hard to calculate in the general case because HTL propagators must be used for the ultrasoft photons and the Coulomb potential must be kept unexpanded in the atom propagator. At these scales, however, the Bose distribution can be expanded, which simplifies somewhat the calculations, and produces the so called Bose enhancement, see (B8) and (B10) below. The dominant contribution arises from the diagrams of fig. 3. We have only been able to work out an analytic expressions for the cases eT ≪ E and eT ≫ E, which will be discussed below, and for its UV behavior. Technical details for the latter are shown in the Appendix C.
Before discussing the two cases which allow to proceed further analytically, we display the energy shift and decay width induced by the temperature-dependent potential (2) and (3) at leading order in perturbation theory, in which r 2 n = n 2 2m 2 µ α 2 [5n 2 + 1 − 3l(l + 1)], n and l being here the principal and angular momentum quantum numbers. The labels n, m, . . . are also used through the paper as a short hand notation for the whole ensemble of quantum numbers of a given Coulomb state, either bound or in the continuum, and φ n (0) is the wave function at the origin. The contributions above are to be added to the ones coming from lower scales, which we discuss below in two particular cases that share the feature that the lower energy scales are hierarchically ordered, and hence the method of the integration-by-regions can be used [36,37] • E ≫ eT In this case the loop integral is dominated by energy and momentum ∼ E for which eT can be treated as a perturbation. At leading order in the HTL expansion we obtain We observe that the leading order infrared divergences appearing at the scale T in (4) and (5) are canceled by the ultraviolet divergences in (6) and (7) respectively. Note that the subleading infrared divergence at the scale T , in the second line of (5), is very much suppressed in this case (m µ α 5 ≫ αm 2 D /m µ ). A similar calculation for the heavy quarkonium case has been presented in [38]. Technical details can be found in the appendix (C 2). We only mention here that a collinear region exists that contributes at this order.
Upon summing up the contributions from both energy regions, namely (4) and (6) for the energy and (5) and (7) for the decay width, finite results are obtained at the desired order, • eT ≫ E, In this case the loop integral is dominated by energy and momentum ∼ eT for which E can be treated as a perturbation. At leading order in the energy expansion we obtain a contribution which is equivalent to adding a new term to the potential that goes like m 2 D r 2 .
We observe that the leading order infrared divergence appearing at the scale T is cancelled by the ultraviolet divergence above. An analogous contribution has also been calculated in the static limit of QCD [6].
At next-to-leading order in the energy expansion we have restricted ourselves to compute the ultraviolet divergence analytically, It cancels the infrared divergence in the second line of (3), as it should (c * is an unknown constant that can be of order 1/α 1/2 because of Bose-enhancement).
Fortunately the contribution of the loop integral for energy and momenta ∼ E is subleading. The calculation at that scale may even require non-perturbative techniques if E gets close to the scale e 2 T [30, 39,40].
Summing up the contributions from the T energy region and from the m D energy region, namely (4) and (10) for the energy and (5) and (11) for the decay width, the leading thermal effects for this situation are obtained, Since m µ ≫ T still holds, our starting point is the NRQED Lagrangian at T = 0 [27] (this is correct up to exponentially small contributions ∼ e −mµ/T ). The thick and the extra-thick solid lines stand for the non-relativistic propagator of the muon and the static propagator of the proton respectively. The remaining lines are as in fig. 2. There is also a leading order correction coming from an extra diagram, which is obtained from the second one by changing the muon line by a proton line.
where ψ and N are the muon and proton Pauli-spinor fields respectively.
Since T ∼ p ≫ E, we can integrate out the scales T and p, which leads to what we call pN RQED T , the suitable effective theory for the scale E, which is similar to the one that has already been introduced in section II.A, with the only difference that now T is of the same order as the cut off of the effective theory. In the photon and electronpositron sectors we have the standard HTL. In the atom sector temperature-dependent potentials are induced. Recall that at the scale T there is no enhancement and vacuum polarization effects due to electron-positron pairs are always suppressed by α. Hence, the leading potential will still be the Coulomb potential but the first α correction to it will already be temperature-dependent. This is given by the diagram in fig. 4.
The temperature-dependent part of the potential induced by the diagram in fig. 4 is both UV and IR safe in momentum space. However, when it is Fourier transformed to coordinate space an IR divergence is encountered. The calculation is rather involved so we only display here the final result, which can be given in terms of one-parameter integrals of special functions. Details are given in the Appendix D (the formulas in that appendix have to be used setting m e = 0 for this case). We obtain, where ρ = 2rT and Si stands for Sine Integral The LO energy correction is obtained by computing the expectation value of the potential of (16) for the desired state and adding the ultrasoft contribution. The calculation in pN RQED T is identical to one carried out in the second case of sect.
Hence the outcome can be directly read off equations (10) and (11). Notice that the infrared divergence in (16) at first order in quantum mechanic perturbation theory, induces a contribution that cancels out the ultraviolet divergence in (11), as it should. Also, the contribution from integrating out the scale m D can be encoded in a correction to the potential, summing it up to δV r the following finite result is obtained, − and hence, at first order in perturbation theory, Since m µ ≫ T still holds, we can also start from NRQED at T = 0. Now we may proceed by sequentially integrating out first the scale T and next the scale p. After integration of the scale T we get an effective theory which consist of HTL contributions in the photon and electron-positron sector, and of NRQED with temperature-dependent matching coefficients in the atom sector. This N RQED >T in the atom sector is identical to the one that we have in the hydrogen atom case [1] , up to order α corrections induced by the electron-positron vacuum polarization.
The next step is to integrate out the energy scale p, namely to match N RQED >T to what will be called pN RQED <T , which is expected to produce temperature-dependent potentials. These potentials must be calculated using HTL photon propagators. Let us separate the two following cases: In this case the computations can be carried out as in section V.B of [1]. The relevant diagram is similar to fig. 4, but instead of a photon propagator with a self-energy insertion we would have to use the tree level HTL photon propagator ( fig. 5). The only difference with respect to [1] is due to the fact that the electron-positron pairs that generated the HTL photon propagators are now taken to be massless, so, in fact, the outcome is simpler: the non-trivial function g(m e β) reduces to πm 2 D 16αT 2 . We obtain then the following leading order potential where This potential coincides with the one first obtained in [2] for QCD (up to trivial changes in color factors made explicit in [1]). As in the hydrogen atom case, we can use this result in order to estimate the dissociation temperature, which is T d ∼ m µ α 2/3 / ln 1/3 α, as anticipated in [1] for the QCD case.
• eT ≫ p, In this case T ≫ T d always holds. Hence, the imaginary part of the potential is bigger than the real one, so it does not make much sense to speak about bound states anymore.
We address now the actual case of a non-vanishing electron mass. Although the real muonic hydrogen is not as close to heavy quarkonium systems as the ideal one (m e = 0), it may still be useful to learn about certain aspects of it. In particular about the rôle of a finite charm mass in the bottomonium system, which is analogous to that of a finite electron mass in muonic hydrogen [41]. This case may then shed light on the effects of the charm quark mass in bottomonium at finite temperature, specially when charm quarks are thermalized. Irrespective of that, muonic hydrogen is a real system that appears in nature, which nowadays is produced in large samples [8], and, therefore, our results may eventually be checked against experiment.
For actual muonic hydrogen m e ∼ p for the lower lying states (n = 1, 2), whereas for the remaining states (n ≥ 3) one may safely consider m e ≫ p [41]. Let us then analyse these two cases separately.
A. Lower lying states (n = 1, 2) As mention above, these two states fulfill p ∼ m e , and hence relativistic electron-positron pairs must be integrated at the same time as the momentum transfer p is. Let us see in the following sections how this is carried out depending on what the temperature is.
Like in the massless case our starting point can be pNRQED. However, now, due to the fact that m e ∼ p (rather than m e = 0), the electron-positron pairs have already been integrated out when calculating the potentials, and hence are not active anymore. The situation is then totally analogous to the hydrogen atom, the only difference being that the potentials get O(α) corrections due to virtual electron-positron pairs, the most important of which is the Uehling potential. In other words, the thermal bath contains neither electrons nor positrons, so the thermal effects are only due to the photons which do not distinguish between electrons and muons. Hence the results concerning this case can be read off section III of ref. [1] by making m → m µ [up to O(α) corrections].

T ∼ p
Again, like in the massless case, we can start with NRQED. The scales T and m e must integrated out at the same time as the energy scale p. In the photon sector, which is not sensitive to the scale p, we get the mass dependent HTL action (see section V.A.1 of [1]). In the electron and positron sectors, which are not sensitive to the scale p either, we get a NRQED T Lagrangian for each of these particles (see section V.A.2 of [1]). In the atom sector, the potentials depend now on both temperature and the electron mass, except for the leading Coulomb potential. The most important correction is a kind of temperature-dependent Uehling potential, which is obtained from the diagram in fig. 4, Further expressions for these functions can be found in the Appendix B of [1]. The computations that lead to (22) are carried out in Appendix D.
Notice that (22) has two infrared divergences, which, as in the m e = 0 case, arise when the Fourier transform of the momentum space potential is taken, in order to get the coordinate space potential. The IR divergence in the fourth line of (22) is similar to the one that appears in equation (16) for the massless case (with T 2 g(m e β) instead of m 2 D ). The IR divergence in the second last line of (22), however, is proportional to m 2 e e βme +1 , and hence distinct of the m e = 0 case. It emerges from a region in which not only the three-momentum transfer is small but also the component of the three-momentum of the electron-positron pair in the loop parallel to the momentum transfer is small. In either case, these IR divergences should cancel against UV divergences in the pNRQED calculation.
The relevant diagram in the pNRQED calculation is again fig. 3, in which the photon line must be understood as the mass-dependent HTL propagator (see Appendix E). In the dominant contribution to this diagram, E − H in the atom propagator can be treated as a perturbation (recall that E − H ∼ m µ α 2 ≪ eT ). Then using the fact that ∆ 11 (k 0 , k) is symmetric with respect to k 0 → −k 0 , we obtain We see that indeed the UV divergences above cancel those of (22), as expected. There is a subtle point in this calculation, however, which we discuss in the Appendix E, that must be correctly dealt with in order to get the UV divergence of the last line (that cancels the IR divergence in the second last line of (22)). As in the massless case, the contribution from the scale m D can be encoded in a correction to the potential, and this summed to δV r , 3πT (e βme +1) Γ(−2/3) + O(α 3 T ). The energy shift and the decay width at first order in perturbation theory can be obtained from (19).

T ≫ p
This case is very similar to the m e = 0 one. We start with NRQED and integrate out the scale T first. Since T ≫ m e , the mass-dependent HTL propagators may be expanded in m e /T , and hence become the usual HTL propagators with next-to-leading order contributions due to the non-vanishing electron mass . Hence the finite mass effects do not affect the gross features of the system. In particular, in the eT ∼ p case the dissociation temperature will be similar to the one in the massless case. For eT ≫ p, like in the massless case, no bound state is expected to survive.

B. Higher energy states (n ≥ 3)
As mentioned before, these states fulfill m e ≫ p ≫ E, and hence relativistic electron-positron pairs may be integrated out before the momentum transfer p is. Let us see in the following sections how this is carried out depending on what the temperature is.

me ≫ T
In this case we can start with a NRQED Lagrangian for the muon in which the relativistic electron-positron pairs have already been integrated out, which gives rise to 1/m 2 e corrections to the Maxwell Lagrangian. Since there are neither electrons nor positrons in the thermal bath, the different situations coincide with those of the hydrogen atom, and the relevant expressions can be read off sections III and IV of [1] by making m → m µ (up to 1/m 2 e corrections).

me T
In this case we can start with a NRQED Lagrangian for the muon keeping the relativistic Dirac Lagrangian for the electron in it. When m e ∼ T both scales must be integrated out at the same time. In the photon sector, a mass-dependent HTL Lagrangian is induced (see section V.A.1 of [1]). In the electron and positron sectors a temperature-dependent NRQED T Lagrangian for each particle is induced (see section V.A.1 of [1]). In the muon sector, a temperature-dependent NRQED >T Lagrangian is also induced. At lower orders, it can be obtained from the diagrams of section IV of [1], together with diagrams containing an electron-positron loop. The most important effect is the appearance of two mass shifts, one ∼ αT 2 /m µ from the diagram (38) of [1] and the other one ∼ α 2 m e from the second diagram in fig. 4. In the proton sector an analogous mass shift ∼ α 2 m e occurs, that is due to a diagram obtained from the previous one by changing the muon line by a proton line.
The next step is to integrate out the scale p, namely matching to pNRQED using the HTL Lagrangian. This has already been done in section V.B of [1]. At leading order, this produces a temperature-dependent potential and further mass shifts, which can be read of a corrected version of formulas (58) and (60) in that reference (making m → m e ). The origin of the corrections is discussed in Appendix E, and boils down to a simple replacement, see (E17). Both the potential and the mass shift contain an imaginary part, which becomes more important than the real part starting at some temperature T d , which we call dissociation temperature. In the following section the dissociation temperatures will be estimate for several states.
For T ≫ m e none of the higher energy states exists anymore, since the dissociation temperatures fulfill m e > T , see next section.

C. Dissociation temperatures, level shifts and decay widths
The dissociation temperatures will be estimated in a similar way as they were in the hydrogen atom case [1]. This is by identifying the momentum scale for which the real and imaginary part of the momentum space potential are equal, p ∼ (16α) 1/3 (g(m e β) + (m e β) 2 n F (m e β)/2) 1/3 T =: m d , and equating it to the typical momentum transfer in the muonic hydrogen atom 2 , p ∼ m µ α/n 2 . The results are displayed in table I. Table II shows the same results for the hydrogen atom. In order to carry out the estimates we have use the formulas of Appendix E, which are meant This result was computed using the assumption that 1 r ≫ mD, so we expect important deviations from the real energy starting at T ∼ 4M eV for higher excited states (T ∼ m e ≫ p), for the case of the lower lying states, instead of the m e = 0 formulas. This is indeed legitimated: We are just not taking advantage of the fact that T d ≫ m e for these states, which produces the simplifications discussed in section III A 3.   From the results in table I, we see that only the lowest lying states n = 1, 2 survive at temperatures of the electronpositron plasma m e T Hence, only transitions between these levels might be observed in an eventual experiment. We shall focus on the experimentally prominent K α transition [17]. We display our results for the energies of the 1S and 2P states, and for the energy of the K α transition as a function of temperature in the range T ∈ (2, 0.05)M eV. in fig. 6, fig. 7 and fig. 8 respectively. The calculations have been carried out numerically, using first order perturbation theory for the potential (27).
Another observable that we can predict with our results is the decay width. In fig. 9 we show the decay width for the 1S state as a function of the temperature. Comparing fig. 6 and fig. 9 it can be seen that T d ∼ 1.7M eV is the temperature that makes the decay width of the same magnitude as the binding energy. In this paper we have discussed the properties of muonic hydrogen in a thermal bath, which may consist not only of blackbody radiation but also of an electron-positron plasma. We have further developed the effective theory techniques for bound state systems at finite temperature initiated in [1], in particular the application of dimensional regularization to the factorization of the various scales in the system. They facilitate enormously the organisation of the calculation. For instance, they make apparent when Coulomb or HTL resummations are necessary and when they are not. In addition, both partial and final results are naturally obtained as a series of small scales over large ones, thus providing a good control on the systematics.
We have discussed two cases. We have first addressed the academic case of muonic hydrogen with a vanishing electron mass, which turns out to be closer to heavy quarkonium states than the actual case with a non-vanishing electron mass, that we have addressed next. All the thermal modifications we have found turn out to be spin independent.
In the zero electron mass case, we have studied how the effects of vacuum polarization modify the picture that we encounter in normal hydrogen [1]. The modifications turn out to be important when the temperature is larger than the binding energy. For instance, they would give the leading order contribution to a hypothetical K α transition for high enough temperatures (2). For temperatures below dissociation, we have presented the leading order, and selected next-to-leading order, thermal corrections to the binding energy and decay width.
In the actual electron mass case, muonic hydrogen behaves very much the same as hydrogen for temperatures below the electron mass. For temperatures larger or of the order of the electron mass the vacuum polarization effects are sizable, and, at some point, make the bound states dissociate. We display in table I the dissociation temperature for the lower laying states. We have also calculated the thermal modifications to a number of observables before dissociation occurred. For instance, we plot the dependence of the K α transition on temperature in fig. 8, which could be tested experimentally in the future [8].
We close with a concrete application to the heavy quarkonium case. As we have mentioned before, the way a finite electron mass affects muonic hydrogen is similar to the way a finite charm quark mass affects bottomonium [41]. Since this should also be the case at finite temperature, we can easily translate to the QCD case the results for the dissociation temperature of muonic hydrogen, which we show in table III.  all the effective field theories that may arise from integrating out the different degrees of freedom. In this paper we have used the following notation. Basically we name the effective field theories as one would do for zero temperature, and we encode the temperature information in a subindex. The subindex T means that the temperature has been integrated out, the subindex m D means that also the scale eT has been integrated out. Since the matching coefficients of the effective field theory that we obtain after integrating out m µ (p) and T is not the same if m µ (p) ∼ T or if m µ (p) >> T we include a symbol <, > or blank depending of the relation between these scales. For example, if we are in m µ (p) ∼ T we will arrive to N RQED T (pN RQED T ), but if m µ (p) >> T we reach N RQED >T (pN RQED >T ) because T is smaller than the energy cutoff of N RQED (pN RQED).

Appendix B: Basic formulas
In this appendix we display a number of formulas of the real-time formalism that are relevant for the paper. Our notation closely follows ref. [31]. Recall that in this formalism there is a doubling of degrees of freedom [30]. Fields are labeled as "1" or "2". Fields "1" ("2") only interact with fields "1" ("2") according to (minus) the original interaction Lagrangian. Fields "1" may be converted to fields "2", and vice versa, through propagation so that propagators become 2 × 2 matrices. For instance, for a free scalar field we have where i∆(K) is the Feynman propagator, and n B (k 0 ) the Bose distribution function.
For the tree-level propagator of the transverse electromagnetic field in the Coulomb gauge, i∆ ij (K), we have The tree-level A 0 -propagator matrix in the Coulomb gauge is diagonal, traceless and the "11" component coincides with the propagator at zero temperature. For the tree-level propagator of a Dirac fermion field, iS(K), we have S(K) ≡ (K / + m)∆(K), where∆(K) follows from ∆(K) by replacing n B by −n F , n F being the Fermi-Dirac distribution.
Since we are always in the case m µ >> T all vertices involving muons or protons will be of type "1". However, vertices involving photons and electrons can be both type "1" and "2". At the order we are calculating, it turns out that we only need propagators of the type "11" for the photons, either at tree level or including one-loop self-energies.
For computations that require loop corrected propagators (for example fig. 3 and fig. 2) it is convenient to use the so called Keldysh representation [31]. In this representation the retarded, advanced and symmetric propagators are defined as, Notice from (B1) that at tree level only ∆ S depends on the temperature. In order to calculate loop corrections to the ∆ 11 propagator in an efficient way, we use the following method [31].
• We compute the ∆ R propagator, using the fact that for this propagator the Dyson equation is of the zero temperature type (this is not so for ∆ S ).
where the self-energy Π R is and ∆ 0 R is the tree-level retarded propagator obtained according to formulas (B1)-(B3). • The advanced propagator ∆ A is the complex conjugate of the retarded one and the symmetric propagator ∆ S reads, in the bosonic case, • Finally, one obtains Note that for |k 0 | ≪ T the size of the symmetric propagator for bosons is larger than what one would expect from naive power counting plus terms suppressed by 1/T . This effect is called Bose-enhancement, and it complicates the power counting of the EFTs at scales lower than the temperature.
Appendix C: Calculations in Section II.A All these computations have in common that the starting point is pNRQED. In this effective theory, and in all the theories that are derived from it by integrating out further degrees of freedom, the leading correction to the Hamiltonian is [6,28] with ∆ ij and ∆ 00 being the "11" transverse and A 0 photon propagators respectively. δH is a potential δV (energyindependent) only if E is much smaller than the scales inside the ∆ propagators. Depending on the concrete calculation that one is doing the explicit expression for ∆ may be different.

Integrating out the T scale
After integrating out the scale T we will reach pN RQED >T . At leading order the difference between pN RQED and pN RQED >T will be a correction on the potential of the type of (C1) where the internal momentum K is of order T .
The LO correction (2) is obtained using the tree level photon propagator in (C1). If one uses the Coulomb gauge thermal effects only appear in transverse photons. This computation is done in detail for different relations between T and E in [1].
The NLO potential comes from including one-loop corrections to the photon propagator (this corresponds to the diagram of fig. 3). This can be done analytically for T >> E. The required expression for the photon propagator are found in section III.B of [6]. Since T >> E the muonic hydrogen propagator can be expanded Let us compute the first term in this expansion Looking at the expressions of [6] and using the relation (B9) one sees that both ∆ ij and ∆ 00 are even functions in k 0 , so one can use for the first term of the expansion So we only need the photon propagator in the limit k 0 → 0 to get the leading order.
In [6] the propagators is given as an integral of a parameter q 0 . The best strategy to perform the calculation is to integrate first over internal momentum and leave the integration of this parameter to the end. This computation was done in section IV.B.1 of [6], here we take the Abelian limit making C F = 1 and C A = 0.
For the next-to-leading order in (3), namely the contributions coming from second term in (C3), we will restrict ourself to the extraction of the leading IR logarithmic behavior. We use the fact that in the infrared the photon self-energy approaches the HTL limit, so instead, we substitute the complete self-energy Π HT L e −β 2 k 2 . The factor e −β 2 k 2 is introduced to regulate UV divergences and does not affect the leading IR behavior we are interested in. The HTL self-energy have the property that can be written as Π HT L = m 2 D f (k 0 /k) with f (x) a non trivial function. Performing the change of variables k 0 = kx the infrared behavior can be easily extracted. First we compute the contribution from the retarded part of the longitudinal photon propagator.
We have used that r i (E − H)r i = − 1 2mµ plus terms that vanish on the physical state. As we are only interested in the leading logarithm behavior, the change Π The equality above is only true as far as the leading IR behavior is concerned. Now with the change k 0 = kx (C10) Ω D−1 is the D − 1 dimensional solid angle. Note that the contribution from the advanced part of the longitudinal photon propagator is obtained by replacing f L R (x) → (f L R ) * (x) in (C10). Since it has all the singularities in the upper half complex plain the corresponding integral gives zero. The integral over x in (C10) can be done using standard techniques from integration in the complex plane, and the result is The contribution from the symmetric part of the longitudinal photon propagator does not produce IR divergences in dimensional regularization. So we do not need to compute it here as we are only interested in the logarithms. Next we proceed analogously for the transverse photon propagator. As in the longitudinal part, only the retarded plus advanced contribution is infrared sensible. We approximate Π T then (C13) The first (second) term in the square brackets have all the singularities in lower (upper) complex k 0 half plane, and hence the whole expression vanishes. The contribution from the symmetric piece of the propagator also vanishes for the same reason. Hence, In order to obtain (6) and (7) the starting point is pN RQED T . In this effective theory the atom self-energy gives the correction to the Hamiltonian (C1). Recall that the photon propagators in (C1) must be taken in the HTL approximation. We use the method of the integration by regions [36,37] in order to evaluate it. In this integral there are three relevant regions: k, k 0 ∼ E with λ = k 0 − k ∼ E (we call it off-shell region), k, k 0 ∼ E with λ ∼ m D (collinear region), and k, k 0 ∼ m D (m D -region). By power counting one can see that the m D -region will be just a higher order correction, so we focus in the regions k, k 0 ∼ E. In the off-shell region the atom propagator can not be expanded, but the HTL photon propagator can. For example, the longitudinal photon retarded propagator can be written as and, after the compulsory expansion The fact that the non-trivial functions appears only in the numerator after the expansion is crucial in order to be able to do the integration analytically. The collinear region does not contribute in the part that is related with longitudinal photons.
The transverse photon retarded propagator is In the off-shell region a expansion like the one for the longitudinal photon propagator has to be made. In the collinear region (that has k, k 0 ∼ E but k 2 0 − k 2 ∼ m 2 D ), the atom propagator can not be expanded also, but the HTL photon propagator has to be expanded around the region |k 0 /k| ∼ 1. As in the previous case, this expansion makes it possible to proceed analytically.
Details of this computation can be found in [38], 3. Integrating out mD for mD >> E After integrating out m D we will arrive to what we will call pN RQED >mD . In the photon sector we will have a non trivial action but we will not need it at the level of precision we are working. In the atom sector there appears a correction of the potential in the matching between pN RQED T (or pN RQED >T ) and pN RQED >mD of the form of (C1) where the internal momentum is of order k ∼ m D .
Because m D >> E we can also put E − H = 0 at leading order in the atom propagator as in section C 1. Hence, we will need the HTL photon propagator in the k 0 → 0 limit(i.e. T >> k >> k 0 ). The HTL photon propagators are very well known and can be found in [6,30,31]. In our case the fact that k 0 → 0 makes Using this in (C1) the results (10) and (11) are obtained. For the next-to-leading order (12) the calculation is very similar to the one we have carried out for the T scale. Like in that case, we will focus on the logarithmic behavior, now in the UV. As an example, we study again the retarded part of the longitudinal photon propagator (note that consistency with eq. (C14) ensures that there will not be a logarithmic contribution from transverse photons) using that scale-less integrals in dimensional regularization are zero Now is when the change k 0 = kx becomes useful If one is only interested in the logarithmic behavior, one can put D − 3 = 1 in the x integration and use (C11). Formula (12) is readily obtained from the expression above.

UV behavior for E ∼ mD
Although in this situation we have not been able to obtain an analytic result for (C1), its UV behavior can be easily isolated as follows.
The piece in brackets in the right hand side of the equation above leads to a ultraviolet finite expression when substituted in (C1). So the ultraviolet divergences arise from the remaining terms in the right hand side of this equation. In fact the computation for these UV divergences is exactly the same as in the case m D >> E, which we have carried in section C 3.
Appendix D: Calculations in Section III.A

Correction to the Coulomb potential in pN RQEDT
In this part we deal with the matching procedure that has to be done for T ∼ p. In perturbation theory the pN RQED potential is related to the Fourier transform of the longitudinal photon propagator in the limit where p 0 → 0 (for P the external momentum of the propagator). For this temperature range the propagator, which can be obtained by the procedure outlined in Appendix B, is needed for T ∼ m e ∼ p >> p 0 . The retarded self-energy reads .
For m e → 0, this self-energy coincides with the Abelian limit of the one found in [6]. In order to obtain the potential the first step is to use formulas (B9) and (B8) to get the corrections to the "11" propagator. Then the propagator is related with to potential plus the self-energy by the following formula at leading order this gives the Coulomb potential V (r) = − α r . For simplicity we define and For the next to leading order we need the corrections to the propagator The contribution of (D7) to (D2) leads to the first term in the first line of (22). In order to calculate the contribution (D8) to (D2), it is convenient to leave the integration over the internal momentum k to the end. Consider then, and use Some terms vanish because of the symmetry λ → −λ. The integral over p in the first and the second term is straight-forward. By using the symmetries in the λ and p variables, the third term can be simplified as follows, . (D15) At this point the integral over p can be done using standard techniques of complex analysis,  4 (D17) After performing the integration in λ we obtain the second and third line in (22). ∆ 3 can be computed in a very similar way, and leads to the second term of the first line in (22). Let us next consider ∆ 4 , .

(D18)
We have separated above the pieces that lead to infrared divergences from the ones that do not. The first term inside the square brackets gives the fourth line in (22), and the second term together with V m gives the fifth line of (22). The rest of terms in the square brackets give the sixth and seventh line of (22). ∆ 5 can be computed in a very similar way, and gives the remaining terms of (22).
Appendix E: Computation of the HTL retarded self-energy for the longitudinal photon in the me = 0 case There are some subtle points in the computation that leads to (26), which do not arise in the massless case and are worth elaborating upon. Let us start from the formula (37) of [31] for the massive case, It is customary to make the change K → −Q in the second term to get a simplified expression that reduces to twice the first term. However, the terms proportional to m 2 e , which do not exist in the massless case, have a stronger IR sensibility than the remaining ones. This leads to ill-defined expressions in the HTL approximation for T ∼ m e . These expressions must be properly defined in order to get consistent results before and after the shift K → −Q has been carried out.
To see this in detail we work out just a small part of the computation which illustrates the point, namely the part of ℑΠ L R (P ) that comes only from the m 2 e term in the numerator. We call this term Π m 2 e (P ), By carrying out the HTL expansion and taking into account that we must expand also p 0 ( p ≫ p 0 in the computation of the potential), we obtain, For illustration purposes, we will make this computation in two ways • Using d D−1 k = dΩ D−2 d(cos θ)k D−2 dk, where θ is the angle between the internal momentum and the external one.
This expression has an end-point singularity, but one can skip it in dimensional regularization by choosing D − 4 > 0. So the result is zero in this way.
• Using d D−1 k = d D−2 k ⊥ dk z We choose z to be the direction parallel to the external momenta. 3 (E8) Notice that the integrand vanishes for all k z except when k z ∼ 0, so one may substitute This simplifies the computation a lot because complex plane integration techniques can be applied, This expression does not have any end-point singularity, and in fact it does not need a regularization anymore, ℑΠ m 2 e = − e 2 m 2 e p 0 2πp(e βme + 1) (E11) This is indeed the expected result, that eventually leads to terms contributing to the second line of (26). We arrive then at the paradoxical situation in which the final result depends on the precise way DR is implemented, either in spherical coordinates or in Cartesian ones. The apparent contradiction is resolved by noticing that DR in spherical coordinates does not allow for the shift K → −Q. In order to show this is the actual reason for it, let us start now from eq. (E2) So far we have used the complete expression for the self-energy. We apply next HTL expansion and also p ≫ p 0 , like above, which agrees with (E11). If the expansion for small p 0 is not carried out the same result is found with the three methods, because then no regularization is needed. The conclusion is that when the formula for the retarded self-energy in the massive case [1], i.e.
is expanded for p ≫ p 0 , DR must be used in Cartesian coordinates in order to properly regulate IR divergences (i.e to be consistent with the shift of momenta carried out at some point in order to get the expression above). This point was overlooked when calculating the potential in formula (57) of [1], and terms analogous to (E14) were missed. The correct formula is obtained by making the following substitution in (57) of [1], g(mβ) → g(mβ) + m 2 β 2 2(e βm + 1) (E17)