  In the recent Letter [1] we considered the approach of nonequilibrium pair
plasma towards thermal equilibrium state adopting a kinetic treatment and
solving numerically the relativistic Boltzmann equations. It was shown that
plasma in the energy range 0.1-10 MeV first reaches kinetic equilibrium, on a
timescale t_{k}<10^{-14} sec, with detailed balance between binary interactions
such as Compton, Bhabha and Moller scattering, and pair production and
annihilation. Later the electron-positron-photon plasma approaches thermal
equilibrium on a timescale t_{th}<10^{-12} sec, with detailed balance for all
direct and inverse reactions. In the present paper we systematically present
details of the computational scheme used in [1], as well as generalize our
treatment, considering proton loading of the pair plasma. When proton loading
is large, protons thermalize first by proton-proton scattering, and then with
the electron-positron-photon plasma by proton-electron scattering. In the
opposite case of small proton loading proton-electron scattering dominates over
proton-proton one. Thus in all cases the plasma, even with proton admixture,
reaches thermal equilibrium configuration on a timescale t_{th}<10^{-11} sec.
We show that it is crucial to account for not only binary but also triple
direct and inverse interactions between electrons, positrons, photons and
protons. Several explicit examples are given and the corresponding timescales
for reaching kinetic and thermal equilibria are determined.
