Saha ionization equation in early universe

TL;DR

This paper examines the ionization equilibrium of hydrogen in the early universe using thermal field theory, comparing ionization-recombination rates with cosmic expansion to understand matter-radiation interactions.

Contribution

It introduces a real-time thermal field theory approach to calculate hydrogen ionization and recombination rates in the early universe, focusing on ground state atoms.

Findings

Calculated ionization-recombination rate differences

Compared rates with universe expansion at various temperatures

Provided insights into hydrogen atom equilibration in the early universe

Abstract

The Saha equation follows from thermal equilibrium of matter and radiation. We discuss this problem of equilibrium in the early universe, when matter consists mostly of electrons, protons and hydrogen atoms. Taking H-atoms in their ground state only and applying the real time formulation of thermal field theory, we calculate the difference of ionization and recombination rates, which controls the equilibration of H-atoms. This rate is compared with the expansion rate of the universe at different temperatures.

Tables

Table 1. Table 1: Reaction rate ( Γ Γ \Gamma ) and expansion rate ( H 𝐻 H ) at different temperatures

$T$(in K)	$\Gamma$ in ${\text{s}}^{-1}$	$H$ in ${\text{s}}^{-1}$
6000	2.6$\times {10}^{-11}$	1.6$\times {10}^{-13}$
5000	1.0$\times {10}^{-13}$	1.2$\times {10}^{-13}$
4000	2.7$\times {10}^{-17}$	8.4$\times {10}^{-14}$
3000	3.4$\times {10}^{-23}$	5.2$\times {10}^{-14}$

Equations

$$\gamma +H\rightleftharpoons p^{+}+e^{-},$$

$${\mathcal{L}}_{\text{mbox}int}=g(\Phi {\overline{\psi }}_p{\gamma }^{\mu }{\psi }_e+h.c.)A_{\mu }.$$

$$g=\frac{e}{m},$$

$$\mathrm{Im}{\overline{\Sigma }}_{(1)}=16\pi g^{2}mM\int \frac{d^3{k}_1}{(2\pi {)}^{3}2{\omega }_1}\frac{d^3{k}_2}{(2\pi {)}^{3}2{\omega }_2}\frac{1}{2{\omega }_3}[n_1(1-{\tilde{n}}_2)(1-{\tilde{n}}_3)-(1+n_1){\tilde{n}}_2{\tilde{n}}_3]\delta (\omega +{\omega }_1-{\omega }_2-{\omega }_3),$$

$$\frac{\mathrm{Im}{\overline{\Sigma }}_{(1)}}{\omega }={\Gamma }_d-{\Gamma }_i,$$

$$n_1=\frac{\exp \left(-\beta {\omega }_1/2\right)}{\exp (\beta {\omega }_1/2)-\exp (-\beta {\omega }_1/2)},1+n_1=\frac{\exp \left(\beta {\omega }_1/2\right)}{\exp (\beta {\omega }_1/2)-\exp (-\beta {\omega }_1/2)}$$

$${\Gamma }_d;{\Gamma }_i=16\pi g^{2}m\frac{M}{\omega }\{\exp [\beta (\omega -{\mu }_2-{\mu }_3)/2];\exp [-\beta (\omega -{\mu }_2-{\mu }_3)/2]\}\times L,$$

$$L=\int \frac{d^3{k}_1}{(2\pi {)}^{3}2{\omega }_1}\frac{d^3{k}_2}{(2\pi {)}^{3}2{\omega }_2}\frac{1}{2{\omega }_3}\frac{\delta (\omega +{\omega }_1-{\omega }_2-{\omega }_3)}{\prod \left\{\exp [\beta ({\omega }_i-{\mu }_i)/2]\mp \exp [-\beta ({\omega }_i-{\mu }_i)/2]\right\}},$$

$$\frac{{\Gamma }_d}{{\Gamma }_i}=\exp [\beta (\omega -{\mu }_2-{\mu }_3)].$$

$$\frac{dn(\omega ,t)}{dt}=(1+n){\Gamma }_i-n{\Gamma }_d,$$

$$n(\omega ,t)$$

$$=\frac{1}{\exp \beta (\omega -{\mu }_2-{\mu }_3)-1}+c(\omega )e^{-\Gamma t},\Gamma ={\Gamma }_d-{\Gamma }_i,$$

$$\mu ={\mu }_2+{\mu }_3,$$

$$\omega =M_H+\frac{q^2}{2M_H},{\omega }_1=k_1,{\omega }_2=m+\frac{k_2^2}{2m},{\omega }_3=M+\frac{k_3^2}{2M}$$

$${\tilde{n}}_2=\exp \left[-\beta \left(\frac{k_2^2}{2m}-{\mu }_2^{\prime }\right)\right],{\tilde{n}}_3=\exp \left[-\beta \left(\frac{k_3^2}{2M}-{\mu }_3^{\prime }\right)\right],n=\exp \left[-\beta \left(\frac{q^2}{2M}-{ϵ}_0-{\mu }_2^{\prime }-{\mu }_3^{\prime }\right)\right],$$

$$N_e=2V\int \frac{d^{3}k}{(2\pi {)}^3}{\tilde{n}}_2(k)=2V{\left(\frac{m}{2\pi \beta }\right)}^{3/2}\exp (\beta {\mu }_2^{\prime }).$$

$$N_p$$

$$N_H$$

$$\frac{N_e{N}_p}{N_H}=V{\left(\frac{m}{2\pi \beta }\right)}^{3/2}\exp (-\beta {ϵ}_0).$$

$$L={\left(\frac{4\pi }{(2\pi {)}^3}\right)}^2\frac{1}{8Mm}{\int }_0^{\infty }d{k}_1{k}_{1}d{k}_2{k}_2^2\exp \left[-\frac{\beta }{2}\left(k_1+\frac{k_2^2}{2m}-{\mu }_2^{\prime }-{\mu }_3^{\prime }\right)\right]\delta \left(-{ϵ}_0+k_1-\frac{k_2^2}{2m}\right).$$

$$L={\left(\frac{4\pi }{(2\pi {)}^3}\right)}^2\frac{1}{8Mm}\exp [-\beta ({ϵ}_0-{\mu }_2^{\prime }-{\mu }_3^{\prime })/2]\frac{\sqrt{\pi }}{4}{\left(\frac{m}{\beta }\right)}^{3/2}\left({ϵ}_0+\frac{3}{4\beta }\right).$$

$$\Gamma =\frac{g^2\sqrt{\pi }}{8{\pi }^{3}M}{\left(\frac{m}{\beta }\right)}^{3/2}{ϵ}_0\exp (-\beta {ϵ}_0)\left(1+\frac{3}{4\beta {ϵ}_0}\right).$$

$$\Gamma \simeq \frac{e^2}{4\pi }\frac{1}{2{\pi }^{3/2}}\frac{m}{M}{\left(\frac{k_{B}T}{m}\right)}^{3/2}{ϵ}_0\exp (-{ϵ}_0/k_{B}T).$$

$$\rho (t)={\rho }_0\left[{\Omega }_M{\left(\frac{T}{T_0}\right)}^3+{\Omega }_R{\left(\frac{T}{T_0}\right)}^4\right],$$

$$H^2=H_0^2\left[{\Omega }_M{\left(\frac{T}{T_0}\right)}^3+{\Omega }_R{\left(\frac{T}{T_0}\right)}^4\right].$$

$$H=7.20\times 10^{-19}{T}^{3/2}{\left({\Omega }_M{h}^2+1.52\times 10^{-5}T\right)}^{1/2}{s}^{-1},$$

$${\Delta }_F(x,x^{\prime })=i⟨0⟩T\varphi (t,\mathbf{x})\varphi (t^{\prime },{\mathbf{x}}^{\prime })0,$$

$$D(x,x^{\prime })=i\frac{\text{mbox}Tr[\exp (-\beta H)\varphi (t,\mathbf{x})\varphi (t^{\prime },{\mathbf{x}}^{\prime })]}{\text{mbox}Tr\exp (-\beta H)}$$

$$D_{ij}(k_0,\mathbf{k})={\int }_{-\infty }^{+\infty }\frac{d{k}_0^{\prime }}{2\pi }\rho (k_0^{\prime },\mathbf{k}){\Lambda }_{ij}(k_0^{\prime },k_0),(i,j=1,2)$$

$$\rho (k_0^{\prime },\mathbf{k})$$

$$=2\pi ϵ(k_0^{\prime })\frac{1}{2\omega }\{\delta (k_0^{\prime }-\omega )+\delta (k_0^{\prime }+\omega )\},\omega =\sqrt{{\mathbf{k}}^2+m^2}.$$

$${\Lambda }_{11}(k_0^{\prime },k_0)=\frac{1+f(k_0^{\prime })}{k_0^{\prime }-k_0-i\eta }-\frac{f(k_0^{\prime })}{k_0^{\prime }-k_0+i\eta },$$

Full text

††journal: ApJ

1 Introduction

A century ago Saha (Saha, 1920) used the thermodynamics of chemical equilibrium to find the degree of ionization of atoms in a thermal bath. Since then it has been vigorously applied to investigate the spectrum of the sun and other stars. More recently, it has been used to find the epoch of hydrogen recombination involving the reaction

[TABLE]

in the early Universe, leading to an understanding of the cosmic microwave background (CMB) radiation observed today (Penzias & Wilson, 1965; Dicke, 1965). (There was an earlier epoch of Helium recombination, which can be treated separately in a first approximation.)

In the early Universe the equilibrium condition is not guaranteed a priori; it depends on the reaction rate and cosmic expansion rate. Also the Saha equation neglects the excited states of the H-atom, which is a very good approximation as long as thermal equilibrium prevails. But away from equilibrium the excited states may be important in the process. Accordingly a number of authors (Peebles, 1968; Zeldovich & Sunyaev, 1969) have investigated the hydrogen ionization and recombination in the realistic case, without assuming equlibrium and including the excited states (2s, 2p) of the H-atom along with the ground (1s) state. The discovery of CMB anisotropies prompted a resurge of the calculation of recombination including non-leading effects. Thus the earlier (effective) three level calculation was replaced with a multi-level one (Seager et al, 1999) including also Helium and their higher excited states. The effect of Raman scattering along with the related two-photon emission (Chluba & Thomas, 2011) was incorporated in the calculation. A review of all such effects is contained in the Karl Schwarzschild lecture by Sunyaev and Chluba (Sunyaev, 2018).

In this note, we attempt an indirect but easy way to study the effect of the excited states in attaining equilibrium. We calculate the ionization and recombination rates in a simplified model, where we exclude the excited states, taking only the H-atom in the ground state. Comparing these reaction rates with the cosmic expansion rate, we may know at which temperature the equilibrium is lost in the simplified model. On the other hand, comparing the fractional hydrogen ionization from the Saha equation with that of the realistic calculation (including the excited states), we can find when equilibrium is lost in the real world. Now comparing the two results for the loss of equilibrium, we may see the role of the exited states of H-atom.

We use the reaction rates to write the Boltzmann equation for an arbitrary (nonequilibrium) distribution of H-atoms. It has a simple analytic solution consisting of two terms (Weldon, 1983). The first term gives the equilibrium distribution, while the second term vanishes exponentially with time. It is the coefficient of time in the exponential of this term, which is identified with the reaction rate tending the distribution to equilibrium. This rate must be compared with the expansion rate of the Universe, given by the-then Hubble parameter.

Though the problem is non-relativistic, we shall use relativistic expressions to evaluate the rates and then apply non-relativistic approximation. We use units with $\hslash(=h/2\pi)$ and $c$ taken to be unity, where $h$ and $c$ are Planck’s constant and the velocity of light respectively. Also we write $\beta=1/k_{B}T$, where $k_{B}$ is the Boltzmann constant and $T$ the temperature.

2 Ionization and recombination rates

In atomic physics the processes represented by (1) are well-known. Recombination (or recapture) is the capture of an electron in the continuum by the atomic nucleus with the emission of photon. Ionization (or photoeffect) is the inverse process, where a photon is absorbed by an atom accompanied by ejection of an electron. These transitions are caused by the electric dipole operator in quantum mechanics (Bethe, 2008; Weinberg, 2013). It is given by the potential function, $V=-(e/m)\bf{A}\cdot\bf{p}$, where $e,m$ and $\bf{p}$ are the electric charge, mass of the electron and its momentum operator and $\bf{A}$ is the radiation field. The transition amplitudes for the above processes are given by the matrix elements of $V$ between two atomic states, where one is discrete and the other in the continuum.

Instead of using the above method to calculate the reaction rates, we use the elegant method of thermal quantum field theory for the problem  (Semenoff, 1983; Niemi, 1984; Mallik & Sarkar, 2016), where we directly get the recombination and ionization probabilities multiplied by appropriate factors involving distribution functions for particles in the medium. We first construct the interaction Lagrangian involving all the particles in (1). Let the photon, electron and proton fields be respectively $A_{\mu}(x)$, $\psi_{e}(x)$ and $\psi_{p}(x)$. We take H-atoms in the ground state only (ignoring its excited states), when it can be represented by an elementary scalar field $\Phi(x)$ (Weinberg, 1995; Dashen, 1974). Then the required effective interaction Lagrangian is

[TABLE]

It describes an electromagnetic interaction and $g$ must have mass dimension $-1$. So we take

[TABLE]

to within some uncertainty.

The calculation of the relevant reaction rate in thermal field theory can be best approached by considering the self-energy graph for the H-atom shown in Fig. 1. It includes in the intermediate state all other particles which appear in (1). The elements of real-time thermal field theory is sketched in Appendix A. These elements are used in Appendix B to calculate this graph, from which we obtain the different imaginary parts, which are collected in (B10) in a compact form. Any particular imaginary part may be isolated by integrating the variables $k_{0}^{\prime}$, $k_{0}^{\prime\prime}$ and $k_{0}^{\prime\prime\prime}$ in (B10) over the appropriate delta functions in the spectral functions, such as in (A3).

Before we write the desired imaginary part, we establish our notation. We shall not use the four-momentum notation anymore; instead, we now denote the magnitudes of three-momenta of H-atom (of mass $M_{H}$), photon, electron (of mass $m$) and proton (of mass $M$) by $q,~{}k_{i}~{}(i=1,2,3)$ and energies by $\omega=\sqrt{q^{2}+M_{H}^{2}},~{}\omega_{1}=k_{1},~{}\omega_{2}=\sqrt{k_{2}^{2}+m^{2}}$ and $\omega_{3}=\sqrt{k_{3}^{2}+M^{2}}$ respectively. We are interested in the imaginary part shown in Fig 2, corresponding to processes in (1) with the photon (and H-atom) incoming and electron and proton outgoing for ionization, when the reverse process (recombination) will automatically be given by the second term in bracket in (B10). Accordingly we choose the delta functions in $k_{0}^{\prime},~{}k_{0}^{\prime\prime}$ and $k_{0}^{\prime\prime\prime}$ variables as $\delta(k_{0}^{\prime}+\omega_{1})$, $\delta(k_{0}^{\prime\prime}-\omega_{2})$ and $\delta(k_{0}^{\prime\prime\prime}-\omega_{3})$. On using (A6) and (A11) to convert $f$ and $\widetilde{f}$ to $n$ and $\widetilde{n}$, we get for the required processes (denoted by subscript 1) from (B10),

[TABLE]

where $n_{1}(\omega_{1}),~{}\widetilde{n}_{2}(\omega_{2})$ and $\widetilde{n}_{3}(\omega_{3})$ are the equilibrium distribution functions for the photon, electron and proton respectively. It resembles the unitarity relation for the S-matrix in vacuum for the two-particle states of photon and H-atom. Compared to that relation, we now have the difference of two terms. Dividing it by $\omega$, we convert it to rates (Weldon, 1983)

[TABLE]

where $\Gamma_{d}$ and $\Gamma_{i}$ are the first and the second term in (4), representing respectively the decay and inverse decay rates of H-atom.

3 Boltzmann equation

The rates $\Gamma_{d}$ and $\Gamma_{i}$ are related. To see this, we write the equilibrium distribution functions explicitly as

[TABLE]

and similarly for factors involving $\widetilde{n}_{2}$ and $\widetilde{n}_{3}$. Taking into account the energy conserving delta function in (4), we get

[TABLE]

where $\mu_{2}$ and $\mu_{3}$ are the chemical potentials for the electron and proton and $L$ is given by

[TABLE]

where the product runs over $i=1,2,3$. For $i=1,~{}\mu_{1}=0$ and we take the upper sign and for $i=2,3$ we take the lower sign. So we get the ratio

[TABLE]

So far we have treated the H-atom as a single particle without any distribution in the medium. Let us now assume an arbitrary (non-equilibrium) distribution $n(\omega,t)$ of these particles. We can write a Boltzmann equation for $n(\omega,t)$, noting that it decreases at the rate $n\Gamma_{d}$ and increases at the rate $(1+n)\Gamma_{i}$ ,

[TABLE]

whose solution is (Weldon, 1983; Le Bellac, 2000)

[TABLE]

where $c(\omega)$ is an arbitrary function and we use (8). If $\mu$ is the chemical potential of the H-atom, the first term in $n(\omega,t)$ satisfies the condition

[TABLE]

which is the condition of chemical equilibrium (Reif, 1985), as can be read off from (1). Observe that this condition arises automatically in our calculation as a result of using the equilibrium thermal propagators. So the distribution function approaches the equilibrium value exponentially in time, irrespective of its initial distribution and its rate is governed by the reaction rate $\Gamma$.

So far the formulae are exact. Let us now make two simplifications appropriate for the problem at hand. First, we write the energies in nonrelativistic approximation:

[TABLE]

with $M_{H}=M+m-\epsilon_{0}$, where $\epsilon_{0}(=13.6~{}eV)$ is the binding energy of the H-atom in the ground state, which we neglect except in the exponential. Our second smplification results from the particle densities being dilute. We replace Fermi-Dirac and Bose-Einstein distributions with those of Maxwll-Boltzmann:

[TABLE]

where we define non-relativistic chemical potentials by $\mu^{\prime}_{2}=\mu_{2}-m$ and $\mu^{\prime}_{3}=\mu_{3}-M$.

The total number of electrons in volume $V$ is

[TABLE]

Similarly, the total number of protons and H-atoms are

[TABLE]

which incorporates the equilibrium condition (11). We then get the Saha equation (Weinberg, 2008)

[TABLE]

4 early Universe

We now examine the equilibrium condition in the early Universe. We first estimate the reaction rate $\Gamma$. As in the previous Section, we reduce the distributions to that of Maxwell-Boltzmann by retaining only the positive exponentials in the products in the denominator of (7). Also we ignore terms of $O(1/M)$ compared to $O(1/m)$ and write $\omega_{2}\omega_{3}\simeq mM$ in this denominator to get $L$ as

[TABLE]

We remove the $k_{1}$ integral with the delta function, when the $k_{2}$ integral reduces to a Gamma function, giving

[TABLE]

As the H-atom concentration is dilute, we set $1+n\simeq 1$ in (9), when the solution (10) shows the equilibrium distribution to be of Maxwell-Boltzmann type and the reaction rate becomes $\Gamma\simeq\Gamma_{d}$. From (6) and (19) we thus get

[TABLE]

Noting (3) we may rewrite it as

[TABLE]

Next, the expansion rate of the Universe is given by the Hubble pamameter $H=\dot{a}/a$, where $a(t)$ is the scale factor in the metric. In the era of interest to us ($T>1000K$), the constant vacuum energy is utterly negligible and we consider the energy density of matter, both non-relativistic ($\rho_{M}$) and relativistic ($\rho_{R}$). Also we do not include the curvature term, which, assuming a prior inflationary epoch, is driven to [math]. Denoting the present values by the subscript (0), we then write the total energy density as

[TABLE]

where $\Omega_{M}$ and $\Omega_{R}$ are fractions of the present critical energy density $\rho_{0}$. Applying Einstein equation with the spatially flat metric, it gives

[TABLE]

Including the contributions of photons and neutrinos in $\Omega_{R}$ and putting in numbers, one gets (Weinberg, 2008)

[TABLE]

where $h$ is the Hubble constant in units of 100km/sec/Mpc. We take $\Omega_{M}h^{2}=0.15$ as in Ref.(Weinberg, 2008).

In Table 1, we show the reaction rate ($\Gamma$) and the expansion rate ($H$) of the Universe at temperatures in the region of interest. It is seen that equilibrium prevails up to about 5000 K, if we include only the ground state of H-atom in calculating the reaction rate.

5 result and conclusion

In this work we do not calculate the fractional hydrogen ionization. We only estimate the role of excited states of H-atom in attaining the equilibrium condition as the temperature falls in the early Universe. We have already estimated above the temparature up to which equilibrium condition prevails in the case of our simplified model with no excited states of H-atom. We now get its estimate from the complete calculation of fractional ionization including the excited states and other physical effects presented in figure 3 of Sunyaev and Chluba (Sunyaev, 2018). Because the Saha equation assumes equilibrium condition, this condition should prevail as long as the complete calculation agrees with the Saha equation. We see that the two calculations start to disagree at redshift z $\sim 1500$ corresponding to T $\sim 4000$ K. Comparing this result with our calculation, we see that the excited states bring down the equilibrium temperature from $5000$ K to $4000$ K

Thus although the equilibrium number density of excited H-atoms is negligible compared to that in the ground state, they provide pathways to facilitate attaining the equilibrium and our calculation gives a quantitative estimate of this effect. Next, for results at low ionization levels, in which the last photon scattering took place, there is appreciable deviation from equilibrium (Peebles, 1968). Still, the Saha equation gives an order of magnitude estimate of the recombination temperature in the early Universe (Kolb, 1989).

Finally we comment on the use of (real time) thermal field theory. The factors involving the distribution functions in reaction probabilities in a medium are known, since Einstein introduced the $A$ and $B$ coefficients by considering detailed balance of equilibrium of atoms in the radiation field  (Einstein, 1917). As is well-known, they appear in matrix elements of creation and destruction operators in quantum field theory. These factors are now put in by hand in all processes taking place in a medium (Uehling, 1933; Weinberg, 1979). Here we show that they arise naturally originating from the thermal propagators. Another advantage of thermal field theory is that the polarization sums are done automatically in reaction probabilities.

6 Acknowledgement

One of us (S.M.) wishes to thank Prof. T. Souradeep for sending him his article entitled ‘Meghnad Saha and the Cosmic Photosphere’. He also thanks Prof. R. Rajaraman for a discussion.

Appendix A Equilibrium thermal field theory

Here we recall the basic elements of equilibrium thermal field theory (Mallik & Sarkar, 2016). Compared to vacuum field theory, it differs essentially in the time path and hence the propagator. To bring out this difference let us take a scalar field $\phi(x)$, with $x^{\mu}=(\tau,{\bf{x}})$, where the time variable $\tau$ may be complex and consider the time ordered propagator, which arises in perturbative calculations. While for the vacuum propagator

[TABLE]

the time variables run over the real time axis, the corresponding thermal propagator

[TABLE]

has the time variables running over an interval in the complex time plane. Here $H$ is the hamiltonian of the system and the thermal trace is over a complete set of states.

The time path at finite temperature may be broadly chosen in two different ways. It may be the imaginary segment from 0 to $-i\beta$ in the complex $\tau$ plane, giving the imaginary time formulation. In the real time formulation, which we shall use here, the time path must traverse the real axis. Then it must end at a point with Im $\tau=-\beta$. There are different ways to complete this path; we shall choose the one shown in Fig 3. Only the two horizontal lines contribute to the propagator, making the propagator a $2\times 2$ matrix.

It is possible to write spectral representation for propagator of fields of any spin. For the problem at hand, we have the scalar field $\Phi(x)$ to represent the H-atom, the vector field $A_{\mu}(x)$ for the photon and the Dirac field $\psi_{e}(x)$ and $\psi_{p}(x)$ for the electron and proton. The form of the spectral representation depends on the bosonic or fermionic nature of the field. The $2\times 2$ matrix propagator for the scalar field is

[TABLE]

where the spectral function is

[TABLE]

Here $\epsilon(k_{0})$ is the sign function defined as $+1$ for $k_{0}>0$ and $-1$ for $k_{0}<0$. As we shall see, only the $11$-component of the propagator appears in our calculation, for which

[TABLE]

where $f(k_{0}^{\prime})$ is a distribution-like function

[TABLE]

We shall use the $\delta$-function in the spectral function to express it in terms of the true distribution function $n(\omega)$,

[TABLE]

For the vector field $A_{\mu}(x)$, the propagator is again given by the one for the scalar propagator with the spectral function,

[TABLE]

For the spin $\frac{1}{2}$ (fermion) propagator, we have a similar representation

[TABLE]

with the spectral function

[TABLE]

and the 11-component of the $\Omega$-matrix is

[TABLE]

Again $\widetilde{f}(p_{0}^{\prime})$ can be written in terms of the fermion distribution function as

[TABLE]

(Also $\widetilde{f}(-\omega)=1-(\exp[\beta(\omega+\mu)]+1)^{-1}$, but we shall not need it.)

The above free propagators and the interacting (complete) ones and hence the self-energies may be diagonalised. Here we are primarily interested in (the imaginary part of) the self-energy of the $\Phi(x)$ field representing H-atom, which diagonalizes as

[TABLE]

where the diagonalizing matrix is

[TABLE]

From (A12) we get

[TABLE]

so that we may evaluate only the 11-component of the $\Sigma$ matrix to get the imaginary part of the diagonalized matrix.

Appendix B Evaluation of self-energy graph

Here we shall evaluate the two-loop thermal self-energy graph (Fig 1) of the H-atom. As the vacuum and medium calculations differ only in the propagators, we can first conveniently find it in vacuum:

[TABLE]

where tr indicates trace over $\gamma$ matrices and $D^{\mu\nu}$, $S^{(e)}$ and $S^{(p)}$ are the vacuum propagators for the photon, electron and proton respectively. Then the 11-component of the thermal self-energy matrix is immediately obtained by replacing the vacuum propagators with the corresponding $11$-component of the thermal propagator matrices:

[TABLE]

where we keep only the independent loop momenta. We write the propagator in their spectral representations (A) and (A8). The tensor and spinor factors in the spectral functions can be collected to give

[TABLE]

Removing these factors, we get the three spectral functions as

[TABLE]

where $k_{3}=q-k_{1}-k_{2}$. Also we segregate the integrals in energy components of $k_{1}$ and $k_{2}$ over the energy denominators of propagators. We thus write (B2) as

[TABLE]

where

[TABLE]

with $f^{\prime}=f(k_{0}^{\prime})$, $\widetilde{f}^{\prime\prime}=\widetilde{f}(k_{0}^{\prime\prime})$, $\widetilde{f}^{\prime\prime\prime}=\widetilde{f}(k_{0}^{\prime\prime\prime})$.

Let us work out the integral over $k_{10}$ first, noting that this variable appears in the first and third factors in (B6). When these two factors are multiplied out, these result four terms. As the integral converges in both the upper and lower half of the $k_{0}$ plane, we can evaluate it closing the integration contour in either half. So only two of these terms, having poles in both the upper and lower halves of $k_{10}$ plane, can contribute to the integral. Thus we evaluate the $k_{10}$ integral in (B6) to get

[TABLE]

Next carry out the $k_{20}$ integral in the same way over the second factor in (B6) and the one just obtained to get

[TABLE]

giving its imaginary part as

[TABLE]

Comparing with (A14) we finally get the imaginary part of the diagonalised self-energy as

[TABLE]

Integrating over $k^{\prime}_{0},k^{\prime\prime}_{0}$ and $k^{\prime\prime\prime}_{0}$ with the delta functions contained in the three spectral functions, we get eight terms corresponding to different particles in the initial and final states. In Section 2 we get one of these terms representing the ionization and recombination probabilities.