Basic Recombination and Formation of

Neutral Hydrogen in Early Universe

A simple study Recombination in Early Universe
Self-Study
by
Kazi Abu Rousan

School of Physical Sciences

Indian Association for The Cultivation of Science,
Kolkata - 700032, India
December 2023
 What this article discuss?
In this article we have studied process of Recombination, i.e., Formation of Neutral Hydrogen.
This study has been done using simplified equilibrium condition and then in Non-equilibrium
condition. For the Non-Equilibrium case we have used the effective threel-level model.

For each case we have studied the how Free Electron Fraction evolve as a function
of temperature and redshift. Finally we calculated the value of the redshift at which the
recombination took place.

i
Contents
1 Introduction 1
1.1 Local Thermal Equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Decoupling and Freeze Out . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Recombination in Equilibrium Condition 5

2.1 Distribution Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Macroscopic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 The Primordial Plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.4 Hydrogen Recombination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Recombination beyond equilibrium 9

3.1 Boltzmann Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Collision Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 Hydrogen Recombination 11
4.1 Selection Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2 Effective Three-Level Atom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5 Conclusion 16

6 Appendix 17
6.1 Hubble’s Constant for Matter-Radiation Equality . . . . . . . . . . . . . . . . . 17
6.2 Codes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1 Introduction
My logic is very simple. If we study something, It must have some goal. This let me ask myself,
Why do we need to invest time in studying Recombination?
The answer for my case is very straight forward. One of the most important and powerful
tool to study early universe is Cosmic Microwave Background Radiation (CMBR) and it’s
fluctuations. Recombination is the exact method which is responsible for CMBR. Hence,
we have to know how recombination happened to understand CMBR. But that’s not all. To
study the fluctuations, we need to know when they entered the horizon. As depending upon
the epoch of that time, the fluctuations increase or decrease in different forms. Hence, to study
the fluctuations, we need to know when this recombination happened. Also, we need to know
what was the electron density at the end of recombination. (Don’t worry if some of the terms
sound like too much, we will discuss those in some later article).
With this, let’s start our study.

1.1 Local Thermal Equilibrium

We know the universe is expanding as time goes on. It was hotter and denser in it’s earlier days.
Particles which were present used to collide frequently and they were in thermal equilibrium
with some temperature T . At these times all the thermodynamical properties were determined
by local equilibrium. Although these equilibrium makes our life simple from calculation point
of view, it is the departure from the equilibrium that makes things interesting and helps in
large structure formations. The Non-equilibrium dynamics allows massive particles to acquire
cosmological abundances.
The most important point in understanding the thermal history of the universe is the com-
parison between the rate of interaction of interest Γ and the rate of expansion of the
universe H. Whenever the time scale of particle interaction is much smaller than the
characteristic expansion time scale, i.e.,
1 1
tc ≡ << tH ≡ (1)
Γ H
local thermal equilibrium is reached. As the universe cooled down due to expansion, tc ∼ tH ,
i.e., interaction rate decreases compared to expansion rate and things go out of equilibrium. At
tc ∼ tH , particles decouple from the thermal bath. As there are different species with different
interaction rates, the decoupling times are also different.
Let’s first review all the important topics which will be needed for our discussion.
Local Thermal Equilibrium (L.T.E) is of the most important keypoints in our study as
It is exactly the point which tells us when things are in equilibrium and also let’s us find the
temperature range.
As an example we can show that the condition eq n −1 is satisfied for Standard Model pro-
cesses at temperatures above a few hundred GeV. We write the rate of particle interactions
as,
Γ = nσv (2)
where n is the number density of particles, σ is their interaction cross section, and v is the
average velocity of the particles. Let’s consider a process of the form 1 + 2 ↔ 3 + 4, we would
write the interaction rate of species-1 as Γ1 = n2 σv, where n2 is the density of the target
species-2 and v is the average relative velocity of 1 and 2.
The interaction rate of species-2 would be Γ2 = n1 σv. We have used the expectation that at
high energies n1 ∼ n2 ≡ n. For T ≥ 100GeV , all known particles are ultra-relativistic, and
hence v ∼ 1. Since particle masses can be ignored in this limit, the only dimensionful scale
is the temperature T . Dimensional analysis then gives n ∼ T 3 . Interactions are mediated by
gauge bosons, which are massless above the scale of electroweak symmetry breaking.

1
Figure 1: Feynman diagram associated with a 2 → 2 scattering process
√ mediated by √ the √
exchange of a gauge
boson. Each vertex contributes a factor of the gauge coupling gA ∝ α. Here σ ∝ ( α × α)2

The cross sections for the strong and electroweak interactions then have a similar dependence,
which also can be estimated using dimensional analysis,
α2
σ ∼ | fig-1|2 ∼ (3)
T2
g2
where α = 4πA is the generalized structure constant associated with the gauge boson-A.
Using this we find,
α2
Γ = nσv ∼ T 3 × 2 × 1 = α2 T (4)
T
We also know that, √
ρ T2
H∼ ∼ (5)
MP l MP l
n n
Taking the ratio of eq −4 and eq −5, we get,
Γ α2 MP l 1016 GeV
∼ ∼ (6)
H T T
Here we have used α ∼ 0.01. This clearly showed that below T ∼ 1016 GeV but above
100GeV , eq n −1 is satisfied.
When particles exchange energy and momentum efficiently they reach a state of
maximum entropy. It is a standard result of statistical mechanics that the number of particles
per unit volume in phase space, i.e., the distribution function takes the form,
1
f (E) = (7)
exp( E−µ
T
) ±1
Here + and − corresponds to fermions and bosons respectively. When the temperature drops
below the mass of the particles i.e., T << m, they become non-relativistic and their dis-
−m
tribution function gets an exponential suppression (f ∼ e T ) . This means that relativistic
particles (mainly radiation) dominate the density and pressure of the primordial plasma. The
total P
energy density is ρr therefore well approximated by summing over all relativistic particles,
ρr ∝ i d3 pfi (p)Ei (p). The result can be written as,

π2
ρr = g∗ (T )T 4
30
where g∗ (T ) is the number of relativistic degrees of freedom. Fig-2 shows the evolution
of g∗ (T ) assuming the particle content of the Standard Model. At early times, all particles
are relativistic and g∗ = 106.75. The value of g∗ decreases whenever the temperature of the
universe drops below the mass of a particle species and it becomes non-relativistic. Today, only
photons and (maybe) neutrinos are still relativistic and g∗ = 3.38.

2
 Figure 2: Evolution of the number of relativistic degrees of freedom assuming the Standard Model.2

1.2 Decoupling and Freeze Out

As we have discussed in the previous section all were in equilibrium but if this equilibrium had
presisted until today then the universe would be mostly photon. Any massive particle species
would be exponentially supressed. Hence to understand how everything is been created it is
very important to understand the deviation from equilibrium that led to the freeze-out of the
massive particles.

Figure 3: A rough sketch of particle freeze-out. For T >> m, the particle abundance follow equilibrium and
for T << m the particles freeze-out and maintain a density which is much larger than suppressed boltzmann
equilibrium abundance.

Using the method used in the previous section, we can again find the temperature of freeze-
out. Below the scale of electroweak symmetry breaking, T ≲ 100GeV , the gauge bosons of the
weak interactions, W ± and Z, recieve masses MW ∼ MZ . The cross section associated with
this process mediated by the weak force becomes,1
σ ∼ G2F T 2
where GF is Fermi’s constant which is GF ∼ Mα2 ∼ 1.17 × 10−5 GeV −2 . As we can see the
W
strength of the weak interactions decreases as the temperature of the universe drops. Using
this,
Γ α2 T 4 MP l  T 2
∼ 2
× 2 ∼ (8)
H MW T 1M eV

3
 Events Redshift(z) Temperature
Inflation ? ?
Baryogenesis ? ?
EW Phase Transition 1015 100 GeV
QCD Phase Transition 1012 150 MeV
Dark Matter Freeze-out ? ?
Neutrino Decoupling 6 × 109 1 MeV
Electron-Positron Annihilation 2 × 109 500 keV
Big Bang Nucleosynthesis 4 × 108 100 keV
Matter-Radiation Equality 3400 0.75 eV
Recombination 1100 − 1400 0.26 − 0.33 eV
Photon Decoupling 1000 − 1200 0.23 − 0.28 eV

Table 1: A rough sketch of important events upto Recombination.

which shows that the ratio drops below 1 at Tdec ∼ 1M eV . Particles that interact with the
primordial plasma only through the weak interaction hence decouple 1M eV . This decoupling
of weak scale interactions has importance for the thermal history of the universe. A rough
sketch of all the important events are shown in the table-1.
In our work we have studied this recombination. In this process Neutral Hydrogen forms
through the reaction e− + p+ → H + γ when the temperature has become low enough that the
reverse reaction is energetically not favoured.
Before recombination happened the strongest coupling between the photons and the rest of
the plasma soup is through Thomson Scattering e− + γ → e− + γ. The drop in the free
electron density after the recombination means this process becomes insufficient and photons
decouple. They stream through the universe freely since then. They are the exact photons we
observe as the CMBR.

4
2 Recombination in Equilibrium Condition
We have observational evidence from the perfect blackbody spectrum of Cosmic Microwave
Background (CMB) that the early universe was in local thermal equilibrium. It should
be noted that the universe can never truly be in equilibrium state as the FRW spacetime
doesn’t posses a time-like killing vector. But if the expansion of the universe is slow enough,
then particles have enough time to settle close to local equilibrium. And as the universe is
homogeneous, the local values of thermodynamics quantities are also global values. Hence,
It is reasonable to study the Recombination in equilibrium condition.
To study this, we will first review some ideas of basic equilibrium thermodynamics and
statistical mechanics.

2.1 Distribution Function

As we have already discussed, the probability that a particle chosen at random has a momentum
p⃗, in general is very a complicated distribution function f (⃗p, t). This function can be a function
of position ⃗x. However, in a homogeneous universe such a dependence is not possible due to
the translational invariance. Also, isotropy requies that the momentum dependence is only in
terms of the momentum magnitude p ≡ |⃗p|.
Now, if we consider time long enough (with respect to (w.r.t) the typical interaction timescale),
then the system will reach equilibrium and is characterised by a time-independent distribution
function given by (eqnn −7),i.e.,
1
f (p, T ) = (9)
exp( E(p)−µ
T
)±1

2.2 Macroscopic Properties

If we want to relate the macroscopic description, we must sum over all possible momentum
states of the particles weighted by their probabilities. To do this using distribution function we
need density of state.
If a particle has g internal degree of freedom then the density of states becomes,
g g
3
= where ℏ = 1 is used (10)
h (2π)3
Using this we can write the number density of particles,
Z
g
n(T ) = d3 pf (p, T ) (11)
(2π)3
In a similar manner energy density is,
Z
g
ρ(T ) = d3 pf (p, T )E(p) (12)
(2π)3
For finding the pressure we have to do a little more work.
Consider a small area element of size dA, with unit normal vector n̂ (see fig-4).

5
 Figure 4: Pressure in a weakly interacting gas of particles.

All particles with velocity v. Striking this area element in the time interval between t and
t + dt, were located at t = 0 in a spherical shell of radius R = vt and width vdt. A solid angle
dΩ2 of this shell defines the volume dV = R2 vdtdΩ2 . Multiplying the phase space density by
dV gives the number of particles in the volume (per unit volume in momentum space) with
energy E(v).
g
dN = 3
f (E)R2 vdtdΩ
(2π)
Not all particles in dV reach the target, only those with velocities directed to the area element.
Taking into account the isotropy of the velocity distribution, we find that the total number of
particles striking the area element dAn̂ with velocity ⃗v = vv̂ is,
v̂ · n̂dAdN gf (E) ⃗v · n̂
dNA = 2
= dAdtdΩ
4πR (2π)3 4π
where ⃗v · n̂ < 0. If these particles are reflected elastically, each transfer momentum 2⃗p · n̂ to the
target. Hence, the contribution to the pressure by the particles of velocity v is,
2⃗p · n̂ p2 p2
Z Z
gf (E) 2 gf (E)
dP (v) = dNA = × cos (θ) sin(θ)dθdϕ = ×
dAdt (2π)3 2πE (2π)3 3E
where v = Ep and integrated over the hemisphere defined by v̂ · n̂ = − cos(θ) < 0.
Hence, the Pressure is,
p2
Z
g 3
p
P (T ) = d pf (p, T ) where E(p) = m2 + p2 (13)
(2π)3 3E(p)
Each particle species i with mi , µi , Ti will have it’s own function fi and hence will have different
density(ni ), energy density(ρi ) and pressure(Pi ). But they can share same temperature Ti = T
if they are in thermal equilibrium. Hence, their densities and pressure can then only differ
because of differences in their masses and chemical potentials.
At early times, the chemical potentials of all particles are much smaller than the
temperature µi << T , and hence can be neglected.

2.3 The Primordial Plasma

Using the previous discussion the state of the early universe in thermal equilibrium. Using
eq n −11, Z ∞
g p2 dp
n= 2 p (14)
2π 0 exp( p2 + m2 /T )exp(− Tµ ) ± 1

6
 ∞
p
p2 p2 + m2 dp
Z
g
ρ= 2 p (15)
2π 0 exp( p2 + m2 /T )exp(− Tµ ) ± 1
Defining new dimensionless variables x ≡ m
T
, ξ ≡ Tp and y ≡ Tµ , this gives us,
gT 3 ∞ ξ 2 dξ
Z
n= p (16)
2π 2 0 exp( x2 + ξ 2 )exp(−y) ± 1
similarly,
∞
p
gT 4 ξ 2 ξ 2 + x2 dξ
Z
ρ= p (17)
2π 2 0 exp( x2 + ξ 2 )exp(−y) ± 1
To go further we will take two cases seperately, i.e., Relativistic and Non-relativistic case.

Relativistic Limit

At higher temperature compared to the mass of the particle, we can use x → 0. This gives(see
fig-10),
gT 3
nf ermi = − 2
× Li3 (−ey ) for Fermions
π (18)
gT 3
nboson = 2 × Li3 (ey ) for Bosons
π
For energy density(see fig-10),
3gT 4
ρf ermi = − × Li4 (−ey ) for Fermions
π2 (19)
3gT 4
ρboson = × Li4 (ey ) for Bosons
π2
For species with 0 chemical potential, i.e., µ = 0 → y = 0, we have,
ζ(3)gT 3 3
n= × for Fermions
π2 4 (20)
ζ(3)gT 3
n= for Bosons
π2
and
π 2 gT 4 7
ρ= × for Fermions
30 8 (21)
π 2 gT 4
ρ= for Bosons
30
Using these values, the total energy density of the early universe ρ is given by summing over all
species. But it is also common to write the density in terms of temperature of the universe T (
which is taken as photon temperature Tγ ). Using these informatiuons, we can define Effective
Number of Degrees of Freedom at temperature T by,
X  Ti 4 7 X  Ti 4
g∗ (T ) ≡ gi + gi (22)
i=b
T 8 i=f
T

Non-relativistic Limit

For temperature smaller than the particle mass, i.e., x >> 1, we will have,
 mT  23 µ − m
n=g exp for both boson and fermion (23)
2π T

7
2.4 Hydrogen Recombination
The Recombination needs two steps. Initially helium atoms formed and after that forms
hydrogen. For simplifying our analysis, we will consider that the universe was filled only with
free electrons, protons and photons. It is a very good approximation as over 90% of the nuclei
are protons. Also, helium recombination was completed before hydrogen recombination, hence
these two can be studied searately.
Hydrogen atoms forms using,
e− + p+ ↔ H + γ (24)
Initially, this reaction keeps the particles in equilibrium and as T < mi for i = e, p and H,
we can use eq n −23. For the given reaction-24, we know that the chemical potentials follow
µp + µe = µH . Using these, we have
 n  gH  mH 2π 3/2 (mp +me −mH )/T
H
= e (25)
ne np eq ge gp me mp T
We will know use mH ≈ mp for the prefactor but can’t do the same for the exponent as exponent
grows much faster and also defining EI ≡ mp + me − mH = 13.6eV , we have,
 n  n   2π 3/2
H H
= = eEI /T (26)
ne np eq n2e eq me T
Here we have used gp = ge = 2 and gH = 4. The spins of the electron and proton in a hydrogen
atom can be aligned or anti-aligned, giving singlet and triplet, so gH = 1 + 3 = 4. Also, as the
universe is neutral, we have used ne = np .
It is easier to use free electron fraction which is defined by,
ne ne
Xe ≡ = (27)
np + nH ne + nH
A funnly ionized universe then corresponds to Xe = 1, while a universe of only neutral atoms
has Xe = 0.
Like already mentioned small helium atom can be neglected. But we can include it’s amount
by considering ne + nH = 0.76nb .5 Then we can use (using eq n −20),
2ζ(3) 3
nb = ne + nH = ηnγ = η T (28)
π2
where η is the baryon to photon ratio. Using this we get,
1 − Xe nH
2
= 2 0.76nb
Xe ne
Using eq n −26, we get Saha Equation,
1 − X  1.52ζ(3)  2πT 3/2 EI /T
e
= η e (29)
Xe2 eq π2 me
The solution is given by,
√
−1 + 1 + 4f 1.52ζ(3)  2πT 3/2 EI /T
Xe = where f (T, η) = η e (30)
2f π2 me

8
 Figure 5: Free Electron Fraction(Xe ) as a function of z for Saha Equation.

The recombination temperature is defined as the temperature at which Xe = 0.5 in

eq n −29. For η = 6 × 10−10 ,we have,

Trec ≈ 0.3217eV
An interesting point to note is that the recombination temperature is significantly below the
binding energy of hydrogen. This is the case as there are many photons for each hydrogen
atom.
Using T = 0.2329 meV (1 + z),2 gives us zrec = 1369.
The details for recombination are little more complex and the moment of recombination is
delayed relative to the Saha prediction, with Xe = 0.5 at zrec ≈ 1270.
But that’s not all. In fig-5, we can see that around z = 1000, Xe = 0, i.e., there should not
be any free electron after z = 1000 but it should be obvious that it’s completely wrong.
Now, let’s see what are the changes for non-equilibrium condition.

3 Recombination beyond equilibrium

To truely understand structue formulasim we need to understand deviations of the systems from
equilibrium. The best and main tool for non-equilibrium study is Boltzmann Equation. We
will start our analysis by studying boltzmann equation first.

3.1 Boltzmann Equation

In the absense of interactions, the number density of the particle species i evolves as,
dni
+ 3Hni = 0 (31)
dt

9
This simply imples that the number of particles in a fixed physical volume (V ∝ a3 ) is
conserved, so that the density decreases as ni ∝ a−3 . To include the effects of interactions, we
add a collision term to the right-hand side of the equation,1
d3 p
Z
dni (t)
+ 3Hni (t) = C[f ] (32)
dt (2π)3
The collision term which include the collisions can change behaviour of the inetgral over the
mometum of the collision. Here C[f ] ≡ Ci [nj ], i.e., collision term depend on the distribution
function. The form of the collision term depends on the specific interaction under consideration.
Let’s try to find the term.

3.2 Collision Terms

Collision include scattering as well as pair creation, annihilation and particle decay. Let’s try
to find the term for 2 → 2 reaction, i.e.,
(1)p + (2)q ↔ (3)p′ + (4)q′ (33)
Here subscripts indicate momentums. These process follow conservation of energy and momen-
tum i.e.,
p⃗ + ⃗q = p⃗′ + ⃗q′ ; E1 (⃗p) + E2 (⃗q) = E3 (⃗p′ ) + E4 (⃗q′ ) (34)
p
where Ei (p) = p2 + m2i denotes the energy-momentum relation for ith particle. Each type
of particle has it’s reprective distribution function fi (⃗x, p⃗, t), i = 1, 2, 3, [Link] collision term for
the species 1 can be written as (see fig-6),
p
⃗+⃗ p′ +⃗
q =⃗ q′
X
C[f1 (p)] = 1
δD (E1 (⃗p) + E2 (⃗q) − E3 (⃗p′ ) − E4 (⃗q′ ))|M|2 (f3 (⃗p′ )f4 (⃗q′ − f1 (⃗p)f2 (⃗q)) (35)
q ′ ,⃗
q⃗,⃗ p′

where |M|2 is the scattering amplitude squared which depends on the details of the interactions,
which can be calculated using feynman’s diagrams. In the eq n −35 we have neglected Pauli
Blocking and Bose enchancement as in our case they are really not that relevent.

Figure 6: Effect of the collision for 2 → 2 [Link] forward reaction removes particles from f1 in proportion
to f1 f2 and increases in proportion to f3 f4 .

10
 Using the previous discussion for continious parameters the collision term will be,
d3 q d3 p′ d3 q ′
Z Z Z
1
C[f1 (p)] = |M|2
2E1 (p) (2π)3 2E2 (q) (2π)3 2E2 (p′ ) (2π)3 2E2 (q ′ )
(36)
×(2π)4 δD 1
(E1 (⃗p) + E2 (⃗q) − E3 (⃗p′ ) − E4 (⃗q′ ))δD
3
(⃗p + ⃗q − p⃗′ − ⃗q′ )
×(f3 (⃗p′ )f4 (⃗q′ − f1 (⃗p)f2 (⃗q))
Combining all these ideas, we can write the Boltzmann Equation for the 2 → 2 process as,
1 d(n1 a3 ) h n n 
1 2
i
= −⟨σv⟩ n 1 n2 − n3 n 4 (37)
a3 dt n3 n4 eq
Here ⟨σv⟩ is the thermally averaged cross section, which is given by,
d3 p d3 q d3 p′ d3 q ′
Z Z Z Z
⟨σv⟩ = e(E1 +E2 )/T
(2π)3 2E1 (p) (2π)3 2E2 (q) (2π)3 2E2 (p′ ) (2π)3 2E2 (q ′ ) (38)
×(2π)4 δD
1
[E1 (⃗p) + E2 (⃗q) − E3 (⃗p′ ) − E4 (⃗q′ )]δD
3
(⃗p + ⃗q − p⃗′ − ⃗q′ )|M|2

Using this for our case of e− + p+ ↔ H + γ, we can show,

 α 2  E 1/2  E 
s I I
⟨σv⟩ ≈ 9.8 ln
me T T
1
where EI = 13.6eV and αs = 137 .
Don’t worry about this as we will directly not use this, rather we will use simple arguments to
get the non-equilibrium rate equation.

4 Hydrogen Recombination
As the free electron fraction drops during recombination, the interation rate also decreases and
can also drop below the expansion rate of the universe. As this happens like discussed before,
our equiulibrium treatment breaks down.
Using eq n −37 for the reaction eq n −24, we get
1 d(ne a3 ) h
2
 n2 
e
i
= ⟨σv⟩ ne − n H (39)
a3 dt nH eq
where we have used ne = np and nγ = neq γ . To go any further we have to study the process by
whihc electron recombine in detail. Below we will do the exact thing in detail. This will help
us how peeble added a simple correction term in the boltzmann equation.
One very important point to mention is that we can just take this equation and solve it
numerically to get some result. But it doesn’t give us correct value although it does remove
some previous errors.
To get the correct answer, we have to study the exact process of recombination, which introduces
peeble’s factor in font of the ⟨σv⟩ term.

4.1 Selection Rule

A hydrogen atom can only make a spontaneous transition from an energy state corresponding
to the quantum numbers n′ , l′ and m′ if the modlulus squared of the associated dipole moment
is nonzero,

d2 = |⟨n, l, m|ex|n′ , l′ , m′ ⟩|2 + |⟨n, l, m|ey|n′ , l′ , m′ ⟩|2 + |⟨n, l, m|ez|n′ , l′ , m′ ⟩|2

11
The spontaneous transition between different energy levels of a hydrogen are only possible
provided,
l′ =l ± 1
(40)
m′ =m ± 0, 1
Similarly, for quadrupole transition the selection rule are,
l′ =l ± 0, 2
(41)
m′ =m ± 0, 1
To from a hydrogen, proton can capture an electron. The photon that is emitted when a free
electron is captured by a proton into the ground state has more than enough energy to ionize
another hydrogen atom, so this process produces no net decrease in ionization. Similarly, the
photon emitted when an electron in a higher orbit of the hydrogen atom with principle quantum
number n ≥ 3 falls into the ground state has more than enough energy to lift an electron in
the ground state of some other hydrogen atom to the n = 2 excited state, so this also produces
no net increase in the number of electrons in the ground state.
To avoid these processes, recombination must first occur to an excited state, which then
decays to the ground state. The photons created in this multistep process have lower energy
and are therefore less likely to ionize the plasma. Peeble first argued that we can thing this
process as an effective three-level atom.6

Figure 7: Diagram for effective three-level atom.

4.2 Effective Three-Level Atom

Introduction

As mentioned Peeble first argued that the hydrogen atom can be treated as an effective three
level atom. The three relevant states are the ground state(1s), the excited state(2s and
2p) and the continuum states of ionized hydrogen (fig-7). The excited states are in thermal
equilibrium with each other (any state for which n ≥ 3) since radiative excitation and decays
are very fast. Hence they all can be considered as a single entity.
The ground state of hydrogen is typically reached by formation of excited states H ∗ in the
reaction e + p → H ∗ + γ, followed by a radiative decays down to the n = 2 excited state.

12
 When an atom in the first excited state decays to the ground state it produces a Lyman-α
photon. This photon has a large probability to be absorbed by a nearby atom, which can then
be ionized. But as the universe expands, the energy of these photons get redshifted. So, if it
doesn’t interact very soon with another atom, then cosmological redshift will take it’s energy
outside the resonant line.

Finding Boltzmann Equation

The analysis of these three level can be done using three basic approximation,5, 6 which were
made by Peeble. We will study them in detail.
1. (Continuum State) Collision between hydrogen atoms and radiative transitions between
the states of these atoms are sufficiently rapid so that all states of the atoms are in
equilibrium with each other at temperature T of the radiation, except for the 1s ground
state, which reaches only through slow or inefficient processes.
This has the consequence that the number density of the state n, l of the hydrogen can be
expressed in terms of the number density of any one state (we will take 2s),
nnl = (2l + 1)n2s exp((E2 − En )/T ) (42)
where En is the ionization energy of state with principle quantum number n. 2s is chosen
as peeble’s effective three-level model assume only that is the distinciable higher state
avalible.
2. The net rate of change in the population of hydrogen atoms in their 1s state is given by
the rate of radiative decay from the 2s and 2p states, minus the rate of excitation of these
states from the 1s state. Following the above point, all other process leading to the ground
state are assumed to be canceled by the reionization or reexcitation of other atoms by the
emitted photon.
Recombnation happens in collisions of electrons and protons, hence it decreases the number
ne a3 of free electron in a co-moving volume a3 at a rate α(T )np ne a3 .
Leaving aside ionization from the ground state which just cancels recombination to the
ground state, the ionization from excited states increases ne a3 at a rate given by a sum of
terms ∝ nnl a3 with n > 1. As eq n −42 gives all the nnl with n > 1 in terms of n2s , hence
ionization increases the number of electrons in co-moving volume a3 at a rate that is given
by β(T )n2s a3 . Hence, we have
d(ne a3 )
− αn2e a3 + βn2s a3
dt
P
Dividing both side by n = np + nH = np + n,l nnl = 0.76nb , we get,
d  ne  αn2 βn2s
=− e + (43)
dt n n n
This must vanish in equilibrium (for the rapid process e + p ↔ 2s) and hence (using
eq n −23),
β(T )  n2e   m T 3/2
e
= = e−E2 /T (44)
α(T ) n2s eq 2π
Now, we have to relate n2s with the total number density n of protons and hydrogen atoms.
3. The total number of excited hydrogen atoms in a co-moving volume changes so much slowly
than individual radiative process that the net increase in this number due to recombination
and reionization of hydrogen is balanced by the net decrease in this number by transitions
to and from the 1s state, i.e.,
αn2e − βn1s = (Γ2s + 3PΓ2p )n2s − εn1s (45)

13
 where Γ2s and Γ2p are the rates for the radiative decay processes 2s → 1s + γ + γ and
2p → 1s + γ. As 2p → 1s + γ + γ is much slower compared to 2s → 1s + γ + γ hence we
have neglected the process. P is the probability that th Lymann α photon emitted in the
decay 2p → 1s + γ will escape to infinity. Also we have used the fact n2p = 3n2s .
Now, we also have nH = n1s + n2s + n2p = n1s + 4n2s . This gives n1s = (nH − 4n2s ). From
these equations, we get,
αn2e + εnH
n2s = (46)
Γ2s + 3PΓ2p + β + 4ε
Using this along with the fact that at equilibrium 1s state would be constant, we have
ε n 
2s
= = e(E2 −EI )/T (47)
Γ2s + 2PΓ2p n1s eq
Using these with the fact that in our temperature range 1 + 4 exp(E2 − EI ) ≈ 1 and also
using Xe = ne /n = 1 − nH /n we have,
dXe Γ2s + 3PΓ2p  nα  m T 3/2
e

= − αnXe2 + (1 − Xe ) e−EI /T (48)
dt Γ2s + 3PΓ2p + β 0.76nb 2π
Writing these in a neat way, we have,
dXe h
2
i
= Cr b(T )(1 − Xe ) − 0.76nb α(T )Xe (49)
dt
where
 α 2  E 1/2  E 
s I I
α(T ) =⟨σv⟩ ≈ 9.8 ln
me T T
 m T 3/2
e
b(T ) =α(T ) e−EI /T
2π (50)
Γ2s + 3PΓ2p
Cr =
Γ2s + 3PΓ2p + β
β =b(T )e3EI /4T
where αs ≈ 1/137. For α(T ) we can use some QFT to find that.1, 2, 5 The eq n −49 is the Boltz-
mann Equation of recombination which is called Peebles Equation.
Hence, we finally have in conclusion,
• There is a small probability that the 2s state will decay to the 1s state through the emission
of two photons. Hence, these two photons don’t have enough energy to excite the atoms
in the ground state back to the first excited state. The rate for this two-photon decay
is,
Γ2s = 8.227s−1 = 5.4134 × 10−15 eV
• As mentioned before Lyman-α photons are emitted in 2p → 1s transition. These has a
large resonant cross section for exciting another hydrogen. But as the universe expands,
the energy of these photons get redshifted. So, if it doesn’t interact very soon with another
atom, then cosmological redshift will take it’s energy outside the resonant line. Here the
rate of recombination via resonance escape is,
8πH
3PΓ2p =
λ3α n1s
where λα ≡ 2π/Eα = 8π/(3EI ) is the lyman-α photon’s wavelength and also n1s ≈
0.76(1 − Xe )nb is the abundance of hydrogen atoms in the 1s state. Using these, we get,
 27E 3  H(T )
I
3PΓ2p = where nbH = 0.76nb (51)
64 nbH (1 − Xe )

14
Solving Peeble’s Equation

Using the evolution of the Hubble Parameter is,

p  1 + z 1/2
H(z) = Ωm H0 (1 + z)3/2 1 + (52)
1 + zeq
Ωm = 0.3089, H0 = 1.5 × 10−33 eV and zeq = 3395. It is convenient to use redshift instead of
time, shown as,
dXe dXe dz dXe
= =− H · (1 + z)
dt dz dt dz
Using this we get (peeble’s equation in terms of z),
dXe Cr (z) 1 h i
= × α(T )nbH Xe2 − b(T )(1 − Xe ) (53)
dz H(z) (1 + z)
As previously mentioned, the equation without Cr is Boltzmann Equation. This doesn’t
take into account the exact process (i.e, the redshift of α photon and the two-photon decay).
The equation is,
dXe 1 h
2
i
= α(T )nbH Xe − b(T )(1 − Xe ) (54)
dz H(z)(1 + z)
Solving these two using a numerical method (see Appendix for code), we get (see fig-8).

Figure 8: Comparison of Saha, Boltzmann and Peeble Equation. Here Boltzmann Equation means the peeble’s
equation without Cr factor in eq n −53.

This shows us the deviation from equilibrium. This also shows us that the free electron
fraction is not zero even at z = 0.
The recombination redshift and temperature for boltzmann eq n −54 is ,
Trec ≈ 0.32eV and z = 1368

And for recombination redshift and temperature for peebles eq n −53 is ,

Trec ≈ 0.298eV and z = 1270

15
This is much more accurate in terms of previous value. Also, we can clearly see that the free
electron fraction is not zero. Rather the electron density freezes out at a non-zero value of
about 2 × 10−4 .
Using the peebles equation solution as the input in Cr ’s expression, we can find the Peebles
Factor. It is shown in fig-9.

Figure 9: Peebles Factor as a function of redshift.

5 Conclusion
When Recombination ocuurs is a very important information we need for further study of
CMBR. As fluctions depend on when they enter the horizon.
In our work as we can see we have found the z of recombination using a simple equilibrium
model. Although it is simple, It is far off. Only using Non-equilibrium model we get the correct
value.

We have used Boltzmann equation for these non-equilibrium study. Although it is not
enough, It does a nicer job in the subject. Finally we have used Effective Three-level
system to study recombination and it’s possible process. That gives us a result which is very
close to the observe value.
In future we would like to study the fluctuation of CMB and the process of reionization to
understand CMB more.

16
 6 Appendix
6.1 Hubble’s Constant for Matter-Radiation Equality
We know,
 H 2  Ω  
r
= Ωr a−4 + Ωm a−3 = Ωm a−3 1 + a−1 (55)
H0 Ωm
where Ωm and Ωr are the critical densities for matter and radiation.
At equality, a = aeq and ΩΩmr = aeq hence,
p
−3/2
 aeq 1/2
H =H0 Ωm a 1+
a
p  1 + z 1/2 (56)
H =H0 Ωm (1 + z)3/2 1 +
1 + zeq
Here we have used a ∝ (1 + z)−1 .

6.2 Codes
For eq n −18 and 19, we have

Figure 10: Pressure in a weakly interacting gas of particles.

For Plotting and solving saha, boltzmann and peeble equations I have used Julia Language.

1 using DifferentialEquations
2 using Plots, LaTeXStrings
3 using Interpolations, Roots
4

5 zeta_3 = 1.20205690316
6 #All in eV
7 EI = 13.6 #Hydrogen Binding,

17
8 Me = 5.10999e5 #Mass of electron
9 zeq = 3395 # Redshift for matter-radiation equality
10 H0 = 1.5e-33
11 fs_c = 1/137 #Fine structure Constant
12 m = 0.3089 #Matter critical density",
13 2s = 8.227*6.58e-16 # Two photon decay rate",
14 z0 = 2.725/(1.160451812e4)
15 = 6e-10 #baryon to photon ratio
16 nb = 0.76*(2**zeta_3)/^2#baryon [Link]
17

18 T(z) = BigFloat(z0*(1 + z))#Conversion of z to T in eV

19

20 #For Equilibrium Part "Saha Equation"

21 pre_fac_saha = nb*(2*/Me)^1.5
22 function Saha_Xe(z)
23 f = pre_fac_saha*T(z)^(3/2)*BigFloat(exp(EI/T(z)))
24 return (-1 + sqrt(1 + 4*f))/(2*f)
25 end
26 #---------------------------------
27 #Finding Where Xe = 0.5
28 soln_saha(z) = Saha_Xe(z) - 0.5
29 sol1 = find_zero(soln_saha, (1300,1400))#Z value where Xe = 0.5 for Saha
30 T(sol1)
31 #---------------------------------
32

33 # For Non-Equilibrium Part

34 function ab_beta(z)#Defining , b and
35 t = T(z); x = EI/t
36 alpha = BigFloat(9.8*(fs_c/Me)^2*sqrt(x)*log(x))
37 b = alpha*(Me*t/(2*))^(1.5)*BigFloat(exp(-x))
38 beta = b*exp(0.75*x)
39 return alpha, b, beta
40 end
41

42 H(z) = H0*sqrt(m)*(1 + z)^(1.5)*sqrt(1 + (1+z)/(1 + zeq))

43

44 # For Boltzmann without peeble factor

45 function boltz_diff(Xe, p, z)
46 prefac = H(z)*(1 + z)
47 alpha, b, beta = ab_beta(z)
48 fact_main = (alpha*nb*Xe^2 *T(z)^3) - b*(1-Xe)
49 return fact_main/prefac
50 end
51 Xe0 = BigFloat(1)#Initial Value of Xe
52 zspan = (BigFloat(1800), BigFloat(0))#Range of soln
53 prob_boltz = ODEProblem(boltz_diff, Xe0, zspan)#Define the problem
54 sol_boltz = solve(prob_boltz, alg_hints = [:stiff])#Solution
55 z_bol = sort(sol_boltz.t)#solution Z values
56 Xe_bol = sort(sol_boltz.u)#solution Xe values
57 #-----------------------------------------------------------------------------------
58 inter_boltz = linear_interpolation(z_bol,Xe_bol)#Linear Interpolation
59 bolt_05(z) = inter_boltz(z) - 0.5 #For our simplification

18
60 sol1 = find_zero(bolt_05, (1100,1500))#Z value at which Xe 0.5 for Boltzmann
61 #-----------------------------------------------------------------------------------
62 T(1368)
63 #Defining peeble factor
64 function Cr_p(z,u)
65 hh = H(z); temp = T(z)
66 alpha, b, beta = ab_beta(z)
67 2p_nu = BigFloat(27*hh*EI^3)
68 2p_di_with_out_xe = BigFloat(64*^2*nb*temp^3)*(1-u)
69 2p = BigFloat((2p_nu/2p_di_with_out_xe))
70 Cr = BigFloat((2s + 2p)/(2s + 2p + beta))
71 return Cr
72 end
73

74 #DifferentialEquation with peeble factor

75 function peeble_diff(Xe, p, z)
76 prefac = H(z)*(1 + z)
77 alpha, b, beta = ab_beta(z)
78 fact_main = (alpha*nb*Xe^2 *T(z)^3) - b*(1-Xe)
79 return Cr_p(z,Xe)*fact_main/prefac
80 end
81 Xe0 = 0.9999# Should be 1 but due to discontinuity 0.9999 is used
82 prob_peeble = ODEProblem(peeble_diff, Xe0, zspan)
83 sol_peeble = solve(prob_peeble, alg_hints = [:stiff])#solution peeble Eqn
84 z_peb = sort(sol_peeble.t)#z values for peeble soln
85 Xe_peb = sort(sol_peeble.u)#Xe values for peeble soln
86 cr_array = broadcast(Cr_p,z_peb,Xe_pebb)#Cr values
87 # Finding z for which Xe = 0.5 for peeble eqn
88 inter_peeb = linear_interpolation(z_peb,Xe_pebb)
89 pebb_05(z) = inter_peeb(z) - 0.5
90 sol3 = find_zero(pebb_05, (1100,1300))#z for which Xe = 0.5
91 T(1270)
92 #---------------------------------------------------------------------------------
93

94 #Plotting Xe for all cases

95 zs = range(0, 1800, length=10_000)
96 saha = Saha_Xe.(zs)
97 plot(yscale=:log10, minorgrid=true)
98 xlims!(0, 1800)
99 ylims!(1e-4, 1.5e+0)
100 xlabel!(L"Redshift ($z$)")
101 ylabel!(L"Free Electron Fraction ($X_e$)")
102 plot!([0,1800],[0.5,0.5],lc=:green,lw=2.5,label="",ls=:dash)#0.5 line
103 annotate!(700, 0.3, L"$X_e = 0.5$")
104 #Saha eqn plot
105 plot!(zs, saha, label="",lc=:black, lw=2)
106 annotate!(1100, 2*1e-4, "Saha", :black)
107 #Boltzmann Eqn plot
108 plot!(z_bol, Xe_bol, label="",lc=:red, lw=2)
109 annotate!(300, 2*1e-4, "Boltzmann", :red)
110 #Peeble eqn plot
111 plot!(z_peb, Xe_peb, label="",lc=:blue, lw=2)

19
112 annotate!(250, 9*1e-4, "Peebles", :blue)
113 # Plotting Cr
114 plot(yscale=:log10,minorgrid=true)
115 xlims!(0, 1800)
116 ylims!(1e-4, 1.5e+0)
117 xlabel!(L"Redshift ($z$)")
118 ylabel!(L"Peeble Constant ($C_r$)")
119 plot!(z_peb,cr_array, minorgrid=true, label="", lw=2.2)
120 #Plotting the line for min Cr
121 c_min = minimum(cr_array)
122 plot!([0,1800],[c_min, c_min],label="",lw=2, ls=:dash)
123 annotate!(1000, c_min - c_min/3, L"$C_r = 1.827\times 10^{-3}$", :black)

I am not sure why Γ symbol is not printing in the document in latex (in code). But I guess It is
understandable where the error is. Here is the github link - [Link]
Online_codes/blob/main/[Link]

References
[1] Scott Dodelson. Modern Cosmology, Academic Press Inc; Second Edition, 2003.
[2] Daniel Baumann Cosmology, Third Edition, Cambridge University Press, 2022.
[3] Barbara Ryden. Introduction to Cosmology, Second Edition, Cambridge University Press,
2016.
[4] Peebles P.J.E. Principles of physical cosmology, Princeton, 1993.
[5] Steven Weinberg. Cosmology, Oxford University Press, 2019.
[6] P.J.E. Peebles. Recombination of the Primeval Plasma, The Astrophysical Journal, Vol.
153, July 1968.
[7] Arfken and Weber. Mathematical Methods for Physicists, 5th Edition, Academic Press,
2002.

20