text:
# A qualitative and numerical study of the matter-radiation interaction in Kantowski-Sachs cosmologies

## I Introduction

The radiative era of cosmic evolution extends from the end of cosmological nucleosynthesis to the decoupling of baryonic matter and radiation, covering the temperature range $`4\times 10^310^6`$ K (roughly between 1 eV and 1 keV). During this period cosmic matter can be describedkotu ; padma ; peacock  as a mixture of two main non–interacting components: one, a non-relativistic and collisionless gas of WIMP’s (cold dark matter, CDM), the other a tightly coupled mixture of non-relativistic baryons, electrons and ultra-relativistic matter (“radiation”, i. e., photons and neutrinos). The standard approach to the radiative era consists of using a Friedmann-Lemaître-Robertson-Walker (FLRW) space-time background kotu ; padma ; peacock  whose sources are described either by equilibrium kinetic theory bern , gauge invariant perturbations efsta , or by hydrodynamical models ct ; rdm ; pm  which in general, fail to incorporate a physically plausible description of the matter-radiation interaction since they assume a full thermodynamical equilibrium throughout the evolution. Since the tight coupling between electrons, baryonic matter and radiation follows from various processes of radiative interaction wei ; sw ; peeb , mostly involving photons and electrons, we can ignore the non-interacting neutrinos and assume the baryon–electron–photon mixture to evolve with a common temperature (local thermal equilibrium) and to behave as a single fluid, the “radiative fluid”. This fluid must be dissipative in order to provide an adequate macroscopic model for these interactions wei . Ideally, all dissipative fluxes (heat flux, bulk and shear viscosities) should be taken into consideration in the study of this tight coupling. However, in order to deal with a mathematically tractable problem, while still aiming at a physically interesting generalization of previous work, we shall study the case in which only shear viscosity is present. Bulk viscosity is not significant in the temperature ranges we are considering wei , and, although neglecting the contribution of heat flux carries physical limitations, this is compensated by the ensuing mathematical simplification of the field and transport equations. This approach has already been tested in various known and new exact solutions susstr .

The simplest class of metrics allowing for anisotropic shear viscous stresses are the Kantowski-Sachs cosmologies ksmh ; krs , characterized by a 4-dimensional isometric group. As the source of space-time we consider a momentum-energy tensor made of CDM (the WIMP gas) and the dissipative radiative fluid whose anisotropic pressure can be identified with shear viscous pressure. Considering the WIMP particle to be the supersymmetric lightest neutralino, with mass $`100`$ GeV Ellis ; Report ; Torrente , we can safely assume that throughout the radiative era these WIMP’s are non-relativistic, collisionless and only interact gravitationally with radiation and baryonic matter. It is also reasonable to assume, for the prevailing temperatures of this era, that the pressure of the WIMP gas and the internal energy of the baryons and electrons are negligible in comparison with the radiation equilibrium pressure. Although the radiative era is dominated by radiation, the rest mass-energy density of the WIMP’s is not negligible and dominates that of the baryons and electrons, hence the full source, CDM plus radiative fluid, can be well approximated by a momentum energy tensor in which CDM provides the bulk of the rest mass energy (“matter”, which ends up dominating the whole dynamics), while the photons (“radiation”) provide the bulk of thermal and dissipative effects. The shear viscosity associated with this source must satisfy appropriate constitutive and transport equations from irreversible thermodynamics that comply with causality and stability hl1 ; hl2 ; wi ; wjs ; these thermodynamical theories are known generically as Extended Irreversible Thermodynamics (EIT) ddj ; jcl ; hl3 . The application of such theories to particular physical systems requires phenomenological coefficients, like the coefficient of shear viscosity, to be provided by Kinetic Theory. In particular, for the tight coupling of electrons, baryonic matter and radiation and its associated photon-electron interaction, the coefficients corresponding to the “radiative gas” model wei ; jcl ; uwi ; dd1 ; dd2  should be employed. The entropy production must be positive definite and the relaxation time of shear viscosity must be a positive and monotonously increasing function, somehow related to the collision times of the radiative processes associated with the radiative era. All these timescales must overtake the Hubble expansion time as baryonic matter and radiation decouple.

The paper is organized as follows: Sections II to IV present and discuss the field equations of Kantowski-Sachs geometry for a mixture of CDM and a radiative fluid, the application of Extended Irreversible Thermodynamics and the appropriate set of equations of state for the models, as well as the evolution equations for the geometric and state variables. The dynamical analysis is carried on in Section V by defining a set of normalized variables, which then leads to a self-consistent and well-behaved autonomous system of ordinary differential (governing evolution) equations. From the qualitative analysis, we identify a saddle point associated with a radiation dominated FLRW cosmology and contained in the invariant set $`\chi =0,`$ associated with the flat Bianchi I model. We argue that initial conditions must be defined near this point. In Section VI we discuss various assumptions on the form of the relaxation time for shear viscosity, these assumptions lead to the identification of a “near-Eckart” and transient regimes, respectively associated with a swift and slow rate of transiency. The effects of using a “truncated” transport equation are discussed qualitatively in Section VII, while Section VIII deals with the numerical analysis of the models bearing in mind the qualitative results obtained in previous sections. The main result that follows from the qualitative and numerical analysis is the fact that a physically plausible evolution is possible only for: (a) initial conditions and evolution close to $`\chi =0;`$ and (b) using the full (not truncated) shear viscosity transport equation of EIT. A detailed discussion and summary of these results is provided in Section IX.

## II Kantowski-Sachs cosmologies with anisotropic stresses

The simplest non-FLRW cosmological metric allowing for anisotropic pressure is that of Kantowski-Sachs (KS) models:

$$ds^2=c^2dt^2+A^2(t)dr^2+B^2(t)\left[d\theta ^2+\mathrm{sin}^2(\theta )d\varphi ^2\right].$$
(1)
For a co-moving 4-velocity $`u^a,`$ the expansion scalar and shear tensor associated with the metric (1) are:

$`\mathrm{\Theta }={\displaystyle \frac{\dot{A}}{A}}+{\displaystyle \frac{2\dot{B}}{B}},\sigma _b^a=\text{diag}[0,2\sigma ,\sigma ,\sigma ],`$
$`\sigma {\displaystyle \frac{1}{3}}\left({\displaystyle \frac{\dot{B}}{B}}{\displaystyle \frac{\dot{A}}{A}}\right),`$ (2)
where a dot denotes the derivative with respect to proper time of fundamental observers, which for the KS metric (1) in co-moving coordinates is given by $`t,`$ (i. e., $`\dot{A}=A_{,t}=u^aA_{,a}`$). We consider as the source of (1), the following stress-energy tensor:

$$T^{ab}=\rho u^au^b+ph^{ab}+\mathrm{\Pi }^{ab},$$
(3)
where $`h^{ab}=c^2u^au^b+g^{ab}`$ and $`\mathrm{\Pi }_b^a`$ is the anisotropic pressure tensor satisfying $`\mathrm{\Pi }_{ab}u^b=\mathrm{\Pi }_a^a=0.`$ The most general form of this tensor for the metric (1) is:

$$\mathrm{\Pi }_b^a=\text{diag}[0,2P,P,P],$$
(4)
where $`P=P(t)`$ is an arbitrary function to be determined by the field equations and subjected to an evolution law for a given physical model associated with (1) and (3). The field equations then become:

$`\kappa \rho =G_t^t={\displaystyle \frac{\dot{B}^2}{B^2}}+{\displaystyle \frac{2\dot{B}}{B}}{\displaystyle \frac{\dot{A}}{A}}+{\displaystyle \frac{1}{B^2}},`$ (5)
$`3\kappa p=2G_\theta ^\theta +G_r^r=`$
$`{\displaystyle \frac{\dot{B}^2}{B^2}}{\displaystyle \frac{2\ddot{A}}{A}}{\displaystyle \frac{4\ddot{B}}{B}}{\displaystyle \frac{2\dot{B}}{B}}{\displaystyle \frac{\dot{A}}{A}}{\displaystyle \frac{1}{B^2}},`$ (6)
$`3\kappa P=G_\theta ^\theta G_r^r=`$
$`{\displaystyle \frac{\dot{B}^2}{B^2}}{\displaystyle \frac{\ddot{A}}{A}}+{\displaystyle \frac{\ddot{B}}{B}}{\displaystyle \frac{\dot{B}}{B}}{\displaystyle \frac{\dot{A}}{A}}+{\displaystyle \frac{1}{B^2}},`$ (7)
where $`\kappa =8\pi G/c^2,`$ while the energy balance is given by:

$$\dot{\rho }+(\rho +p)\mathrm{\Theta }+6\sigma P=0.$$
(8)
### II.1 Mixture of cold dark matter and a radiative fluid

We will assume that the stress-energy tensor (3) corresponds to a mixture of a non-relativistic gas of WIMP’s and a radiative fluid with shear viscosity corresponding to a “tightly coupled” mixture of photons, electrons and baryons sharing a common temperature $`T.`$ Hence, $`\rho `$ and $`p`$ in eqn. (3) are the total mass-energy density and equilibrium pressure given by:

$`\rho =\rho _\text{w}+\rho _\text{b}+\rho _\text{e}+\rho _\gamma ,`$
$`p=p_\text{w}+p_\text{b}+p_\text{e}+p_\gamma ,`$ (9)
$`p_\gamma =\rho _\gamma /3,`$
with the subindices $`\gamma ,w,b,`$ and $`e`$ denoting photons, WIMP’s, baryons and electrons, respectively. The three latter components satisfy each the equation of state of a non-relativistic ideal gas:

$`\rho _\text{w}=(m_\text{w}c^2+3k__BT_\text{w}/2)n_\text{w},p_\text{w}=n_\text{w}k__BT_\text{w},`$
$`\rho _\text{b}=(m_\text{b}c^2+3k__BT/2)n_\text{b},p_\text{b}=n_\text{b}k__BT,`$ (10)
$`\rho _\text{e}=(m_\text{e}c^2+3k__BT/2)n_\text{e},p_\text{e}=n_\text{e}k__BT,`$
where $`m_\text{w},m_\text{b},m_\text{e}`$ are the respective particle masses of the WIMP’s, the baryons (a proton mass) and the electrons, $`k__B`$ is Boltzmann’s constant, $`T_\text{w}`$is the temperature of the WIMP gas and $`T`$ is the common temperature of the radiative mixture. During the radiative era creation/annihilation processes cease to be significant and so the particle number densities, $`n_\text{w},n_\text{b},n_\text{e},`$ satisfy conservation laws of the form:

$$\dot{n}+n\mathrm{\Theta }=0,\text{with}n=n_\text{w},n_\text{b},n_\text{e}$$
(11)
which can be integrated, leading to

$$n=\frac{N}{AB^2},\text{with}N=N_\text{w},N_\text{b},N_\text{e}$$
(12)
where $`N_\text{w},N_\text{b},N_\text{e}`$ are the constant (i. e., conserved) number of WIMP’s, baryons and electrons, respectively.

The radiation component of the radiative fluid can be given either in terms of:

($`i`$) the Stefan-Boltzmann law:

$$\rho _\gamma ^{_{sb}}=aT^4,p_\gamma ^{_{sb}}=aT^4/3,$$
(13)
where $`a\pi ^2k__B^4/(15\mathrm{}^3c^3)`$ is Stefan-Boltzmann constant, or

($`ii`$) an ultra-relativistic ideal gas:

$$\rho _\gamma ^{_{ig}}=3n_\gamma k__BT,p_\gamma ^{_{ig}}=n_\gamma k__BT,$$
(14)
where $`n_\gamma `$ is the number density of ultra-relativistic particles, subjected to a balance law analogous to (11) and given by an expression similar to (12) with the conserved photon number $`N_\gamma `$.

In order to simplify eqn. (9), we can examine the ratios of particle numbers and rest mass densities of the different particle components. Considering the currently estimated padma  ratio of photons per baryon, we have:

$$\nu _\text{b}\frac{N_\text{b}}{N_\gamma }=\frac{n_\text{b}}{n_\gamma }2.67\times 10^8\mathrm{\Omega }_\text{b}h^2,$$
(15)
where $`\mathrm{\Omega }_\text{b}`$ is the baryon abundance today (roughly $`0.04\pm 0.01`$) and $`h0.7`$ is the adimensional Hubble factor lan . Regarding the WIMP gas, if we assume that it is made up of the lightest supersymmetric neutralinos with $`m_\text{w}100`$ GeV Ellis ; Report ; Torrente , we have that Report :

$$\nu _\text{w}\frac{N_\text{w}}{N_\gamma }=\frac{n_\text{w}}{n_\gamma }2.82\times 10^8\mathrm{\Omega }_\text{w}h^2.$$
(16)
where $`\mathrm{\Omega }_\text{w}0.3\pm 0.1`$ is the neutralino (CDM) abundance today. Using eqns. (10), (12) and (14), we can rewrite (9) as:

$`\rho =m_\text{w}c^2n_\text{w}\left(1+{\displaystyle \frac{m_\text{b}\nu _\text{b}}{m_\text{w}\nu _\text{w}}}+{\displaystyle \frac{m_\text{e}\nu _\text{e}}{m_\text{w}\nu _\text{w}}}\right)+`$
$`3n_\gamma k__BT\left(1+{\displaystyle \frac{\nu _\text{b}}{2}}+{\displaystyle \frac{\nu _\text{e}}{2}}+{\displaystyle \frac{\nu _\text{w}T_\text{w}}{2T}}\right),`$
$`p=n_\gamma k__BT\left(1+\nu _\text{b}+\nu _\text{e}+{\displaystyle \frac{\nu _\text{w}T_\text{w}}{T}}\right),`$ (17)
Bearing in mind that for electrons $`N_\text{e}N_\text{b},`$ so that $`\nu _\text{e}\nu _\text{b},`$ while $`m_\text{e}m_\text{b}1`$ GeV, $`\nu _\text{w}1,`$ and $`\nu _\text{e}\nu _\text{b}1,`$ then,

$`{\displaystyle \frac{m_\text{b}\nu _\text{b}}{m_\text{w}\nu _\text{w}}}10^2{\displaystyle \frac{\mathrm{\Omega }_\text{b}}{\mathrm{\Omega }_\text{w}}}1,{\displaystyle \frac{m_\text{e}\nu _\text{e}}{m_\text{w}\nu _\text{w}}}10^5{\displaystyle \frac{\mathrm{\Omega }_\text{b}}{\mathrm{\Omega }_\text{w}}}1,`$
while for the temperature range $`4\times 10^3`$ K $`<T<10^6`$ K, we have:

$$0.013\frac{m_\text{w}c^2n_\text{w}}{n_\gamma k__BT}3.25,$$
(18)
showing that the radiative era is initially radiation dominated but rest mass energy density is not negligible and ends up becoming dominant. Therefore, even if $`T_\text{w}/T`$ in eqn. (17) is not negligible, we have:

$`\rho m_\text{w}c^2n_\text{w}+3n_\gamma k__BT,pn_\gamma k__BT,`$ (19)
The same type of approximation can be obtained if we use the Stefan-Boltzmann law (13), since the ratio of pressures in eqns. (10) and (13): $`p_\gamma ^{_{sb}}/p_\text{b}=aT^3/[3n_\text{b}k__B],`$ is proportional to $`\nu _\text{b}`$ (likewise for WIMP’s). Therefore, the mixture of a gas of WIMP’s and a radiative fluid can be accurately described in the desired temperature range by the approximated equation of state:

$$\rho =mc^2n^{(m)}+\rho ^{(r)},p=p^{(r)}=\rho ^{(r)}/3$$

with

$$mm_\text{w},n^{(m)}n_\text{w},p^{(r)}p_\gamma ,$$
(20)
where $`\rho ^{(r)}`$ follows from either one of eqns. (13) or (14), hence $`n^{(r)}=n_\gamma `$ and $`\rho ^{(m)}=m_\text{w}c^2n_\text{w}.`$ For the remaining of this paper, the superindices $`^{(r)}`$ and $`^{(m)}`$ will refer to quantities associated with photons (“radiation”) and WIMP’s (“matter”), respectively. From eqns. (5) and (6), the equation of state (20) can be given as the following constraint:

$$\frac{\ddot{A}}{A}+\frac{2\ddot{B}}{B}+\frac{2\dot{A}}{A}\frac{\dot{B}}{B}+\frac{\dot{B}^2}{B^2}+\frac{1}{B^2}\frac{\kappa mc^2n^{(m)}}{2}=0.$$
(21)
## III Extended Irreversible Thermodynamics

If the source (3) is meant to describe a mixture of non-relativistic CDM and a radiative fluid (as argued in previous sections), the anisotropic pressure must be identified with a shear viscous stress of the latter fluid and must be compatible with a suitable thermodynamical formalism. We shall consider the so called “Extended Irreversible Thermodynamics” (EIT) wi ; wjs ; ddj ; jcl , a theory complying with causality and stability requirements hl3  and supported by the Kinetic Theory of gases, Information Theory and by the Theory of Hydrodynamical Fluctuations jcl . When shear viscosity is the only dissipative agent, the corresponding generalized entropy current, $`S^a,`$ obeying the usual balance law with non-negative divergence, and up to second order in $`\mathrm{\Pi }^{ab},`$ takes the form:

$$S^a=nSu^a,S=S_{(e)}\frac{\alpha \mathrm{\Pi }_{cd}\mathrm{\Pi }^{cd}}{2nT},$$
(22)
where $`S`$ is the entropy per particle, $`n=n^{(r)}+n^{(m)}=(1+\nu _\text{w})n^{(r)}`$ is the total particle number density ($`n^{(r)}`$), $`\alpha `$ is a phenomenological coefficient to be specified later and $`S_{(e)}`$ is defined by the equilibrium Gibbs equation:

$$nT\dot{S}_{(e)}=\dot{\rho }(\rho +p^{(r)})\frac{\dot{n}}{n}=3\dot{p}+4p\mathrm{\Theta }=\sigma _{ab}\mathrm{\Pi }^{ab},$$
(23)
where we have used eqns. (11) and (20) to eliminate $`\dot{\rho }.`$ Fulfillment of the second law of thermodynamics requires $`S_{;a}^a0,`$ which from the definition of $`S^a`$ and $`S`$ in (22) leads to:

$$\dot{S}0,$$
(24)
together with

$$\alpha =\frac{\tau _{\text{rel}}}{2\eta },$$
(25)
and the evolution equation of the viscous pressure, i. e., the transport equation hl3 ; matr :

$$\tau _{\text{rel}}\dot{\mathrm{\Pi }}_{cd}h_a^ch_b^d+\mathrm{\Pi }_{ab}\left[1+\frac{ϵ_0\eta T}{2}\left(\frac{\tau _{\text{rel}}}{T\eta }u^c\right)_{;c}\right]+2\eta \sigma _{ab}=$$

$$\tau _{\text{rel}}\dot{\mathrm{\Pi }}_{cd}h_a^ch_b^d+\mathrm{\Pi }_{ab}\left[1\frac{ϵ_0\tau _{\text{rel}}}{2}\left(\frac{\dot{T}}{T}+\frac{\dot{\eta }}{\eta }\frac{\dot{\tau }_{\text{rel}}}{\tau _{\text{rel}}}\mathrm{\Theta }\right)\right]+$$

$$2\eta \sigma _{ab}=0,$$
(26)
where $`\eta ,T,\tau _{\text{rel}},\sigma _{ab},\mathrm{\Pi }_{ab}`$ are the coefficient of shear viscosity, the temperature, the relaxation time, and the shear and shear viscosity tensors, respectively. The parameter $`ϵ_0`$ can take only two values: $`ϵ_0=1`$ (“full” transport equation) and $`ϵ_0=0`$ (“truncated” transport equation, also known as the Israel-Stewart equation wi ; wjs ; jcl ), while Eckart’s non-causal transport equation follows by setting $`\tau _{\text{rel}}=0.`$ The coefficient of shear viscosity as well as other related quantities can be obtained by a variety of means wei  including Kinetic Theory, Statistical Mechanics or both uwi ; dd1 ; dd2 ; zim . Unless specifically stated otherwise, we shall consider only the full transport equation $`ϵ_0=1.`$ We will discuss the implications of the truncated equation ($`ϵ_0=0`$) in Section VII. Evaluating $`\dot{S}`$ from (22) and using eqns. (23), (25), and (26) we obtain:

$$\dot{S}=\frac{\mathrm{\Pi }_{ab}\mathrm{\Pi }^{ab}}{2\eta nT}=\frac{2[S_{(e)}S]}{\tau _{\text{rel}}}0,$$
(27)
where we have used eqn. (22) and $`ϵ_0=1`$ in eqn. (26). Notice how $`\dot{S}`$ takes a very simple form, illustrating the role of $`\tau _{\text{rel}}`$ as the reference timescale associated with the entropy change from $`S`$ to $`S_{(e)}>S.`$ Also, the second law of thermodynamics, (i. e., (24)) is fulfilled if $`\eta 0`$ or, equivalently, if $`\tau _{\text{rel}}>0`$ and $`S_{(e)}S`$ hold. Another important requirement that follows from the second law of thermodynamics and the stability of equilibrium states is that $`S`$ be a convex functional, i. e., $`\delta ^2S<0.`$ For the KS models all quantities depend only on time and so a necessary (but not sufficient) condition is given by $`\ddot{S}<0,`$ which leads to:

$$\left(1+\frac{\dot{\tau }_{\text{rel}}}{2}\right)\mathrm{\Pi }_{ab}\mathrm{\Pi }^{ab}+2\eta \sigma _{ab}\mathrm{\Pi }^{ab}>0.$$
(28)
For the applicability of the general relations, eqns. (22), (24), (25), (26), (27) and (28) to the models we are concerned, we must impose the equation of state (20), with either one of the choices eqn. (13) or eqn. (14). As a physical reference to infer the form that the coefficients $`\tau _{\text{rel}},\eta ,\alpha `$ may take for the radiative fluid, consider the“radiative gas” model associated with the photon-electron interaction wei ; ddj ; jcl ; uwi , characterized by $`p,\rho `$ complying with the equations of state discussed in section III. For the radiative gas the forms of $`\eta ,\alpha `$ in terms of the relaxation time of the dissipative process, $`\tau _{\text{rel}}`$, are:

$$\eta _{_{(rg)}}=\frac{4}{5}p^{(r)}\tau _{\text{rel}},\alpha _{_{(rg)}}=\frac{5}{8p^{(r)}},$$
(29)
where $`p^{(r)}`$ is either $`aT^4/3`$ or $`n^{(r)}k__BT`$ (eqns. (13) or (14), respectively) and the subscript “$`(rg)`$” emphasizes that these quantities are specific to the radiative gas. Applying eqn. (29) into eqns. (22) and (23), we get for the entropy per particle:

$$S=S_{(e)}\frac{15P^2}{8p^{(r)}n^{(r)}T},$$
(30)
$$\dot{S}_{(e)}=\frac{3\sigma P}{n^{(r)}T},$$
(31)
where we have neglected the entropy of the non-relativistic matter, so that $`n^{(m)}+n^{(r)}=(1+\nu _\text{w})n^{(r)}`$ $`n^{(r)},`$ and $`S`$ is approximately the photon entropy (this is justified because the number of photons is so much larger than that of WIMP’s, baryons and electrons). The transport equation (26) becomes:

$$\dot{P}+\frac{8}{5}p^{(r)}\sigma +\left(\frac{4}{3}\mathrm{\Theta }+\frac{1}{\tau _{\text{rel}}}\right)P+\frac{\lambda \sigma }{p^{(r)}}P^2=0,$$
(32)
where $`\lambda 1+1/(2\lambda _0)`$ and $`\lambda _0=1/2`$ or $`2`$ (for $`p^{(r)}`$ given by the Stefan-Boltzmann law (13) or the ideal gas law (14), respectively); also condition (27) takes the form:

$$\dot{S}=\frac{15P^2}{4p^{(r)}n^{(r)}T\tau _{\text{rel}}}=\frac{2[S_{(e)}S]}{\tau _{\text{rel}}}0.$$
(33)
Notice that $`S<S_{(e)}`$ and $`\tau _{\text{rel}}>0`$ must hold in order for $`\dot{S}>0`$ to be satisfied, while (23) implies that we must have $`\sigma P<0`$ in order that $`\dot{S}_{(e)}>0.`$ Regarding the interpretation of $`S_{(e)},`$ if we assume $`\dot{S}_{(e)}>0`$<sup>1</sup><sup>1</sup>1Notice that $`\dot{S}_{(e)}>0`$ is a sufficient but not necessary condition for $`\dot{S}>0`$. Under the framework of Extended Irreversible Thermodynamics, the latter is the physicaly relevant condition. and the Steffan-Boltzmann law (13), eqn. (23) then yields:

$$S_{(e)}=\frac{4aT^3}{3n^{(r)}},$$
(34)
a function that is only constant (equilibrium photon entropy) if $`P=0.`$ This constant is given explicitly by the black body formulae lan :

$$n^{(r)}|_{_{P=0}}=\frac{30\zeta (3)aT^3}{k__B\pi ^4}$$

$$S_{(e)}|_{_{P=0}}=\frac{2\pi ^4k__B}{45\zeta (3)}3.60k__B.$$
(35)
where $`\zeta `$ is the Riemann zeta function. However, if we characterize the evolution to equilibrium as $`P0,`$ then, as this evolution proceeds, $`S_{(e)}`$ in eqn. (34) must tend to the constant entropy given in (35). Hence, we can identify eqn. (35) as the equilibrium state associated with eqn. (27) and eqn. (33), attained as $`P0`$ and $`SS_{(e)}`$ (i. e. as the radiation relaxes) in the timescale provided by the relaxation time $`\tau _{\text{rel}}.`$ For the ideal gas law, equation (23) does not yield eqn. (34), but $`S_{(e)}k__B\mathrm{ln}(T^3/n^{(r)}),`$ an expression that coincides with eqns. (34) and (35) only in equilibrium (if $`\dot{S}_{(e)}=0`$). However, since both EIT and Eckart’s theory assume near equilibrium states, quantities like $`P^2`$ and $`\sigma P,`$, appearing in eqns. (30) to (33) must be small, hence the ratio $`T^3/n^{(r)}`$ is nearly constant and so we can also assume that $`S_{(e)}`$ given by eqn. (34) is approximately valid for the ideal gas law.

Conditions (33) and $`\delta ^2S<0`$ associated with eqn. (28) must be satisfied by any self consistent thermodynamical system. The importance of these conditions will become evident when discussing the numerical integration of the evolution equations. The relaxation time, $`\tau _{\text{rel}},`$ is qualitatively analogous to and larger than the mean collision time between particles and it may, in principle, be estimated by collision integrals provided the interaction potential is known. Since we are concerned with mixtures of baryons, electrons and photons evolving in the temperature range $`4\times 10^3`$ K $`<T<10^6`$ K (from the end of cosmic nucleosynthesis to decoupling), convenient references for comparing $`\tau _{\text{rel}}`$ are the collision times associated with Compton and Thompson scatterings zp :

$$t_c=\frac{m_\text{e}c^2}{k__BT}t_\gamma ,$$
(36)
$$t_\gamma =\frac{1}{2c\sigma __Tn_\text{b}}\left[1+\left(1+\frac{4h^3n_\text{b}e^{B_0/k__BT}}{\left(2\pi m_\text{e}k__BT\right)^{3/2}}\right)^{1/2}\right],$$
(37)
where $`\sigma __T,B_0,m_\text{e},`$ and $`h`$ are the Thompson scattering cross section ($`6.65\times 10^{25}`$ cm<sup>2</sup>), the hydrogen atom binding energy ($`13.6`$ eV), the electron mass and Planck’s constant, respectively. Equation (37) is obtained from the number density of free electrons provided by the Saha equation. Notice that we are using the baryon number density, $`n_\text{b},`$ and not the number density of WIMPS, $`n^{(m)}`$. For higher temperatures in the range of interest, Compton scattering is the most efficient radiative process keeping baryonic matter and radiation tightly coupled, though it is no longer effective in lower and intermediate temperature ranges ($`T<10^4`$ K). The photon-electron interaction of the radiative era requires that microscopic collision times $`t_\gamma ,t_c,`$ as well as $`\tau _{\text{rel}},`$ be much smaller than the timescale of cosmic expansion given by the Hubble time, approximately $`t__H3/\mathrm{\Theta }.`$ For the lower temperature range of the radiative era, just before recombination, Thomson scattering becomes the dominant radiative process, so that the decoupling of baryonic matter and radiation can be associated with the condition $`t_\gamma =t__H,`$ which should be approximately equivalent to $`\tau _{\text{rel}}=t__H.`$

## IV Evolution Equations

Since we need to determine a self-consistent set of ordinary differential equations governing the KS models, it is convenient to express the field equations and eqn. (21) in terms of $`\mathrm{\Theta }`$ and $`\sigma `$ by eliminating $`\dot{A},\dot{B}`$ from (2):

$$\frac{\dot{A}}{A}=\frac{\mathrm{\Theta }}{3}2\sigma ,\frac{\dot{B}}{B}=\frac{\mathrm{\Theta }}{3}+\sigma ,$$
(38)
which leads to:

$`\kappa \rho ={\displaystyle \frac{\mathrm{\Theta }^2}{3}}3\sigma ^2+{\displaystyle \frac{1}{B^2}},`$ (39)
$`\kappa p^{(r)}={\displaystyle \frac{\mathrm{\Theta }^2}{3}}{\displaystyle \frac{2\dot{\mathrm{\Theta }}}{3}}3\sigma ^2{\displaystyle \frac{1}{3B^2}},`$ (40)
$`\kappa P=\dot{\sigma }+\sigma \mathrm{\Theta }+{\displaystyle \frac{1}{3B^2}},`$ (41)
while, using eqns. (39) and (40), the equation of state (21) becomes:

$$\dot{\mathrm{\Theta }}+\frac{2\mathrm{\Theta }^2}{3}+3\sigma ^2+\frac{1}{B^2}\frac{\kappa mc^2n^{(m)}}{2}=0.$$
(42)
We can eliminate $`B`$ from the equations above with the help of equation (39) (which becomes a constraint) and the second of equations (38) (the evolution equation for $`B`$). Using equation (20), equations (41) and (42) then become:

$$\dot{\sigma }+\sigma ^2+\sigma \mathrm{\Theta }\frac{\mathrm{\Theta }^2}{9}+\frac{\kappa }{3}mc^2n^{(m)}+\kappa p^{(r)}\kappa P=0,$$
(43)
$$\dot{\mathrm{\Theta }}+\frac{\mathrm{\Theta }^2}{3}+6\sigma ^2+\frac{\kappa }{2}mc^2n^{(m)}+3\kappa p^{(r)}=0,$$
(44)
which are the evolution equations for $`\sigma `$ and $`\mathrm{\Theta }.`$ Since we are assuming particle number conservation, an evolution equation for non-relativistic particle number density (the WIMP’s) follows from equation (11):

$$\dot{n}^{(m)}+n^{(m)}\mathrm{\Theta }=0,$$
(45)
a conservation law satisfied also by $`n^{(r)}=n^{(m)}/\nu _\text{w}`$ (if using the ideal gas law (14)). Another evolution equation is provided by (8), which applied to the equation of state (20) and using $`\sigma _{ab}\mathrm{\Pi }^{ab}=6\sigma P`$ yields:

$$\dot{p}^{(r)}+\frac{4}{3}p^{(r)}\mathrm{\Theta }+2\sigma P=0,$$
(46)
becoming the evolution equation for $`p^{(r)}`$. The transport equation (32) derived in the previous section, namely:

$$\dot{P}+\frac{8}{5}p^{(r)}\sigma +\left(\frac{4\mathrm{\Theta }}{3}+\frac{1}{\tau _{\text{rel}}}\right)P+\frac{\lambda \sigma }{p^{(r)}}P^2=0,$$
(47)
is the evolution equation for the shear stress $`P.`$ In addition to these evolution laws, eqn. (33) can be thought of as an evolution equation for $`S,`$ while we can transform eqn. (46) into an evolution equation for $`T`$ by using either (13) or (14), leading to:

$$\frac{\dot{T}}{T}+\frac{\mathrm{\Theta }}{3}+\frac{\sigma P}{\lambda _0p^{(r)}}=0.$$
(48)
Equations (43), (44), (45), (46), and (47) represent a self-consistent and closed system of first order ODE’s for $`n^{(m)},P,p^{(r)},\sigma ,\mathrm{\Theta }.`$ Notice that this set of evolution equations is fully determined if $`\tau _{\text{rel}}`$ is known. In the cosmological context $`\tau _{\text{rel}}`$ might be proportional to the timescale defined by the expansion scalar (approximately the Hubble time):

$$\tau _{\text{rel}}t__H3/\mathrm{\Theta }.$$
(49)
Alternatively, and depending on the temperature and energy range one is considering, $`\tau _{\text{rel}}`$ could be identified as proportional to a collision time (say, Thomson or Compton scattering) given by equations (36) or (37), i. e.,

$$\tau _{\text{rel}}t_\gamma (n,T),\tau _{\text{rel}}t_c(n,T),$$
(50)
Since we can assume in equation (47) two different equations of state for the radiation component, strictly speaking, equations (43)-(47) constitute two different systems of evolution equations parameterized by the two possible values of $`\lambda _0`$ and $`\lambda `$: –Stefan-Boltzmann law:

$$3p^{(r)}=aT^4,\text{with}\lambda _0=2\text{and}\lambda =5/4,$$

–Ideal gas law:

$$p^{(r)}=n^{(r)}k__BT,\text{with}\lambda _0=1/2\text{and}\lambda =2.$$
(51)
## V Dynamical Analysis

### V.1 The governing equations.

Let us define:

$$Q\frac{P}{p^{(r)}},$$
(52)
a ratio that is related (from eqn. (30)) to the deviation from equilibrium:

$$S_{(e)}S=\frac{15}{8}\lambda _0Q^2,$$

where:

$$\lambda _0=\frac{p^{(r)}}{n^{(r)}T}=\{\begin{array}{c}\frac{3}{4}S_{(e)},\text{Stefan-Boltzmann law}\\ \\ k__B,\text{Ideal Gas law}\end{array}$$
(53)
with $`S_{(e)}`$ given by eqn. (34). Therefore, $`Q`$ must be a small quantity for the cosmic times we are interested in. Note that when $`P`$ (and hence $`Q`$) is zero, the shear vanishes and the Kantowski-Sachs model reduces to an isotropic FLRW model. In principle, $`P`$ (and hence $`Q`$) can be positive or negative, but since $`P`$ represents a viscous pressure it should be negative in the expanding regime. We shall focus henceforth on the case in which $`Q`$ is negative. The energy conditions imply that $`1Q1/2.`$ In addition, on physical grounds we expect the second term in equation (8) to dominate the third term in this equation, which is satisfied whenever $`\sigma ^2Q^2/(\mathrm{\Theta }/3)^2<4.`$ However, since physically we expect $`Q^2`$ to be small, these constraints are satisfied handily.

We introduce now the new normalized variables:

$`\mathrm{\Sigma }{\displaystyle \frac{\sigma }{\mathrm{\Theta }/3}},\mathrm{\Omega }_{(m)}{\displaystyle \frac{\kappa mc^2n^{(m)}}{\mathrm{\Theta }^2/3}},`$
$`\mathrm{\Omega }_{(r)}{\displaystyle \frac{\kappa p^{(r)}}{(\mathrm{\Theta }/3)^2}},`$ (54)
a definition that is motivated by the fact that $`\mathrm{\Theta }/3`$ is approximately the Hubble expansion factor $`H=\mathrm{\Theta }/3+\sigma _{ab}n^an^b`$ for unit space-like vectors defined by $`n^an_a=1,u^an_a=0,`$ hence $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ are approximately equivalent to the observational parameters (the $`\mathrm{\Omega }`$’s) for the CDM and radiation contents of the mixture. We also define (for $`\mathrm{\Theta }>0`$) the new independent variable:

$$\frac{d}{dt}=\frac{\mathrm{\Theta }}{3}\frac{d}{d\tau }\tau =\frac{\mathrm{\Theta }}{3}𝑑t,$$
(55)
the evolution equations (43) to (47) become:

$$\mathrm{\Omega }_{(m)}^{}=\mathrm{\Omega }_{(m)}\left[14\mathrm{\Sigma }^2\mathrm{\Omega }_{(m)}2\mathrm{\Omega }_{(r)}\right],$$
(56)
$$Q^{}=\frac{8}{5}\mathrm{\Sigma }\frac{Q}{\tau _{\text{rel}}\mathrm{\Theta }/3}(\lambda 2)\mathrm{\Sigma }Q^2.$$
(57)
$$\mathrm{\Omega }_{(r)}^{}=\mathrm{\Omega }_{(r)}\left[2+2\mathrm{\Sigma }(Q2\mathrm{\Sigma })\mathrm{\Omega }_{(m)}2\mathrm{\Omega }_{(r)}\right],$$
(58)
$`\mathrm{\Sigma }^{}=\left(12\mathrm{\Sigma }\right)\left(1\mathrm{\Sigma }^2\right)\mathrm{\Omega }_{(m)}\left(2\mathrm{\Sigma }\right)/2`$
$`\mathrm{\Omega }_{(r)}\left(1Q\mathrm{\Sigma }\right),`$ (59)
where we have used

$$\frac{3\dot{\mathrm{\Theta }}}{\mathrm{\Theta }^2}=\frac{\mathrm{\Theta }^{}}{\mathrm{\Theta }}=12\mathrm{\Sigma }^2\frac{\mathrm{\Omega }_{(m)}}{2}\mathrm{\Omega }_{(r)},$$
(60)
and a prime denotes differentiation with respect to $`\tau .`$ Note that this last equation implies that $`\mathrm{\Theta }`$ is monotonically decreasing. The constraint (42) becomes:

$$\chi 1\mathrm{\Sigma }^2\mathrm{\Omega }_{(m)}\mathrm{\Omega }_{(r)}=\frac{3}{B^2\mathrm{\Theta }^2}.$$
(61)
We note that:

$`\chi ^{}=\chi \left[2\mathrm{\Sigma }(2\mathrm{\Sigma }1)+\mathrm{\Omega }_{(m)}+2\mathrm{\Omega }_{(r)}\right].`$
Clearly $`\chi =0`$ is an invariant set of the above differential equations, which corresponds to the Bianchi I (zero curvature) sub-case.

Eventually, the models re-collapse and $`\mathrm{\Theta }`$ changes sign. At the point of maximum expansion (when $`\mathrm{\Theta }=0`$) the variables above diverge and the normalized equations are no longer valid. However, for the times we are interested in, in the expanding phase far from re-collapse, the above variables and equations are valid. Indeed, in principle we can use the above system to follow the evolution of the models all the way back to the big bang. From equation (44) we can see that in this regime the curvature is small and that the variables $`\mathrm{\Sigma }^2,\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)}`$ are well-behaved. Compact variables can be defined by normalizing with $`\mathrm{\Theta }^2+B^2`$ (instead of $`\mathrm{\Theta }^2`$) that are valid for all times CG ; however, the physical assumptions used here are not valid at later times.

In addition to equations (56) to (59), evolution equations for $`n,T,S_{(e)}`$ and $`S`$ follow by using the variables defined in eqns. (54) and (55) in equations (45), (48), (31) and (33):

$`n^{}=3n,`$
$`n(\tau )=n_ie^{3\tau },\text{where:}n=n^{(m)}\text{or}n^{(r)},`$ (62)
$$T^{}=T\left[1+\frac{\mathrm{\Sigma }Q}{\lambda _0}\right],$$
(63)
$$S_{(e)}^{}=3\lambda _0\mathrm{\Sigma }Q,$$
(64)
$$S^{}=\frac{2t__H}{\tau _{\text{rel}}}\left[S_{(e)}S\right]=\frac{15t__H}{4\tau _{\text{rel}}}\lambda _0Q^2,$$
(65)
where $`\lambda _0`$ is given by eqn. (LABEL:Sdef4), the subindex $`i`$ denotes evaluation at an initial time $`\tau =\tau _i,`$ and $`t__H=3/\mathrm{\Theta },`$ is approximately the Hubble time which follows from eqn. (60) as:

$$t__H=\frac{3}{\mathrm{\Theta }_i}\mathrm{exp}\left\{\left[1+2\mathrm{\Sigma }^2+\frac{\mathrm{\Omega }_{(m)}}{2}+\mathrm{\Omega }_{(r)}\right]𝑑\tau \right\},$$
(66)
while the relation between physical time $`t`$ and $`\tau `$ follows from (55) and (60):

$$t=t__H𝑑\tau ,$$
(67)
where $`t__H`$ is given by (66) above.

### V.2 Qualitative Properties

Consider the dynamical implications of assuming that $`\tau _{\text{rel}}`$ is given by eqn. (49), namely:

$$\tau _{\text{rel}}=\frac{\gamma _0}{\mathrm{\Theta }/3}=\gamma _0t__H,$$
(68)
where $`\gamma _0>0`$ is a constant. For the range we are interested in, $`\gamma _01.`$ Equation (57) becomes:

$$Q^{}=\frac{8\mathrm{\Sigma }}{5}\frac{Q}{\gamma _0}(\lambda 2)\mathrm{\Sigma }Q^2,$$
(69)
and so equations (56), (58), (59) and (69) now constitute a closed four-dimensional system of first order autonomous differential equations for $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)`$, a system that depends on the value of the constant parameter $`\gamma _0.`$ Moreover, from equation (61) and the above discussion, in the regime we are interested in $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},`$ and $`\mathrm{\Sigma }`$ are bounded and physical conditions imply that $`1Q1/2.`$ Consequently a local analysis of the stability of the equilibrium points of this system will provide useful dynamical information.

Setting the right-hand-side of equation (69) to zero we obtain:

$$(\lambda 2)\mathrm{\Sigma }Q^2+\frac{Q}{\gamma _0}+\frac{8\mathrm{\Sigma }}{5}=0,$$
(70)
which is a quadratic equation for $`Q`$ if $`\mathrm{\Sigma }`$ is given and $`\lambda 2`$. Stefan-Boltzman law. ($`\lambda =5/4`$):

The equilibrium points at finite values are (note that all such points are given below; however, $`(\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},\mathrm{\Sigma }^2)`$ are not necessarily bounded by unity):

1. $`\mathrm{\Omega }_{(m)}=0=\mathrm{\Omega }_{(r)}`$, six points:

(ia) $`\mathrm{\Sigma }=1/2`$ and

$`Q_\pm =4[3\gamma _0\pm \sqrt{9\gamma _0^2+2/15}]/(9\gamma _0^2);`$

(ib) $`\mathrm{\Sigma }=1`$ and

$`Q_\pm =2[3\gamma _0\pm \sqrt{9\gamma _0^2+8/15}]/(9\gamma _0^2);`$

(ic) $`\mathrm{\Sigma }=1`$ and

$`Q_\pm =2[3\gamma _0\pm \sqrt{9\gamma _0^2+8/15}]/(9\gamma _0^2).`$
2. $`\mathrm{\Omega }_{(m)}=0,\mathrm{\Omega }_{(r)}0,`$ four points:

(iia) $`\mathrm{\Omega }_{(r)}=(12\gamma _0+17/15)/(9\gamma _0^2),`$ and

$`\mathrm{\Sigma }_\pm =Q_\pm =\pm 2\sqrt{3\gamma _0+8/15}/(3\gamma _0);`$

(iib) $`\mathrm{\Omega }_{(r)}{}_{\pm }{}^{}=\frac{(4+3\gamma _0)[2+\gamma _0(1\pm u)]}{\gamma _0[2+3\gamma _0(1u)]},`$

$`\mathrm{\Sigma }_\pm =(4+3\gamma _0)/[2+3\gamma _0(1u)],`$ and

$`Q_\pm =(2\pm 3\gamma _0u)/(4+3\gamma _0),`$ where

$`u=2\sqrt{1+(8\gamma _0/5)(4+3\gamma _0)}/(3\gamma _0);`$

$`\mathrm{\Sigma }=0=Q,\mathrm{\Omega }_{(r)}=1`$ is a particular solution.
3. $`\mathrm{\Omega }_{(r)}=0,\mathrm{\Omega }_{(m)}0,`$ one particular solution:

$`\mathrm{\Sigma }=0=Q,`$ and $`\mathrm{\Omega }_{(m)}=1.`$
4. $`\mathrm{\Omega }_{(m)}0,\mathrm{\Omega }_{(r)}0`$, two points:

$`\mathrm{\Omega }_{(m)}=(40+79/\gamma _0)/(64\gamma _0),`$

$`\mathrm{\Omega }_{(r)}=15(8+3\gamma _0)/(128\gamma _0),`$

$`\mathrm{\Sigma }=\pm \sqrt{5(8+3\gamma _0)/(2\gamma _0)}/8,`$ and

$`Q_\pm =4\sqrt{2(8+3\gamma _0)/(45\gamma _0)}.`$

The two equilibrium points given by (ia), can be shown to lead to a value for $`Q`$ which is unphysical in the sense that its magnitude is much too large (in fact, when $`\mathrm{\Omega }_{(m)}=0=\mathrm{\Omega }_{(r)},Q`$ would be expected to vanish). The equilibrium points (iv) are also unphysical: since clearly $`\mathrm{\Sigma }^21/4`$ and we get from (iv) that $`\mathrm{\Sigma }^2Q^2=1/4,`$ this would imply that $`Q^21`$ ! At the equilibrium points (ii) when $`\mathrm{\Omega }_{(r)}1`$ we have that $`\mathrm{\Sigma }=20Q/[\gamma _0(15Q^232)],`$ which leads to a quartic equation for Q; however, for physical values of the parameter $`\gamma _0,`$ this equation has no real roots and hence no solutions of physical interest. The two equilibrium points given by the particular solution of (ii) and by (iii), corresponding to FLRW models, namely $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)`$ given by $`(1,0,0,0)`$ and $`(0,1,0,0)`$, can easily be shown to be saddles. The equilibrium points (ib,c), namely $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(0,0,\pm 1,Q_\pm ),`$ which belong to the invariant set $`\chi =0`$ (i. e., correspond to a Bianchi I model with no curvature) have eigenvalues $`3,1/\gamma _0+3\lambda ^2\mathrm{\Sigma }Q/2,22\mathrm{\Sigma }Q,4\mathrm{\Sigma }.`$. In this case, the eigenvalue $`1/\gamma _0+3\lambda ^2\mathrm{\Sigma }Q/2`$ can only be positive if we take the ’positive square root’ of eqn. (37), (i. e., $`Q=Q_+`$). However, for $`\{\lambda ^2,\gamma _0\}1,`$ i. e., $`\lambda ^2\gamma _01,`$ it follows that the eigenvalue $`22\mathrm{\Sigma }Q_+`$ can never be positive. Consequently, this equilibrium point can not be a source.

We note that there are no sinks at finite values. However, this is to be expected since the models evolve toward maximum expansion at which the variables become unbounded. The models subsequently re-collapse.

Ideal Gas case ($`\lambda =2`$):

Equations (56), (58) and (59) remain unchanged and the governing equation (69) becomes:

$$Q^{}=\frac{8\mathrm{\Sigma }}{5}\frac{Q}{\gamma _0}.$$
(71)
At an equilibrium point we immediately see that:

$$Q=\frac{8\gamma _0\mathrm{\Sigma }}{5}.$$
(72)
However, there are no major qualitative changes in the analysis. In particular, there are fewer equilibrium points:

1. $`\mathrm{\Omega }_{(m)}=0=\mathrm{\Omega }_{(r)}`$, three points:

(ia) $`\mathrm{\Sigma }=1/2`$ and $`Q=4\gamma _0/5;`$

(ib) $`\mathrm{\Sigma }=1`$ and $`Q=8\gamma _0/5;`$

(ic) $`\mathrm{\Sigma }=1`$ and $`Q=8\gamma _0/5.`$
2. $`\mathrm{\Omega }_{(m)}=0,\mathrm{\Omega }_{(r)}0`$:

$`\mathrm{\Sigma }=0=Q,\mathrm{\Omega }_{(r)}=1`$ is a particular solution; otherwise,

A quadratic equation for $`\mathrm{\Sigma }`$ which gives another two equilibrium points.
3. $`\mathrm{\Omega }_{(r)}=0,\mathrm{\Omega }_{(m)}0`$, one particular solution:

$`\mathrm{\Sigma }=0=Q,`$ and $`\mathrm{\Omega }_{(m)}=1.`$
4. $`\mathrm{\Omega }_{(m)}0,\mathrm{\Omega }_{(r)}0`$, two points that imply:

$`\mathrm{\Sigma }^2=5/(16\gamma _0)5/16`$ (for $`\gamma _01`$), which leads to $`\mathrm{\Omega }_{(m)}+2\mathrm{\Omega }_{(r)}=15/(4\gamma _0)<0`$.

## VI Assumptions on the relaxation time

In the general dynamical system, eqns. (56) – (59), the relaxation time $`\tau _{\text{rel}}`$ needs to be specified in order for the governing system of equations to be closed. For the qualitative dynamical analysis we have assumed (68), so that $`1/(\tau _{\text{rel}}\mathrm{\Theta })`$ is a constant, since otherwise we would either not be dealing with an autonomous system, or would be looking at an autonomous but much more difficult dynamical system. However, equation (68) is a simplifying assumption that cannot be supported by thermodynamical arguments, perhaps a more realistic assumption would be to consider instead:

$$\tau _{\text{rel}}=\frac{\gamma (\tau )}{\mathrm{\Theta }/3}=\gamma (\tau )t__H$$
(73)
so that equations (57) and (65) become

$$Q^{}=\frac{8\mathrm{\Sigma }}{5}\frac{Q}{\gamma (\tau )}(\lambda 2)\mathrm{\Sigma }Q^2.$$
(74)
$$S^{}=\frac{2}{\gamma (\tau )}\left[S_{(e)}S\right],$$
(75)
where $`\gamma =\gamma (\tau )`$ is a function that could be suitably adjusted so that $`\tau _{\text{rel}}`$ has a form that is qualitatively analogous to that of microscopic timescales like eqn. (36) or (37), timescales that are physically relevant for the matter source under consideration. The ratio $`3/(\tau _{\text{rel}}\mathrm{\Theta })=t__H/\tau _{\text{rel}}=1/\gamma `$ should provide a comparison of the timescale for the relaxation (transient) effects in the radiative fluid with the timescale of cosmic expansion. Hence this ratio should approach unity as baryonic matter and radiation decouple, so that $`\gamma 1`$ should be a consistent choice for near decoupling conditions, while $`\gamma 1`$ or $`\gamma 1,`$ correspond to after decoupling (late times) and much before decoupling (earlier times). Ideally, we should obtain $`\tau _{\text{rel}}`$ from collision integrals associated with each of the various radiative processes occurring in the radiative era, but such an endeavor would merit a separate paper by itself and will not be attempted here. Instead, we will consider $`\tau _{\text{rel}}`$ as an “effective” relaxation time, encompassing the different radiative processes. We will examine the non-transient limit (or near-Eckart regime) and use the dynamical equations themselves in order to suggest a suitable form for $`\tau _{\text{rel}}`$. The discussion of this section will be complemented and tested numerically in section VIII (see also figures 5 and 6).

### VI.1 The “near-Eckart” and the transient regimes.

In order to examine the relaxation process as $`SS_{(e)},`$ we will assume that $`\gamma (\tau )`$ in (73) is a smooth function so that we can always expand it in the form $`\gamma \gamma (0)+\gamma ^{}(0)\tau +\gamma ^{\prime \prime }(0)\tau ^2/2.`$ Therefore, at early times $`\tau 0`$ we can always associate the constant $`\gamma _0`$ in eqn. (68) as $`\gamma _0=\gamma (0),`$ so that the corresponding form of $`\tau _{\text{rel}}`$ is approximately correct at least near $`\tau =0.`$ Also, at early times we must have $`\gamma (0)<1,`$ and so $`1/\gamma _01,`$ while $`\mathrm{\Sigma }(0)`$ and $`Q(0)`$ are necessarily small quantities. Hence the terms $`\mathrm{\Sigma }`$ and $`\mathrm{\Sigma }Q^2`$ in eqn. (74) will be much smaller than $`Q(0)/\gamma _0`$ and so we have that near $`\tau =0,Q^{}Q/\gamma _0,`$ so that:

$`QQ(0)\mathrm{exp}\left({\displaystyle \frac{\tau }{\gamma _0}}\right),`$
$`SS_{(e)}{\displaystyle \frac{15}{8}}\lambda _0Q^2(0)\mathrm{exp}\left({\displaystyle \frac{2\tau }{\gamma _0}}\right),`$ (76)
where $`\lambda _0`$ is given by eqn. (LABEL:Sdef4). Since the process of relaxation to equilibrium can be characterized as the decay of the dissipative flux $`Q0`$ as $`S`$ grows and asymptotically approaches the equilibrium state given by eqn. (35), the numerical value of $`\gamma _0`$ in (76) may be interpreted as a measure of the “rate of transiency” in terms of how fast or slow the system accomplishes this relaxation in comparison with the timescale provided by $`t__H.`$ We can then identify two possible situations

* The near-Eckart regime. If $`\gamma _01,`$ then eqn. (76) indicates a very fast relaxation with a very abrupt decay of $`Q`$ and $`S`$ to $`S_{(e)}.`$ In the very limit $`\gamma _00`$ (so that $`\tau _{\text{rel}}0`$ as well) we have $`QQ(0)\delta (\tau ),`$ so that the relaxation is infinitely fast, in agreement with the non-causal nature of Eckart’s classical theory. The fast relaxation associated with a very small $`\gamma _0`$ implies a very short duration of the relaxation process (small $`\tau _{\text{rel}}`$), indicating that the approximation $`\gamma =\gamma _0,`$ as well as the expressions in eqn. (76) for $`Q(\tau )`$ and $`S(\tau ),`$ are approximately valid for the whole evolution time. Hence, since $`Q`$ practically vanishes very quickly, we have for most of the evolution time that $`SS_{(e)},`$ agreeing with the fact that in Eckart’s theory the entropy is rigorously and unambiguously given by its equilibrium form (the hypothesis of “local equilibrium”, see section 1.3.1 of jcl ). The near-Eckart regime is appropriate to describe a given radiative process for which microscopic timescales are much smaller than the cosmological expansion timescale $`t__H.`$ See figures 5.
* The transient regime. If $`\gamma _0<1`$ but $`\gamma _0O(10^1)O(1),`$ then $`Q`$ and $`S`$ also decay but the relaxation process is much slower, hence the term “transient”. In this case, the expressions in eqn. (76) and the assumption $`\gamma =\gamma _0`$ are only good approximations for $`\gamma ,Q(\tau )`$ and $`S(\tau )`$ near $`\tau =0.`$ In general, we must use:

$`QQ(0)\mathrm{exp}\left[{\displaystyle \frac{d\tau }{\gamma (\tau )}}\right],`$
$`SS_{(e)}{\displaystyle \frac{15}{8}}\lambda _0Q^2(0)\mathrm{exp}\left[{\displaystyle \frac{2d\tau }{\gamma (\tau )}}\right],`$ (77)
leading to eqn. (75) and to:

$$S^{\prime \prime }S_{(e)}^{\prime \prime }\frac{2(1+\gamma ^{})(S_{(e)}S)}{\gamma ^2},$$
(78)
so that the fulfillment of $`S^{}>0`$ and $`S^{\prime \prime }<0`$ can be examined in terms of $`\gamma (\tau )`$ and $`S_{(e)}.`$ Sufficient (but not necessary) conditions follow by demanding that $`\gamma `$ is a monotonically increasing function ($`\gamma ^{}>0`$) and $`S_{(e)}^{}>0,S_{(e)}^{\prime \prime }<0.`$ Another condition on $`\gamma `$ is furnished by eqn. (28), which together with eqns. (29), (52), (54), (55), (60), and (73) yield

$`0<2+\gamma ^{}Q^2+`$
$`\gamma \left[{\displaystyle \frac{16}{5}}\mathrm{\Sigma }Q\left(1+2\mathrm{\Sigma }^2+{\displaystyle \frac{\mathrm{\Omega }_{(m)}}{2}}+\mathrm{\Omega }_{(r)}\right)Q^2\right],`$ (79)
a condition that should be tested numerically for any given choice of $`\tau _{\text{rel}}`$.

The relaxation time $`\tau _{\text{rel}}`$ can be approximated by a given $`\gamma =\gamma _0,`$ but as the decoupling era is reached, $`\gamma `$ must increase to $`O(1)`$ allowing for $`\tau _{\text{rel}}`$ to overtake $`t__H`$.

### VI.2 Dynamical relaxation times.

Let us consider the following ansatz:

$$\frac{1}{\gamma }=\frac{8}{5}\frac{\mathrm{\Sigma }}{Q}\left[1+\zeta (\tau )\right],$$
(80)
which provides an exact relation for the relaxation time in eqn. (73). In regimes in which the Eckart theory is a good approximation, we can assume $`\gamma \gamma _01,`$ so that eqn. (80) implies $`Q\mathrm{\Sigma }`$ but $`|\mathrm{\Sigma }/Q|1`$ and $`\zeta `$ is ’small’. Such regimes would correspond to a radiative process that takes place in timescales smaller than a mean collision time, thus decaying very fast to equilibrium. In case we wish to consider processes taking place on timescales comparable and larger than main collision times, then transient effects are important and the near-Eckart regime is no longer appropriate. In particular, we can construct an expression for $`\tau _{\text{rel}}`$ that acts as an ‘effective’ relaxation time that encompasses the relaxation times for the main radiative processes acting in the radiative era under consideration. Having this idea in mind, a reasonable and more general expression for $`\zeta `$ can be obtained from the dynamical equations themselves. Consider the conditions discussed above for a near-Eckart regime, and assume that:

$$Q=\mu _0\mathrm{\Sigma },$$
(81)
where $`\mu _0`$ is a non-negative constant. We obtain an expression analogous to eqn. (80) by substituting eqn, (81) into the (full EIT) evolution equations (57), and (59), leading to the consistency requirement:

$$\frac{1}{\gamma }=\frac{1}{\gamma _0}+a_0\mathrm{\Omega }_{(m)}b_0\mathrm{\Sigma }^2+\left(1\mu _0\frac{1}{\mathrm{\Sigma }}\right)\chi ,$$
(82)
where

$`{\displaystyle \frac{1}{\gamma _0}}=1+\mu _0+{\displaystyle \frac{8}{5\mu _0}}1+4\sqrt{{\displaystyle \frac{2}{5}}}>1,`$
$`a_0{\displaystyle \frac{1}{2}}\mu _0,b_0=1+(\lambda 3)\mu _0,`$ (83)
and $`\chi `$ is the curvature term $`\chi 1\mathrm{\Omega }_{(m)}\mathrm{\Omega }_{(r)}\mathrm{\Sigma }^2`$ given by eqn. (61), a term that is very small and can be neglected. This is an extremely simple dynamical relation for the relaxation time, and for a wide range of conditions it might be a very good approximation. In addition, it has some important physical properties. For early times in which $`\chi ,\chi /\mathrm{\Sigma },\mathrm{\Omega }_{(m)},\mathrm{\Sigma }^2`$ are very small, we have a near-Eckart regime associated with $`\gamma \gamma _01,`$ as expected. This is appropriate for the Compton scattering, the dominant radiative process in the early part of the period under consideration, a process that quickly thermalizes and ceases to be effective. At later times, as $`\mathrm{\Omega }_{(m)}`$ increases toward a value of order unity at recombination (and to a lesser extent, the curvature term also grows), the constants $`\gamma _0,a_0,b_0`$ in eqn. (83), as well as the initial conditions, can be selected in such a way that $`\gamma ,`$, given by (82), increases sufficiently as to allow $`\tau _{\text{rel}}`$ to overtake $`t__H`$ and to approach the characteristic timescale of Thomson scattering in eqn. (37). We show in section VIII that adequate parameters and initial conditions can be found so that $`\tau _{\text{rel}}`$ associated with eqn. (82) has the expected behavior (see figures 6).

## VII Truncated Theory

The discussion so far has been based on the full transport equation. In order to appreciate the effect of considering the truncated transport equation, it is useful to rewrite equation (26) with $`\eta `$ and $`\tau _{\text{rel}}`$ given by eqns. (29) and (73) in terms of $`Q`$ and $`\mathrm{\Sigma },`$ thus allowing the full and truncated equations to appear jointly. This yields:

$$Q^{}=\frac{8\mathrm{\Sigma }}{5}\left[\frac{1}{\gamma (\tau )}4(1ϵ_0)\right]Q(\lambda ϵ_02)\mathrm{\Sigma }Q^2.$$
(84)
where the full and truncated theories are respectively given by $`ϵ_0=1`$ and $`ϵ_0=0.`$ Assuming a transient regime, instead of eqn. (76) we have near $`\tau =0`$:

$`QQ(0)\mathrm{exp}\left[{\displaystyle \frac{\tau }{\gamma _0}}+4(1ϵ_0)\tau \right],`$
$`SS_{(e)}{\displaystyle \frac{15}{8}}\lambda _0Q^2(0)\mathrm{exp}\left\{2\left[{\displaystyle \frac{\tau }{\gamma _0}}+4(1ϵ_0)\tau \right]\right\},`$ (85)
so that using the truncated transport equation ($`ϵ_0=0`$) introduces a large linear term ($`4\tau `$) that is absent in the full theory. This linear term could change dramatically the form of $`Q`$ and the relaxation of $`S`$ to $`S_{(e)}.`$ As we show below, the truncated transport might lead to $`Q\mathrm{\Sigma }`$ being positive (implying that $`\dot{S}_{(e)}<0`$), so even at this level (early times) there might be problems of consistency in the truncated approach of EIT.

A comprehensive analysis in the case of the truncated theory, similar to that presented in Section V, can be undertaken. The evolution equations (56), (58), and (59) remain unchanged, while the evolution equation (57) for $`Q`$ must be replaced by the truncated equation that follows by setting $`ϵ_0=0`$ in eqn. (84)

$$Q^{}=\frac{8\mathrm{\Sigma }}{5}\left(\frac{1}{\gamma _0}4\right)Q+2\mathrm{\Sigma }Q^2.$$
(86)
At an equilibrium point:

$$2\mathrm{\Sigma }Q^2\left(\frac{1}{\gamma _0}4\right)Q\frac{8\mathrm{\Sigma }}{5}=0.$$
(87)
However, we immediately note from equation (87) that close to an equilibrium point for small $`\mathrm{\Sigma }`$ and $`Q`$,

$$\mathrm{\Sigma }\frac{5}{8}\left(4\frac{1}{\gamma _0}\right)Q,$$
(88)
so that for $`\gamma _0\frac{1}{4}`$ we have that $`\mathrm{\Sigma }`$ and $`Q`$ have the same sign (unlike the non-truncated theory case), and hence it is immediately clear that there will be a different qualitative dynamical behavior in the truncated theory.

In particular, close to the FLRW equilibrium point $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(1,0,0,0)`$ (which we will use to determine the initial conditions in our numerical analysis in the following section), from a calculation of the corresponding eigenvalues we have that:

$$\mathrm{\Omega }_{(m)}e^\tau ,(1\mathrm{\Omega }_{(r)})e^{2\tau },$$
(89)
$$\mathrm{\Sigma },Qe^{\alpha _+\tau }+\beta e^{\alpha _{}\tau },$$
(90)
from equations (56) - (59) in the non-truncated case, where $`\beta `$ is a constant, and

$$\alpha _\pm =\frac{1}{2\gamma _0}\left[\gamma _01\pm \sqrt{\left(\gamma _01\right)^24\gamma _0\left(1+8\gamma _0/5\right)}\right].$$
(91)
For physical values $`\gamma _01,`$ both $`\alpha _\pm `$ have negative real parts and as noted earlier, this FLRW equilibrium point is a saddle $`(\text{e. g., for}\gamma _0=1/2,\alpha _\pm =1/2\pm i\sqrt{67/20})`$.

In the truncated case we have that eqns. (89) and (90) are satisfied close to the FLRW equilibrium point, but now with

$$\alpha _\pm ^{\text{tr}}=\frac{1}{2\gamma _0}\left[5\gamma _01\pm \sqrt{\left(5\gamma _01\right)^2+4\gamma _0\left(12\gamma _0/51\right)}\right].$$
(92)
We first note that for $`\gamma _01/5,`$ at least one of the $`\alpha ^{\text{tr}}`$ has a positive real part which leads to a change of stability (indeed, for $`1/5\gamma _05/12,`$ this equilibrium point is a source !). Physically, this means that in the truncated case $`\mathrm{\Sigma }`$ and $`Q`$ in equation (90) have a growing mode, and hence their magnitudes increase leading to a breakdown in the physical model (and the time period for which the assumptions are valid). This breakdown is seen in the numerics (figures 7). As a comparison, for $`\gamma _0=1/2,\alpha _\pm \frac{3}{2}(1\pm 1.1);`$ the growing mode which evolves approximately as $`e^{3\tau },`$ leads to a rapid increase in the magnitudes of $`\mathrm{\Sigma }`$ and $`Q`$ and the models fail after a relatively short cosmological time (see figures 7).

## VIII Numerical Integration

From the dynamical analysis carried out in section V, there are no sources (at finite values) in the physical regime; that is, during the time period for which the various physical assumptions used here are valid there are no past attractors. It is reasonable to assume as conditions prevailing at the beginning of the regime we are interested in, that the universe is approximately isotropic and spatially homogeneous (e. g., almost FLRW with $`|Q|1`$ and $`\mathrm{\Sigma }1`$) and that the radiation component is dominant. These assumptions are consistent with current observations, and we are also assuming that some mechanism (such as, for example, inflation) at early times has driven the universe toward this configuration. Thus the model is close to the particular solution mentioned in (ii) with $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)`$ given by $`(1,0,0,0).`$ As noted above, this equilibrium point is a saddle, however it is a ’stronger’ attractor than other saddles in that it has more positive eigenvalues. Also, this equilibrium point lies in the invariant subspace $`\chi =0,`$ associated with the zero curvature Bianchi I sub-case. Since $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(1,0,0,0)`$ is an equilibrium point lying in an invariant set, if a model is driven toward this point in phase space it can stay an arbitrarily long time close to this point (i. e., the universe can spend an extended period close to this FLRW model). The universe will then eventually begin to evolve away from this configuration and from the invariant set as time evolves. Therefore, we shall assume initial conditions for our numerical integration based on the fact that the universe starts evolving from situations close to this equilibrium point.

### VIII.1 Initial conditions.

Since we must have (for physical reasons) $`Q<0`$ and $`Q\mathrm{\Sigma }<0,`$ initial conditions close to $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(1,0,0,0)`$ can be given by:

$`\mathrm{\Omega }_{(r)}(0)=1ϵ,\mathrm{\Omega }_{(m)}(0)=ϵ\delta >0,`$
$`\mathrm{\Sigma }(0)=\mathrm{\Sigma }_i>0,Q(0)=Q_i<0,`$ (93)
where $`ϵ,\delta ,\mathrm{\Sigma }_i,Q_i`$ are real small constants. The initial values of “total Omega”: $`\mathrm{\Omega }_{(tot)}=\mathrm{\Omega }_{(m)}+\mathrm{\Omega }_{(r)}`$ and $`\chi `$ given by eqn. (61) are:

$$\mathrm{\Omega }_{(tot)}(0)=1\delta ,\chi (0)=\delta (\mathrm{\Sigma }_i)^2,$$
(94)
so that the value of $`\delta `$ reflects the deviation of $`\mathrm{\Omega }_{(tot)}`$ from unity, while the deviation from the invariant set $`\chi =0`$ depends on $`\delta `$ and $`(\mathrm{\Sigma }_i)^2`$.

* The factor $`\gamma `$.

As discussed earlier, during the interactive range we are interested in, the various timescales must satisfy

$$t__H>\tau _{\text{rel}}>t_c>t_\gamma ,$$
(95)
where $`t__H=3/\mathrm{\Theta },\tau _{\text{rel}}`$ follows from eqn. (73) and $`t_c,t_\gamma `$ are given by eqns. (36) and (37). Considering that $`(\mathrm{\Theta }_i/3)^2(1/3)\kappa aT_i^4`$ and inserting the constants appearing in eqns. (36) and (37), we arrive at the following initial values:

$`\mathrm{log}_{10}(t__H)|_i8.8,\mathrm{log}_{10}(t_c)|_i7.6,`$
$`\mathrm{log}_{10}(t_\gamma )|_i3.4,`$ (96)
thus suggesting the range

$$0.06<\gamma _i<1$$
(97)
for the initial value $`\gamma (0)=\gamma _i`$<sup>2</sup><sup>2</sup>2We distinguish between the constant initial value $`\gamma _i`$ and the case in which $`\gamma =\gamma _0`$ for all the evolution period..

If assuming that $`\gamma =\gamma _0`$ for all times, then the near-Eckart regime can be associated with $`\gamma _0O(10^2),`$ while a transient regime follows by setting $`\gamma _0O(1)<1`$.
* The constants $`\mathrm{\Sigma }_i`$ and $`Q_i`$.

From evaluation of eqn. (LABEL:Sdef4) at $`\tau =0,`$ the constant $`Q_i^2`$ is proportional to the initial deviation of the photon entropy from its equilibrium value:

$$\frac{15\lambda _0}{8S_{(e)}}Q_i^2Q_i^21\frac{S(0)}{S_{(e)}(0)},$$
(98)
where we have used eqns. (34) and (35) so that $`\lambda _0/S_{(e)}O(1),`$ with $`\lambda _0`$ given by eqn. (LABEL:Sdef4). Both numbers $`\mathrm{\Sigma }_i,Q_i`$ are initial ratios of off-equilibrium and anisotropic variables ($`\mathrm{\Pi }_{ab},\sigma _{ab}`$) with respect to equilibrium and isotropic variables ($`p^{(r)},\mathrm{\Theta }/3`$). Since we are choosing initial conditions close to a saddle point associated with near FLRW conditions and we must assume near thermal equilibrium, then $`\mathrm{\Sigma }_i,|Q_i|`$ must necessarily be small numbers ($`1`$). A maximal bound on $`|Q_i|`$ and $`|\mathrm{\Sigma }_i|`$ can be fixed from CMB observations mes  making it reasonable to take

$$|Q_i|<0.010.1,|\mathrm{\Sigma }_i|<0.0010.01,$$
(99)
However, we will comment further ahead on the sensitivity of the functions to these initial values.
* The constants $`ϵ`$ and $`\delta `$.

The values for $`ϵ`$ and $`\delta `$ are restricted by the ratio of photons to WIMP’s. Considering the neutralino as the WIMP particle with $`m100`$ GeV, using the ideal gas law and eqn. (18) yields:

$$\frac{ϵ\delta }{1ϵ}=\frac{\mathrm{\Omega }_{(m)}(0)}{\mathrm{\Omega }_{(r)}(0)}=\frac{mc^2}{3k__BT_i}\nu _\text{w}0.013,$$
(100)
leading to the following constraint

$$\delta 1.013ϵ0.013,$$
(101)
that must be satisfied by all initial conditions compatible with $`T_i10^6`$ K and with the observational constraints $`\mathrm{\Omega }_\text{w}0.3\pm 0.1`$ and $`h0.7.`$ Further restrictions on $`ϵ`$ and $`\delta `$ follow by demanding that $`\mathrm{\Omega }_{(r)}`$ decreases and $`\mathrm{\Omega }_{(m)}`$ increases at the initial time $`\tau =0.`$ From the expressions for $`\mathrm{\Omega }_{(m)}^{}`$ and $`\mathrm{\Omega }_{(r)}^{}`$ in the differential equations, these conditions imply:

$$14\mathrm{\Sigma }_i^2<ϵ\delta <2\mathrm{\Sigma }_i|Q_i|4\mathrm{\Sigma }_i^2,$$

so that $`ϵ+\delta >0`$ must hold, leading to the following minimal values of $`ϵ`$ and $`\delta :`$

$$ϵ>0.0064,\delta >0.0065.$$
(102)
The condition that $`\mathrm{\Omega }_{(r)}`$ decreases at $`\tau =0,`$ together with $`Q_i<0`$ and $`\mathrm{\Sigma }_i>0,`$ imply:

$$0<\mathrm{\Sigma }_i<\frac{1}{4}\sqrt{Q_i^2+4(ϵ+\delta )}\frac{|Q_i|}{4}.$$
(103)

### VIII.2 A physically plausible evolution and the range of validity of the models.

It is important to specify a criterion in order to distinguish a physically plausible evolution for the models. We define such an evolution by the following conditions that must hold all along the range of validity discussed previously:

1. $`\mathrm{\Omega }_{(m)}`$ increases while $`\mathrm{\Omega }_{(r)}`$ decreases. The transition from a radiation- to a matter-dominated epoch occurs within the radiative era. However, the ratio $`\mathrm{\Omega }_{(m)}/\mathrm{\Omega }_{(r)}`$ must remain finite in all the validity range.
2. $`S`$ must be an increasing and convex function, tending asymptotically to the equilibrium photon entropy given by eqn. (35).
3. $`\dot{S}_{(e)}\mathrm{\Sigma }Q`$ must be very small. Idealy, we should have $`\dot{S}_{(e)}>0`$, though this condition might fail to hold as long as $`\dot{S}>0`$ holds (see jcl  for examples).
4. Initially, we have eqns. (95), (96), and (97), but then at later stages $`t_c`$ (the Compton scattering timescale) is no longer relevant, while $`t_\gamma `$ and $`\tau _{\text{rel}}`$ should overtake $`t__H,`$ so that the baryon–photon decoupling is defined as $`t_\gamma =t__H`$ and should occur at $`\tau =\tau __D`$ such that $`T(\tau __D)=T__D4\times 10^3`$ K.

We will not be concerned with the evolution of the models after the radiative era, since the assumptions regarding a hydrodynamic description of the radiative fluid break down. After the baryon–photon decoupling, an appropriate treatment of cosmic matter requires a different theoretical framework based on Kinetic Theory ellisvan .

Consider a transient regime ($`\gamma \gamma _0O(1)<1`$), together with “test” initial conditions given by eqns. (93) with $`\delta =0,`$ satisfying eqns. (99), (101), (102), and (103), hence $`ϵ=0.0128.`$ For the time being, $`\mathrm{\Sigma }_i`$ will be taken to be two orders of magnitude smaller than $`Q_i.`$ This yields the following initial conditions lying very near $`\chi (0)=0`$:

$`\mathrm{\Omega }_{(m)}(0)=0.0128,\mathrm{\Omega }_{(r)}(0)=0.9872,`$
$`\mathrm{\Sigma }(0)=0.001,Q(0)=0.1,`$
$`\text{so that:}\mathrm{\Omega }_{(tot)}(0)=1,\chi (0)=10^6.`$ (104)
In order to illustrate how the different variables should behave in a physical evolution taking place in the appropriate timescale, we integrate the system of governing equations (56), (58), (59), and (69) for (104) and $`\gamma =\gamma _0=0.7.`$ The results are displayed in figures 1 and 2. As shown by figures 1a, 1b, 1c and 1d, the functions $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},\mathrm{\Sigma },Q,S_{(e)}`$ and $`S`$ comply with conditions 1, 2 and 3 of the physical evolution mentioned above (we shall discuss condition 4 in section VIII-E). Figure 2a depicts joint logarithmic plots of the various timescales $`t_\gamma ,t_c,t__H,`$ and physical time $`t,`$ respectively given by eqns. (36), (37), (66) and (67), while the radiation temperature $`T`$ is displayed in figure 2b. By comparing figures 1 and 2, it is evident that the range of validity of the models is roughly $`0\tau <6,`$ corresponding to $`10^6`$ K $`>T>10^3`$ K, with $`t10^5`$ years (the physical time for the radiative era), while the transition from radiation to matter dominance taking place at about $`\tau =4`$ ($`T10^4`$ K). Numerical integration of the governing equations for initial conditions different from (104) might yield important qualitative changes in the state variables plotted in figures 1, like $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)}`$ or $`S,`$ but not of those plotted in figures 2, such as $`T`$ or the timescales (36), (37), (66), $`\tau _{\text{rel}}`$ or physical time $`t.`$ Therefore, figures 2 provide a general estimation of the range of validity of the models for a wide range of initial conditions. We will consider more general initial conditions in the following subsection.

### VIII.3 Sensitivity to deviations from $`\chi =0`$.

We will test the effect of initial deviations from the invariant set $`\chi =0`$ on state functions obtained by numerical integration of the governing equations. We consider initial conditions as in (93), keeping $`\mathrm{\Sigma }_i,`$ and $`Q_i`$ fixed, but now we take $`\delta 0`$:

$`\begin{array}{cc}\mathrm{\Omega }_{(m)}(0)=0.0128\delta ,& \mathrm{\Omega }_{(r)}(0)=0.9872,\\ \mathrm{\Sigma }(0)=0.001,& Q(0)=0.1,\\ \mathrm{\Omega }_{(tot)}(0)=1\delta ,& \chi (0)=\delta 10^6,\end{array}`$ (105)
Since $`\delta `$ can be given in terms of $`ϵ`$ by eqn. (101) and $`ϵ=0.0128`$ corresponds to $`\delta =0,`$ testing values of $`ϵ`$ near $`ϵ=0.0128`$ determines the initial deviation from $`\mathrm{\Omega }_{(tot)}=1`$ and $`\chi =0.`$ Notice that $`\delta `$ can be positive or negative, respectively, for $`ϵ>0.0128`$ or $`ϵ<0.0128,`$ so that $`\chi (0)=0`$ if $`\delta =\mathrm{\Sigma }_i^2=10^6,`$ while $`\chi (0)`$ is positive/negative if $`\delta >\mathrm{\Sigma }_i^2=10^6`$ or $`\delta <\mathrm{\Sigma }_i^2=10^6`$ (though, from eqn. (61), curvature has opposite sign to $`\chi `$). We integrate numerically the governing equations for initial conditions (105), $`\gamma =\gamma _0=0.7`$ (transient regime) and for various values of $`\delta .`$ The resulting forms of $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(tot)},`$ and $`S`$ are respectively plotted in figures 3a to 3d (for positive curvature $`\chi <0`$) and in figures 4a to 4d (for negative curvature $`\chi >0`$). These figures clearly show that a physical evolution is only possible for initial conditions that deviate very slightly from $`\chi (0)=0`$ (less than $`10^{5.5}`$), leading to orbits that remain very close to the invariant set $`\chi =0.`$ If we fix $`ϵ`$ and $`\delta `$ (for any combination of values compatible with (101), (102) and (103)) and vary $`\mathrm{\Sigma }_i,`$ so that deviation from $`\chi (0)=0`$ is governed by $`\mathrm{\Sigma }_i,`$ we obtain exactly the same behavior displayed by figures 3 and 4, leading to the same conclusion: a physical evolution is only possible for initial conditions for which $`|\chi (0)|10^{5.5}`$.

### VIII.4 The near-Eckart regime

We examine the near-Eckart regime by assuming $`\gamma =\gamma _0=0.0011,`$ together with initial conditions (104). As shown in figure 5a, the functions $`Q(\tau )`$ and $`S(\tau )`$ clearly have the form (76), indicating a quick relaxation in terms of an abrupt exponential decay (in about $`\tau 0.01`$). Figure 5b shows how the equilibrium entropy $`S_{(e)}`$ tends to a constant value in about $`\tau 0.02,`$ a longer time than the relaxation of $`S/S_{(e)},`$ thus indicating that the effective relaxation time is provided by $`S_{(e)},`$ and $`S_{(e)}^{}>0`$ holds for all the evolution, in agreement with the “local equilibrium hypothesis” (see section VI.1). Figure (5c) reveals how $`(5/8)Q/\mathrm{\Sigma }\gamma _0=0.001,`$ so that Eckart’s transport equation ($`P+2\eta \sigma =0`$ with $`\eta `$ given by eqn. (29)) is approximately valid once the quick relaxation is over. A good approximation to the relaxation time in a near-Eckart regime is given by the assumption (73) with $`\gamma =(5/8)\gamma _0(Q/\mathrm{\Sigma })`$ and $`\gamma _0=0.001.`$ Figure 5d illustrates that this is a correct assumption, since the obtained $`\tau _{\text{rel}}`$ overtakes $`t__H`$ in the very short period that coincides with the duration of the relaxation process ($`\tau 0.02`$).

The plots of the functions $`\mathrm{\Omega }_{(m)}(\tau )`$ and $`\mathrm{\Omega }_{(r)}(\tau )`$ are identical to those that would have resulted in the transient regime had we chosen initial conditions (104) and a much larger value of $`\gamma _0.`$ The function $`\mathrm{\Sigma }(\tau )`$ is affected, becoming almost constant for very small $`\gamma _0.`$ This is reasonable, since eqn. (56) does not contain $`Q,`$ while eqns. (58) and (59) do contain this function, but $`|Q(\tau )|`$ and $`Q^{}`$ decay very fast becoming almost zero for most of the time range and so the differential equations for the functions $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ (but not $`\mathrm{\Sigma }`$) are practically unaffected. The effect of varying $`\gamma _0`$ (as long as we have small values, $`<0.01`$) is simply to make the decay of $`Q`$ and $`S_{(e)}^{}`$ slightly more or less abrupt (depending on whether $`\gamma _0`$ is smaller or larger than $`0.001`$) and has no noticeable effect on numerical curves of other functions.

### VIII.5 Testing the relaxation times numerically.

The assumption $`\gamma =\gamma _0<1`$ leads to a reasonable relaxation time only in the earlier stages of the evolution, while a choice $`\gamma =\gamma _0>1`$ might work for later conditions (near matter-radiation decoupling) but not for earlier times. This can be appreciated in figure 2b, since the relaxation times obtained for different values of $`\gamma =\gamma _0`$ would be curves parallel to $`t__H.`$ As we mentioned before, $`\gamma _0`$ controls the rate of decay of $`Q`$ and $`S,`$ but it is still interesting to check if other functions are sensitive to changes in the numerical constant value of $`\gamma _0.`$ Assuming initial conditions (104) and integrating the dynamical system for various values of $`\gamma _0>0.1`$ (thus excluding the non-transient zone) shows that the other functions (such as $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)}`$ or $`\mathrm{\Sigma }`$) are essentially identical to those shown in figures 1 to 4, being thus unaffected by the numerical value of $`\gamma _0`$.

Since $`\gamma `$ constant is not very realistic, we test now the expression for $`\tau _{\text{rel}}`$ that follows from equations (82) and (83) in section VI.2. Ignoring the curvature term, we have that $`\mu _0/\gamma `$ equals

$$\frac{8}{5}+\mu _0(1+\mu _0)+\frac{\mu _0(12\mu _0)\mathrm{\Omega }_{(m)}}{2}\mu _0\left[1+(\lambda 3)\mu _0\right]\mathrm{\Sigma }^2,$$
(106)
It is evident that choosing a sufficiently small $`\mu _0`$ yields $`\gamma _i5\mu _0/81,`$ an initial value characteristic of a near-Eckart regime. As expected, numerical tests with $`\mu _01`$ lead to practically identical curves as those corresponding to the near-Eckart regime (figures 5a, 5b and 5c). However, for $`\mu _0>1/2,`$ the factor multiplying $`\mathrm{\Omega }_{(m)}`$ in the denominator of (106) can become negative for some $`\tau ,`$ thus opening the possibility that the denominator might become small and so $`\gamma `$ might increase to $`O(1)`$ as $`\tau `$ reaches later times related to the baryon–photon decoupling era. In order to explore this possibility, we integrated the dynamical system for initial conditions (104) and under the assumption of $`\gamma `$ given by eqn. (106), with various values of $`\mu _0>1/2.`$ As shown by figure 6a, the choice $`\mu _0>15`$ leads to $`\gamma `$ diverging as (approximately) $`\tau 7,`$ this leads to $`\tau _{\text{rel}}`$ overtaking $`t__H`$ (figure 6b) around $`\tau 6.8.`$ It would have been nicer to have $`\tau _{\text{rel}}`$ overtaking $`t__H`$ at an earlier time (say $`\tau 55.5`$) as required by condition 4 for a physically plausible evolution, but we feel that the form of $`\tau _{\text{rel}}`$ associated with (106) is a reasonable approximation to a physical relaxation time that acts as an “effective” relaxation time for the radiative era. Finally, another consequence of dealing with a more reasonable form for $`\tau _{\text{rel}}`$ is the fact that $`\dot{S}_{(e)}\mathrm{\Sigma }Q`$ being positive and very small (condition 3 of a physically plausible evolution) is better satisfied than in the case of constant $`\gamma _0`$ (compare figures 1b and 6c).

### VIII.6 The truncated equation.

Considering initial conditions (104) and $`\gamma =\gamma _0=0.7,`$ the integration of the governing equations (56), (58), (59), and (86) yield the curves for $`Q,\mathrm{\Sigma },\mathrm{\Omega }_{(m)},`$ and $`\mathrm{\Omega }_{(r)}`$ depicted in figures 7a, 7b and 7c. These figures reveal that a physical evolution fails to occur, as the growing modes predicted by the qualitative analysis clearly emerge, making all these functions undertake an unphysical growth. We tested other values of $`\gamma _0`$ and $`\gamma `$ given by eqn. (106), as well as other initial conditions and obtained very similar curves to those of figures 7a-c, all failing to comply with the criteria for a physically plausible evolution.

### VIII.7 A baryonic scenario without CDM

If we had considered only the radiative fluid (baryons, photons and electrons) as the matter source of the models, then the bulk of the rest mass energy density would have been due to the baryons, so that $`n^{(m)}`$ would have been identified with $`n_\text{b}`$ instead of $`n_\text{w}`$ and $`\mathrm{\Omega }_{(m)}`$ would correspond to the baryonic $`\mathrm{\Omega }_\text{b}.`$ It is interesting to examine numerically the consequences of this “baryonic scenario”. Considering the baryon–photon number ratio in eqn. (15), the baryonic scenario implies replacing eqns. (100) and (101) by

$$\frac{ϵ\delta }{1ϵ}=\frac{m_\text{b}c^2}{3k__BT_i}\nu _\text{b}0.002,\delta =1.002ϵ0.002,$$
(107)
while initial conditions are then given by eqns. (104) and (105), but $`\delta =0`$ now corresponds to $`ϵ=0.00199`$ instead of $`ϵ=0.0128.`$ Intuitively, we do not expect a mayor qualitative change in the resulting graphs, though it is reasonable to expect that $`\mathrm{\Omega }_{(m)}`$ will be smaller and $`\mathrm{\Omega }_{(r)}`$ larger, since baryons have less rest mass density (by one or two orders of magnitude) than WIMP’s and so it should take longer for baryons to dominate over radiation. The numerical curves that result are as expected intuitively, with $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(tot)},`$ and $`S/S_{(e)}`$ having very similar forms as the curves of figures 3 and 4, with $`\mathrm{\Omega }_{(r)}`$ decreasing slightly slower than in the case with WIMP’s. Since the obtained curves for $`\mathrm{\Omega }_{(m)}`$ in the baryonic scenario are so close to those obtained in figures 3 and 4 for the case with WIMP’s, these curves yield baryon abundances that are clearly incompatible with the bounds placed by cosmic nucleosynthesis ($`\mathrm{\Omega }_\text{b}10^2`$).

## IX Discussion and Conclusion.

We have studied a class of dissipative Kantowski-Sachs models describing the cosmological evolution during the radiative era characterized by radiative processes involving baryons, electrons and photons, considered as a single dissipative “radiative” fluid. We also assumed the presence of CDM, in the form of a non-relativistic gas of WIMP’s (lightest neutralinos). Although this gas does not interact with the radiative fluid, it provides the bulk of the rest mass energy density and thus it strongly influences the dynamics of the models and the resulting values of cosmological timescales, such as $`t__H.`$ On the other hand, the radiative fluid provides the bulk of thermal and dissipative effects, related to the rate of change and relaxation of the radiation entropy to its equilibrium value.

After defining new normalized variables $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},`$ $`\mathrm{\Sigma },Q,`$ a set of evolution equations has been derived based upon appropriate thermodynamical laws and equations of state. The qualitative and numerical study of these evolution equations has clarified various aspects of the dynamical behavior of the models, their physical viability, as well as a peculiar sensitivity to certain initial conditions related to deviation from the invariant set $`\chi =0.`$ We discuss below the main features emerging from previous sections.

The definition of the phase space variables $`\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},`$ $`\mathrm{\Sigma },Q`$ leads in a natural way to expressing the relaxation time $`\tau _{\text{rel}}`$ as proportional to $`t__H=3/\mathrm{\Theta }`$ (see equations (68) and (73)). The understanding of the relaxation process can be accomplished by studying the effect of different choices of the proportionality factor, $`\gamma (\tau )>0,`$ on the exponential decay of the dissipative stress (related to $`Q`$) and of the photon entropy $`S`$ to its equilibrium value $`S_{(e)}.`$ We have identified a “near-Eckart” regime if this decay is abrupt ($`\gamma 1,`$ and becoming instantaneous in the limit $`\gamma 0,`$ so that $`\tau _{\text{rel}}0`$), while a “transient regime” ($`\gamma O(1)<1`$) can be associated with a slower decay. Both the near-Eckart and the transient regimes are compatible with a physically plausible evolution. The difference between the two regimes is the timescale of their relaxation process: for the transient regime this timescale can be comparable with the duration of the radiative era, for the near-Eckart regime it is much shorter (about eight orders of magnitude in physical time). This is well illustrated by the differences in evolution timescales between figures 5 and figures 1 to 4. Comparing figures 2b and 5d, it is evident that the relaxation timescale of the near-Eckart regime ($`10^8`$ sec) is much shorter than that of the Compton scattering ($`10^{11}`$ sec). Therefore the near-Eckart regime yields a relaxation that is too swift and so it is inadequate to examine the two main radiative processes of the radiative era: the Compton scattering and (more so) the Thomson scattering. It is important to mention this fact in view of recent claims geroch  that a transient theory of irreversible thermodynamics is not really necessary (see pavher  for a comprehensive discussion). It is evident that the cosmological study of radiative processes in pre-decoupling times needs to be accomplished with a transient regime.

As revealed by figures 1, 2, 3 and 4, a physically plausible evolution is possible for all the duration of the radiative era for initial conditions given by an initial state very close to $`\chi (0)=0,`$ hence lying very close to the invariant set of zero curvature $`\chi =0.`$ All models complying with a physical evolution begin their evolution near the equilibrium point $`(\mathrm{\Omega }_{(r)},\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(1,0,0,0),`$ a saddle with positive eigenvalues (i. e., stable) associated with the FLRW sub-case of the models, and the proceeding evolution remains very close to the invariant set $`\chi =0.`$ This is an extremely interesting feature of these models, as it relates a geometric property of KS solutions with the constraints imposed by the physics and by observational evidence, since recent data from CMB observations indicates $`\mathrm{\Omega }_{(tot)}1.`$ By looking at the curves of $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(tot)}`$ with $`\delta 10^5`$ in figures 4a and 4c (negative curvature), it is evident that for $`\tau >6`$ these curves decrease from their values $`\mathrm{\Omega }_{(m)}1`$ and $`\mathrm{\Omega }_{(tot)}1`$ around $`\tau =6.`$ Had we plotted these curves for larger values of $`\tau ,`$, extending to the present era ($`\tau 15`$), we would have obtained $`\mathrm{\Omega }_{(m)}\mathrm{\Omega }_{(tot)}0.3,`$ in agreement with the currently accepted value of $`\mathrm{\Omega }_{(m)},`$ but not of $`\mathrm{\Omega }_{(tot)}.`$ Of course, the estimated contribution to $`\mathrm{\Omega }_{(tot)},`$ today, for non-relativistic matter (CDM plus baryons) is only $`0.3,`$ with the remaining two thirds of the critical density probably related to a $`\mathrm{\Lambda }`$-type “dark energy” interaction whose precise nature and properties are still uncertain. However, this discrepancy with regards to $`\mathrm{\Omega }_{(tot)}`$ today is not surprising since we did not consider any $`\mathrm{\Lambda }`$-type interaction, and so it does not affect our results. Also, the models we are considering are only valid for a specific range of cosmological times: $`10^3z10^6,`$ in which this “dark energy” would likely not have been dominant. Still, the close link between a physically plausible evolution and $`\mathrm{\Omega }_{(tot)}`$ near unity is remarkable.

It is interesting to compare our results to those reported previously weber  dealing with the perfect fluid sub-case of the KS models examined in this paper (though, these models did not consider a CDM component). As reported in weber , there are numerical solutions in which both matter and radiation normalized densities, $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},`$ decay to zero as the models re-collapse and approach a crunch singularity. By looking at the forms of $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ in figures 4, it is evident that such evolution is similar to that depicted by curves associated with initial conditions $`\delta 10^3.`$ However, the evolution that results from these initial conditions fails to comply with our physical criteria, since the entropy $`S`$ is no longer a convex function for all of the time range (see figure 4d) and starts decreasing at too early times. In the perfect fluid case, these examples satisfy an appropriate equation of state and all of the energy conditions and also the photon entropy is simply $`S_{(e)}`$ and is constant for all times, hence there is no physical reason to discard these curves (other than remarking that such behavior of $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)}`$ is not observed in the real universe). However, for the dissipative source under examination here, $`S`$ and $`S_{(e)}`$ are not constant and the conditions for a physical evolution are more stringent, hence we can apply clear physical criteria to discard these perfect fluid “strange” cases.

For initial conditions near the equilibrium point $`(\mathrm{\Omega }_{(r)},`$ $`\mathrm{\Omega }_{(m)},\mathrm{\Sigma },Q)=(1,0,0,0)`$ (as in equation (104)), the numerical curves of the “equilibrium variables” $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ are not affected by the choice of $`\gamma =\gamma _0,`$ the constant proportionality factor between $`\tau _{\text{rel}}`$ and $`t__H.`$ This is a consequence of the fact that for these initial conditions the values of $`\mathrm{\Sigma },Q`$ and $`SS_{(e)}`$ remain small during all of the evolution. Therefore, the evolution equations for $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ (equations (56) and (58)) are practically unaffected by the presence of $`\mathrm{\Sigma }`$ and $`Q,`$ and so are insensitive to the rate of transiency given by $`\gamma `$ constant. Since $`\mathrm{\Sigma }`$ and $`Q`$ govern the deviation from the FLRW equilibrium point, this decoupling of $`\mathrm{\Omega }_{(m)}`$ and $`\mathrm{\Omega }_{(r)}`$ from $`\mathrm{\Sigma }`$ and $`Q`$ is then a consequence of the evolution of the system always remaining close to thermal equilibrium. We tested this behavior for a particular form of $`\gamma `$ (eqn. (106), section VI.2): the curves for $`\mathrm{\Omega }_{(m)},\mathrm{\Omega }_{(r)},\mathrm{\Sigma },`$ and $`Q`$ are practically identical with those that follow from choosing $`\gamma =0.7`$ in figures 1 to 4. However, as shown in figures 6a and 6b, for large values of $`\mu _0`$ defined by (83), the obtained relaxation time $`\tau _{\text{rel}}`$ behaves similarly to what one would expect of a relaxation parameter for the radiative era. Although it has become common practice to simply equate $`\tau _{\text{rel}}`$ with a microscopic interaction time, like $`t_\gamma `$ or $`t_c,`$ the relaxation time is not a microscopic but a mesoscopic or even macroscopic quantity (though it must be qualitatively analogous to interaction times pavher ). Since it can be extremely cumbersome to evaluate $`\tau _{\text{rel}},`$ it is useful to have a concrete example where this relaxation parameter can be adequately approximated by the same dynamical equations associated with the models.

Finally, by means of qualitative arguments supported by the numerical analysis, we have shown in section VII (figures 7), that the truncated equation does not comply with a physically plausible evolution. This is an important result, since we have found a concrete example in which a truncated transport equation leads to unphysical evolution of dissipative fluxes. Although this conclusions strictly applies to the KS models under consideration, we must point out that one should be very cautious when applying these equations to other models and other equations of state.

A possible extension of this work would be to consider instead of CDM other forms of dark matter, such as “warm” dark matter (WDM) or axions. Another possibility is to include, together with dark matter, a scalar field associated with “dark energy”. Another route to generalize the present work is use a class of metrics associated with a geometry that is less restrictive than KS, for example the non-static spherical symmetry (perhaps under the assumption of self–similarity). We regard the present analysis of the Kantowski-Sachs models as a first step toward an understanding of the dynamics of cosmic matter in more general and physically motivated inhomogeneous models.

## Acknowledgements

This work has been partially supported by the Universidad Nacional Autónoma de México (UNAM) under grant DGAPA-IN-122498.