The phenomenon of Bose-Einstein condensation is traditionally associated with and experimentally verified for low temperatures: either of nano-Kelvin scale for alkali atoms [1-3] or room temperatures for quasi-particles [4,5] or photons in two dimensions [6]. Here we demonstrate out of first principles that for certain initial conditions non-equilibrium plasma at relativistic temperatures of billions of Kelvin undergoes condensation, predicted by Zeldovich and Levich in their seminal work [7]. We determine the necessary conditions for the onset of condensation and discuss the possibilities to observe such a phenomenon in laboratory and astrophysical conditions.

The non-equilibrium electron-positron-photon plasma thermalization process is studied using relativistic Boltzmann solver, taking into account quantum corrections both in non-relativistic and relativistic cases. Collision integrals are computed from exact QED matrix elements for all binary and triple interactions in the plasma. It is shown that in non-relativistic case (temperatures $k_B T\leq 0.3 m_e c^2$) binary interaction rates dominate over triple ones, resulting in establishment of the kinetic equilibrium prior to final relaxation towards the thermal equilibrium, in agreement with the previous studies. On the contrary, in relativistic case (final temperatures $k_B T\geq 0.3 m_e c^2$) triple interaction rates are fast enough to prevent the establishment of kinetic equilibrium. It is shown that thermalization process strongly depends on quantum degeneracy in initial state, but does not depend on plasma composition.

We present a new efficient method to compute Uehling-Uhlenbeck collision integral for all two-particle interactions in relativistic plasma with drastic improvement in computation time with respect to existing methods. Plasma is assumed isotropic in momentum space. The set of reactions consists of: Moeller and Bhabha scattering, Compton scattering, two-photon pair annihilation, and two-photon pair production, which are described by QED matrix elements. In our method exact energy and particle number conservation laws are satisfied. Reaction rates are compared, where possible, with the corresponding analytical expressions and convergence of numerical rates is demonstrated.