# Backreaction effects of dissipation in neutrino decoupling

## I Introduction

Non-equilibrium processes in the early universe are typically associated with dynamical transitions or particle decouplings. In the case of neutrino decoupling, the standard approach is to treat the process as adiabatic (see e.g. ). The small non-equilibrium effects are thus usually neglected, which provides a reasonable approximation. However, given the increasing accuracy of cosmological observations and theoretical modeling, it is worthwhile revisiting the standard equilibrium models of processes such as neutrino decoupling, in order to see whether non-equilibrium corrections can lead to observable consequences. Recently, non-equilibrium corrections in neutrino decoupling have been calculated in a number of papers, using complicated kinetic theory and numerical computations (see  for a short review). The corrections are very small, as expected. For example, in  it was found that non-equilibrium effects lead to a small change in the decoupling temperature for neutrinos. Spectral distortions have also been analyzed , showing the remarkable fact that they amount to as much as 1% or more for the higher-energy side of the spectrum. Although these corrections in the spectrum, energy density and temperature of the neutrino component have hardly any effect on primordial helium synthesis, yielding a change in the mass fraction of $`10^4`$, they can lead to other effects that may be observable. Thus it is shown that the non-equilibrium increase in neutrino temperature, which leads to an extra injection of energy into the photon spectrum, leads to a shift of equilibrium epoch between matter and radiation which, in turn, modifies the angular spectrum of fluctuations of the cosmic microwave background radiation .

Despite the accuracy of these models in obtaining corrections to the decoupling temperature and distribution function due to non-equilibrium effects, they still make use of the standard Friedman equations for a perfect (i.e non-dissipative) fluid. This leads to the physically inconsistent situation in which, say, the energy density and expansion evolve in time like a radiative fluid in equilibrium. One expects that small distortions in the particle equilibrium distribution function should be reflected in the macroscopic (i.e fluid) description, as given by the stress-energy tensor, by adding a bulk viscous pressure to the equilibrium one. Here we consider an alternative thermo-hydrodynamic model of dissipative effects in neutrino decoupling, simple enough to produce analytic solutions for the backreaction effects on the universal scale factor, and estimates for the entropy production due to dissipation. As explained above these effects are not the focus of recent papers, which use sophisticated kinetic theory models focusing on the neutrino temperature. Our simplified approach cannot compete with these models for accuracy and completeness, but it has the advantage of simplicity, allowing for a qualitative understanding of effects not previously investigated in detail. A similar approach has previously been developed in  to the reheating era that follows inflation.

The thermo-hydrodynamic model is based on an approximation to kinetic theory which respects relativistic causality. This approximation is the Grad moment method, leading to the causal thermodynamics of Israel and Stewart  in the hydrodynamic regime (see also  for an alternative but equivalent approach). This causal theory is a generalization of the more commonly used relativistic Navier-Stokes-Fourier theory. The latter, due to Eckart , may be derived via the Chapman-Enskog approximation in kinetic theory. The resulting theory is quasi-stationary and noncausal, and suffers from the pathologies of infinite wavefront speeds and instability of all equilibrium states . The main new ingredient in the causal transport equations is a transient term which contains the relaxation time. Our simple model is based on a one-component fluid. In , relaxation time processes are incorporated in a two-fluid model. In this setting, electrons and positrons on the one side and neutrinos and antineutrinos on the other side, are found to be in two different equilibrium states with slightly different temperatures. The system evolves towards a state of thermal equilibrium in a characteristic relaxation time.

Dissipative effects in the decoupling of a given species of particles arise from the growing mean free path of the decoupling particles in their weakening interaction with the cosmic fluid. Eventually the mean collision time exceeds the gravitational expansion time, and decoupling is complete. A hydrodynamic model may be used to cover the early stages of the decoupling process, but it will eventually break down when the mean collision time becomes large enough .

In the conditions prevailing at the time of neutrino decoupling, it is reasonable to neglect sub-horizon metric fluctuations and treat the spacetime as a Friedmann model. (The incorporation of perturbations in our model would use the covariant formalism for dissipative fluids developed in .) The dynamical effects of spatial curvature and any surviving vacuum energy will be negligible, so that we can reasonably assume a spatially flat geometry. Furthermore, we assume that the average 4-velocities of the neutrinos (regarded as massless) and of the photon-electron-positron gas are the same. With all these assumptions, only scalar dissipation is possible. Dissipation during neutrino decoupling arises because the falling temperature lowers the interaction rate with leptons as the lepton mass can no longer be ignored relative to the thermal energy. Thus dissipation is directly reflected in a deviation of the equation of state from the thermalized radiation form $`p=\frac{1}{3}\rho `$. Within a hydrodynamic one-fluid model, such dissipation is described via bulk viscosity, which vanishes in the $`p=\frac{1}{3}\rho `$ limit, but is nonzero otherwise. We will use the full (i.e. non-truncated) version of the causal transport equation for bulk stress.

## II Causal transport equation for bulk stress

The particle number 4-current and the energy-momentum tensor are

$`N^a=nu^a,T^{ab}=\rho u^au^b+(p+\mathrm{\Pi })h^{ab},`$

where $`\rho `$ is the energy density, $`p`$ is the equilibrium (hydrostatic) pressure, $`n`$ is the particle number density, $`\mathrm{\Pi }`$ is the bulk viscous pressure, and $`h^{ab}=g^{ab}+u^au^b`$ is the projector into the comoving instantaneous rest space. Particle and energy-momentum conservation

$`_aN^a=0,_bT^{ab}=0,`$

lead to the equations

$`\dot{n}+3Hn=0,`$ (1)   
$`\dot{\rho }+3H(\rho +p+\mathrm{\Pi })=0,`$ (2)   
where $`H`$ is the Hubble expansion rate. The specific entropy $`s`$ and the temperature $`T`$ are related via the Gibbs equation

$$nTds=d\rho \frac{\rho +p}{n}dn.$$
(3)   
Then it follows that

$$nT\dot{s}=3H\mathrm{\Pi },$$
(4)   
where $`\mathrm{\Pi }`$ is always non-positive. The Grad moment approximation in kinetic theory (or phenomenological arguments) leads to the full causal transport equation  for $`\mathrm{\Pi }`$:

$$\tau \dot{\mathrm{\Pi }}+\mathrm{\Pi }=3\zeta H\frac{1}{2}\tau \mathrm{\Pi }\left[3H+\frac{\dot{\tau }}{\tau }\frac{\dot{\zeta }}{\zeta }\frac{\dot{T}}{T}\right],$$
(5)   
where $`\tau `$ is the relaxation time scale, which allows for causal propagation of viscous signals, and $`\zeta 0`$ is the bulk viscous coefficient as given below. Quasi-stationary, noncausal theories have $`\tau =0`$, which reduces the evolution equation (5) to an algebraic equation $`\mathrm{\Pi }=3\zeta H`$. This leads to instantaneous propagation of viscous signals. Note also that the causal relaxational effects lead to a small increase in the sound speed over its adiabatic value :

$$c_\mathrm{s}^2c_\mathrm{s}^2+c_\mathrm{b}^2\text{ where }c_\mathrm{b}^2=\frac{\zeta }{(\rho +p)\tau }.$$
(6)   
This result, which is not well known, is derived in the appendix.

The approximation used in deriving the transport equation (also in the quasi-stationary case) requires that $`|\mathrm{\Pi }|\rho `$, which is reasonable for most dissipative processes (see  for a nonlinear generalization of the causal transport equation.)

Equation (5) as it stands is known as the full or non-truncated transport equation for bulk viscous pressure . When the term containing the square bracket on the right is neglected, we get the truncated equation which is usually used. Under many conditions, truncation leads to a reasonable approximation. We will use the full equation.

Taking $`n`$ and $`\rho `$ as independent variables, the Gibbs equation (3) leads to the integrability condition

$$n\left(\frac{T}{n}\right)_\rho +(\rho +p)\left(\frac{T}{\rho }\right)_n=T\left(\frac{p}{\rho }\right)_n,$$
(7)   
and together with the energy conservation equation (2) this gives the temperature evolution equation

$$\frac{\dot{T}}{T}=3H\left(\frac{p}{\rho }\right)_n\frac{1}{T}\left(\frac{T}{\rho }\right)_n3H\mathrm{\Pi }.$$
(8)   
The first term on the right accounts for adiabatic cooling due to expansion, whereas in the second term, viscosity contributes to heating of the fluid (note that $`\mathrm{\Pi }`$ is always non-positive).

Using equations (1) and (2), the Gibbs equation takes the form

$$n^2Tds=\left[\frac{n3H\mathrm{\Pi }}{3H(\rho +p)+3H\mathrm{\Pi }}\right]d\rho +(\rho +p)\left(\frac{n}{p}\right)_\rho \left[\frac{\dot{p}}{\dot{\rho }}d\rho dp\right].$$
(9)   
As expected we learn from the last equation that when the fluid is perfect ($`\mathrm{\Pi }=0`$), the specific entropy is conserved along the flow lines ($`\dot{s}=0`$). Furthermore, if a barotropic equation of state for $`n`$ holds, i.e. $`n=n(\rho )`$, then $`ds=0`$ so that $`s`$ is a universal constant, the same on all flow-lines, and the fluid is called isentropic.<sup>*</sup><sup>*</sup>* The same reasoning applies when the temperature is barotropic.  Yet, as Eq. (9) shows, this is no longer true in the presence of dissipation, i.e. a barotropic particle number density no longer forces $`ds`$ to vanish.

For simplicity, we assume the linear barotropic equation of state

$$p=(\gamma 1)\rho ,$$
(10)   
where $`\gamma `$ is constant and we are interested in the case $`\gamma \frac{4}{3}`$. The adiabatic speed of sound $`c_\mathrm{s}`$ is given by

$`c_\mathrm{s}^2=\left({\displaystyle \frac{p}{\rho }}\right)_s,`$

which for a perfect fluid (either barotropic or not) becomes

$`c_\mathrm{s}^2={\displaystyle \frac{\dot{p}}{\dot{\rho }}}.`$

When Eq. (10) holds then $`c_\mathrm{s}=\sqrt{\gamma 1}`$. Using Eq. (10) and the integrability condition (7), we find

$$T=\rho ^{(\gamma 1)/\gamma }F\left(\frac{\rho ^{1/\gamma }}{n}\right),$$
(11)   
where $`F`$ is an arbitrary function which satisfies $`\dot{F}=0`$. If $`T`$ is barotropic, then $`F`$ is constant and we have a power-law form with fixed exponent for the temperature 

$$T\rho ^{(\gamma 1)/\gamma }.$$
(12)   
In the non-dissipative case, these barotropic equations for $`p`$ and $`T`$ are compatible with the ideal gas law

$$p=nT,$$
(13)   
but in the presence of dissipation this is no longer true. In effect, equations (10), (12) and (13) imply $`n\rho ^{1/\gamma }`$, i.e.

$`{\displaystyle \frac{\dot{n}}{n}}={\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\dot{\rho }}{\rho }},`$

which implies, by using Eq. (2), that $`\mathrm{\Pi }=0`$. We shall drop in the sequel a barotropic equation of state for the temperature in favour of the more physically appealing equation of state (13) together the $`\gamma `$-law in (10).

## III Dissipation in neutrino decoupling

A hydrodynamic approach in the expanding universe requires a particle collision time $`t_\mathrm{c}`$ short enough to adjust to the falling temperature. As the natural time-scale for the expanding universe is $`H^1`$, we have

$`t_\mathrm{c}<H^1.`$

If $`t_\mathrm{c}H^1`$, then an equilibrium state can in principle be attained. Dissipative phenomena could play a prominent role for $`t_\mathrm{c}H^1`$.

We learn from kinetic theory that $`t_\mathrm{c}`$ is determined by

$$t_\mathrm{c}=\frac{1}{n\sigma v},$$
(14)   
where $`n`$ is the number density of the target particles with which the given species is interacting, $`\sigma `$ the cross-section and $`v`$ the mean relative speed of interacting particles. For the decoupling of massless neutrinos in the early universe, $`v=1`$, the target number density is that of electrons, and 

$`\sigma G__FT^2,`$

where $`G__F`$ is the Fermi coupling constant. At the neutrino decoupling temperature $`T_\mathrm{d}`$, we have $`m_\mathrm{e}/T_\mathrm{d}\frac{1}{2}`$, so that the rest mass energy $`m_\mathrm{e}`$ of electrons starts to become important. Since the electron number density in the radiation dominated era evolves as $`n_\mathrm{e}a^3`$, where $`a`$ is the scale factor, we have from Eq. (14) that

$$t_\mathrm{c}\frac{a^3}{T^2}.$$
(15)   
Dissipation due to massless particles with long mean free path in a hydrodynamic fluid is described by the radiative transfer model. The bulk viscous coefficient takes the form 

$$\zeta =4rT^4\mathrm{\Gamma }^2t_\mathrm{c},$$
(16)   
where $`r`$ is $`\frac{7}{8}`$ times the radiation constant and $`\mathrm{\Gamma }`$ measures the deviation of $`p/\rho `$ from its pure-radiation value:

$$\mathrm{\Gamma }=\frac{1}{3}\left(\frac{p}{\rho }\right)_n,$$
(17)   
where $`p`$ and $`\rho `$ refer to the pressure and energy density of the radiation/matter mixture as a whole. Since we assume the linear equation of state (10), it follows that $`\mathrm{\Gamma }`$ is a perturbative constant parameter in our simple model:

$`\mathrm{\Gamma }=\frac{4}{3}\gamma 1.`$

The assumption that $`\mathrm{\Gamma }`$ is constant relies on the assumption that decoupling takes place rapidly. Since standard adiabatic treatments of decoupling  assume instantaneous decoupling, this assumption should be a reasonable first approximation.

We may neglect the $`3\zeta H`$ term on the right of the transport equation (5), since it is $`O(\mathrm{\Gamma }^2)`$. Note that our simple model would thus break down in the quasi-stationary Eckart theory, since it would immediately lead to $`\mathrm{\Pi }=O(\mathrm{\Gamma }^2)`$. The relaxation timescale $`\tau `$ in causal radiative transfer  is given by $`\tau =t_\mathrm{c}`$. The term $`\dot{\zeta }/\zeta `$ on the right of Eq. (5) becomes

$`{\displaystyle \frac{\dot{\zeta }}{\zeta }}=H+O(\mathrm{\Gamma }),`$

on using equations (8) and (15). The full transport equation (5) becomes, to lowest order

$$\tau \dot{\mathrm{\Pi }}+\mathrm{\Pi }=4\tau H\mathrm{\Pi }.$$
(18)   
(We can think of the right hand side as an effective source term relative to the truncated transport equation.) We can rewrite this in the standard truncated form as

$$\tau _{}\dot{\mathrm{\Pi }}+\mathrm{\Pi }=0,$$
(19)   
where the effective relaxation time acquires an expansion correction:

$$\tau _{}=\frac{\tau }{1+4\tau H}.$$
(20)   
The amount of reduction depends on the size of $`\tau =t_\mathrm{c}`$ relative to $`H`$. The hydrodynamical description requires $`\tau H<1`$. If $`\tau H1`$, then $`\tau _{}\tau `$. But if $`\tau H`$ is close to 1, the reduction could be significant.

The Friedmann equation

$$\rho =3H^2,$$
(21)   
together with Eq. (2) leads to

$$\mathrm{\Pi }=2\dot{H}(43\mathrm{\Gamma })H^2.$$
(22)   
On using equation (22) we get from (18) the evolution equation for $`H`$

$$\ddot{H}+H\dot{H}(83\mathrm{\Gamma }+N)+H^3(2\frac{3}{2}\mathrm{\Gamma })(N+4)=0,$$
(23)   
where

$$N=(\tau H)^1,$$
(24)   
which is of the order of the number of interactions in an expansion time. Now, from equations (10), (13), (15) and (24) we have

$$N=\left(\frac{Ha}{H_\mathrm{d}a_\mathrm{d}}\right)^3,$$
(25)   
where the expression $`na^3`$ has been used and $`a_\mathrm{d}`$ and $`H_\mathrm{d}=H(a_\mathrm{d})`$ are the values at which $`N=1`$, so that $`a_\mathrm{d}`$ is determined by the equation

$$t_\mathrm{c}(a_\mathrm{d})H(a_\mathrm{d})=1.$$
(26)   
Changing the independent variable to the scale factor $`a`$, developing equation (23) and collecting the previous results, yields

$`a^2HH^{\prime \prime }+a^2H^2+aHH^{}\left[93\mathrm{\Gamma }+\left({\displaystyle \frac{Ha}{H_\mathrm{d}a_\mathrm{d}}}\right)^3\right]`$ (27)   
$`+\left(2\frac{3}{2}\mathrm{\Gamma }\right)H^2\left[4+\left({\displaystyle \frac{Ha}{H_\mathrm{d}a_\mathrm{d}}}\right)^3\right]=0,`$ (28)   
where a prime denotes $`d/da`$. We expand $`H`$ as

$$H=\overline{H}+\delta H\text{ where }\delta H=\mathrm{\Gamma }h+O(\mathrm{\Gamma }^2).$$
(29)   
The equilibrium Hubble rate $`\overline{H}`$ corresponds to the thermalized radiation state $`p=\frac{1}{3}\rho `$, so that $`\mathrm{\Gamma }=0`$, and Eq. (29) becomes

$`a^2\overline{H}\overline{H}^{\prime \prime }+a^2\overline{H}^2+9a\overline{H}\overline{H}^{}+8\overline{H}^2+\left(a\overline{H}\overline{H}^{}+2\overline{H}^2\right)\left({\displaystyle \frac{\overline{H}a}{\overline{H}_\mathrm{d}a_\mathrm{d}}}\right)=0.`$

The unique power-law solution is the well-known perfect radiative solution

$$\overline{H}=H_0\left(\frac{a_0}{a}\right)^2=\frac{1}{2t},$$
(30)   
where $`a_0`$ marks the start of the dissipative decoupling process, so that $`H=\overline{H}`$ for $`a<a_0`$.

Substituting Eq. (29) into (28) and using the fact that

$`{\displaystyle \frac{H_0a_0}{H_\mathrm{d}a_\mathrm{d}}}={\displaystyle \frac{a_\mathrm{d}}{a_0}}+O(\mathrm{\Gamma }),`$

we find that to $`O(\mathrm{\Gamma })`$:

$$a^2h^{\prime \prime }+a\left[5+\left(\frac{a_\mathrm{d}}{a}\right)^3\right]h^{}+\left[4+2\left(\frac{a_\mathrm{d}}{a}\right)^3\right]h=\frac{3}{2}H_0\left(\frac{a_0}{a_\mathrm{d}}\right)^2\left(\frac{a_\mathrm{d}}{a}\right)^5.$$
(31)   
Defining $`\alpha =a/a_\mathrm{d}`$, we can rewrite this as

$$\frac{d^2h}{d\alpha ^2}+\left[\frac{5}{\alpha }+\frac{1}{\alpha ^4}\right]\frac{dh}{d\alpha }+\left[\frac{4}{\alpha ^2}+\frac{2}{\alpha ^5}\right]h=\left(\frac{3}{2}H_0\alpha _0^2\right)\frac{1}{\alpha ^7}.$$
(32)   
Now we use the following general result : if $`\phi `$ is a solution of

$`y^{\prime \prime }+f(x)y^{}+g(x)y=k(x)`$

when $`k=0`$, then the general solution is

$`y=C_1\phi +C_2\phi {\displaystyle \frac{dx}{\phi ^2E}}+\phi {\displaystyle \frac{1}{\phi ^2E}\left(\phi Ek𝑑x\right)𝑑x},`$

where $`E=\mathrm{exp}f𝑑x`$. By inspection, a solution of the homogeneous equation (32) is $`1/\alpha ^2`$. It follows that the general solution is

$$h(a)=H_0\left(\frac{a_0}{a}\right)^2\left\{c_1+c_2\mathrm{Ei}\left[\frac{1}{3}\left(\frac{a_\mathrm{d}}{a}\right)^3\right]+\frac{3}{2}\mathrm{ln}\left(\frac{a}{a_\mathrm{d}}\right)\right\},$$
(33)   
where $`c_1`$ and $`c_2`$ are arbitrary integration constants and Ei is the exponential-integral function 

$`\mathrm{Ei}(x){\displaystyle _{\mathrm{}}^x}{\displaystyle \frac{e^v}{v}}𝑑v=𝒞+\mathrm{ln}x+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^k}{k!k}},`$

with $`𝒞`$ denoting Euler’s constant.

By equations (22) and (33), the bulk stress to first order is

$$\mathrm{\Pi }=(3\overline{H}^24\overline{H}h2h^{}Ha)\mathrm{\Gamma },$$
(34)   
This expression holds for $`a>a_0`$, where $`a_0`$ marks the onset of dissipative evolution. Thereafter, the bulk stress decays according to the causal law (19). In order to relate the constants $`c_1`$ and $`c_2`$, we require, according to standard matching conditions, that $`H`$ is continuous. Thus $`h(a_0)=0`$, which fixes $`c_1`$:

$$c_1=c_2\mathrm{Ei}\left[\frac{1}{3}\left(\frac{a_\mathrm{d}}{a_0}\right)^3\right]\frac{3}{2}\mathrm{ln}\left(\frac{a_0}{a_\mathrm{d}}\right).$$
(35)   
Thus, using Eq. (33), we see that the backreaction of the dissipative decoupling process on the expansion of the universe is given by

$`\delta H`$ $`=`$ $`\overline{H}\left\{c_2\left(\mathrm{Ei}\left[{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{a_\mathrm{d}}{a}}\right)^3\right]\mathrm{Ei}\left[{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{a_\mathrm{d}}{a_0}}\right)^3\right]\right)+{\displaystyle \frac{3}{2}}\mathrm{ln}\left({\displaystyle \frac{a}{a_0}}\right)\right\}\mathrm{\Gamma }`$ (37)   
$`+O(\mathrm{\Gamma }^2).`$   
Substituting Eq. (35) into Eq. (34), we find that the bulk stress becomes

$$\mathrm{\Pi }=\overline{\rho }\left\{2c_2\mathrm{exp}\left[\frac{1}{3}\left(\frac{a_\mathrm{d}}{a}\right)^3\right]\right\}\mathrm{\Gamma }+O(\mathrm{\Gamma }^2),$$
(38)   
where $`\overline{\rho }=3\overline{H}^2`$ is the equilibrium energy density. Since $`\mathrm{\Pi }<0`$, we require $`c_2<0`$. Below we find a prescription for $`c_2`$ in terms of physical parameters.

## IV Conclusion

In order to complete the model, we need to determine the remaining arbitrary constant $`c_2`$ in terms of physical parameters. A rough estimate, which is consistent with the simplicity of the model, arises as follows. We estimate the duration of the dissipative process as

$$\mathrm{\Delta }aa_\mathrm{d}a_0,$$
(39)   
i.e. we assume that the process ends at $`a_\mathrm{d}`$. Then by Eqs. (8) and (13), the fractional viscous rise in temperature due to decoupling is approximately

$$\frac{\mathrm{\Delta }T}{T}\frac{\mathrm{\Pi }(a_0)}{\overline{\rho }(a_0)}\frac{\mathrm{\Delta }a}{a_0}.$$
(40)   
We can consider the fractional temperature increase as an input from previous kinetic-theory investigations (as described in the introduction), which typically predict it to be $`O(10^3)`$. Note that this small temperature increase is due to dissipative heating, and is not to be confused with the larger temperature increase arising from electron-positron annihilation, which occurs after neutrino decoupling. Our model does not consider the annihilation process.  Then equations (38)–(40) and (34) allow us to estimate the constant $`c_2`$ in terms of the physical parameters $`a_\mathrm{d}/a_0`$, $`\mathrm{\Delta }T/T`$ and $`\mathrm{\Gamma }`$ as:

$$c_2\mathrm{\Gamma }\frac{1}{2}\frac{\mathrm{\Delta }T}{T}\left\{\frac{\mathrm{exp}\left[\frac{1}{3}\left(\frac{a_d}{a_0}\right)^3\right]}{\left(\frac{a_d}{a_0}\right)1}\right\}.$$
(41)   
Finally, we can also estimate the entropy production due to decoupling. By Eqs. (4) and (40), the viscous increase in entropy per particle is approximately

$$\mathrm{\Delta }s3\frac{\mathrm{\Delta }T}{T}.$$
(42)   
Our model describes the response of the cosmic fluid to a bulk stress, which is a very simple thermo-hydrodynamic approximation to more realistic kinetic theory models of neutrino decoupling, but which nevertheless accommodates the dissipative effects and respects relativistic causality. The simplicity of our model allows us to derive analytic forms for the dynamical quantities and the backreaction effects, but it does not incorporate a mechanism for bringing the dissipative process to an end.

Acknowledgements:

This work was partially supported by a European Science Exchange Programme grant.

## A Characteristic velocities for bulk viscous perturbations

Following , we derive equation (6) for the dissipative contribution to the sound speed. The full analysis of the causality and stability of the Israel-Stewart theory was performed in a series of papers by Hiscock and Salmonson . They showed that both issues are closely related and obtained general expressions for the characteristic velocities for dissipative perturbations. Here we extract from their general expressions specific results for the case in which only bulk viscosity is present.

The purely bulk viscous case stems from the general expressions of  by setting all the coefficients coupled to heat flux and shear viscosity to zero. This yields for the speed of propagating transverse modes

$`v__T^2={\displaystyle \frac{(\rho +p)\alpha _1^2+2\alpha _1+\beta _1}{2\beta _2[\beta _1(\rho +p)1]}}0,`$

which is what one expects for scalar sound-wave perturbations. Equation (128) of  governing the speed $`v=v__L`$ of propagating longitudinal modes becomes, on dividing by $`\beta _0\beta _2`$ and setting $`\alpha _0=\alpha _1=0`$,

$`[\beta _1(\rho +p)1]v^4`$   
$`+\left[{\displaystyle \frac{2n}{T}}\left({\displaystyle \frac{T}{n}}\right)_s{\displaystyle \frac{(\rho +p)}{nT^2}}\left({\displaystyle \frac{T}{s}}\right)_n\beta _1\left\{(\rho +p)\left({\displaystyle \frac{p}{\rho }}\right)_s+{\displaystyle \frac{1}{\beta _0}}\right\}\right]v^2`$   
$`+{\displaystyle \frac{1}{nT^2}}\left({\displaystyle \frac{T}{s}}\right)_n\left[(\rho +p)\left({\displaystyle \frac{p}{\rho }}\right)_s+{\displaystyle \frac{1}{\beta _0}}\right]\left[{\displaystyle \frac{n}{T}}\left({\displaystyle \frac{T}{n}}\right)_s\right]^2=0.`$   
Dividing by $`\beta _1`$ and taking $`\beta _1\mathrm{}`$, we have

$$v^2=\left(\frac{p}{\rho }\right)_s+\frac{1}{(\rho +p)\beta _0}.$$
(A1)   
The first term on the right is the adiabatic contribution $`c_\mathrm{s}^2`$ to $`v^2`$, and the second term is the dissipative contribution $`c_\mathrm{b}^2`$, which, requiring $`v^21`$, leads to

$$c_\mathrm{b}^2\frac{\zeta }{(\rho +p)\tau }1c_\mathrm{s}^2.$$
(A2)   
We also learn from  that causality and stability require

$$\mathrm{\Omega }_3(\lambda )(\rho +p)\left\{1\lambda ^2\left[\left(\frac{p}{\rho }\right)_s+\frac{1}{(\rho +p)\beta _0}\right]\right\}0,$$
(A3)   
for all $`\lambda `$ such that $`0\lambda 1`$. This condition is seen to be hold on account of the inequality (A2).

The expression for $`c_\mathrm{b}`$ refines and corrects the statement in  (the first paper to apply causal bulk viscosity in cosmology) that $`\zeta /\rho \tau =1`$ is required by causality.