Microphysical manifestations of viscosity and consequences for anisotropies in the very early universe

It has been known that a non-perfect fluid that accounts for dissipative viscous effects can evade a highly anisotropic chaotic mixmaster approach to a singularity. Viscosity is often simply parameterised in this context, so it remains unclear whether isotropisation can really occur in physically motivated contexts. We present a few examples of microphysical manifestations of viscosity in fluids that interact either gravitationally or, for a scalar field for instance, through a self-coupling term in the potential. In each case, we derive the viscosity coefficient and comment on the applicability of the approximations involved when dealing with dissipative non-perfect fluids. Upon embedding the fluids in a cosmological context, we then show the extent to which these models allow for isotropisation of the universe in the approach to a singularity. We first do this in the context of expansion anisotropy only, i.e., in the case of a Bianchi type-I universe. We then include anisotropic 3-curvature modelled by the Bianchi type-IX metric. It is found that a self-interacting scalar field at finite temperature allows for efficient isotropisation, whether in a Bianchi type-I or type-IX spacetime, although the model is not tractable all the way to a singularity. Mixmaster chaotic behaviour, which is well known to arise in anisotropic models including anisotropic 3-curvature, is found to be suppressed in the latter case as well. We find that the only model permitting an isotropic singularity is that of a dense gas of black holes.

It has been known that a non-perfect fluid that accounts for dissipative viscous effects can evade a highly anisotropic chaotic mixmaster approach to a singularity. Viscosity is often simply parameterised in this context, so it remains unclear whether isotropisation can really occur in physically motivated contexts. We present a few examples of microphysical manifestations of viscosity in fluids that interact either gravitationally or, for a scalar field for instance, through a self-coupling term in the potential. In each case, we derive the viscosity coefficient and comment on the applicability of the approximations involved when dealing with dissipative non-perfect fluids. Upon embedding the fluids in a cosmological context, we then show the extent to which these models allow for isotropisation of the universe in the approach to a singularity. We first do this in the context of expansion anisotropy only, i.e., in the case of a Bianchi type-I universe. We then include anisotropic 3-curvature modelled by the Bianchi type-IX metric. It is found that a self-interacting scalar field at finite temperature allows for efficient isotropisation, whether in a Bianchi type-I or type-IX spacetime, although the model is not tractable all the way to a singularity. Mixmaster chaotic behaviour, which is well known to arise in anisotropic models including anisotropic 3-curvature, is found to be suppressed in the latter case as well. We find that the only model permitting an isotropic singularity is that of a dense gas of black holes.

I. INTRODUCTION
While the currently observable universe is isotropic to a very high degree, this is not a generic feature of spacetimes near singularities -rather the opposite. In fact, within general relativity and under certain assumptions in the matter sector, the most generic approaches to singularities (such as cosmological big bang or big crunch singularities, but also black hole singularities) are highly anisotropic. They display a chaotic mixmaster behaviour [1], with infinite chaotic oscillations on a finite time interval. This behaviour proceeds in epochs, the beginning and end of which are well approximated by the Kasner metric [2]. This is known as the Belinski-Khalatnikov-Lifshitz (BKL) singularity [3] (see, e.g., [4][5][6] for reviews).
In the context of the current paradigm of very early universe cosmology, anisotropies from the big bang initial singularity would be quickly washed away by a period of accelerated expansion (i.e. inflation). However, our knowledge of possible pre-inflationary physics is very scarce, and the question of what happened near the big bang remains of fundamental interest. In particular, semi-classical general relativity most likely does not hold anymore at such high energy scales, and what was the nature of the initial big bang singularity (if there was one) remains an open question, especially whether it was of BKL type or rather isotropic. Within string theory, chaotic anisotropies are expected (e.g., [7][8][9][10]), while * c.ganguly@damtp.cam.ac.uk † jerome.quintin@aei.mpg.de some semi-classical higher-derivative theories of gravity have stable isotropic cosmological singularities (e.g., [11]) or can limit the growth of anisotropies [12,13]. Other theories of modified gravity can similarly bound shear anisotropies [14] or screen them [15,16]. Moreover, in a gravitational ultraviolet-complete theory such as quadratic gravity [17], requiring semi-classical cosmological transition amplitudes from the big bang to today to be well defined and finite severely constrains anisotropic singularities [18,19].
In the context of alternative very early universe scenarios such as models of bouncing cosmology, the question of the evolution of anisotropies also plays an important role, even well before the approach to the high-curvature big crunch/bounce singularity (or before a non-singular bounce occurs). For instance, in a Bianchi type-I universe, the contribution from shear anisotropies to the total energy density budget is proportional to 1/a 6 , where a is the spatially averaged scale factor. While this decays very rapidly in an expanding universe, it conversely grows much faster than for other known matter types (e.g., pressureless dust and radiation), thus representing an instability to standard isotropic contracting models. This is not an issue for ekpyrotic cosmology since the isotropic background scaling solution arises from a scalar field with energy density growing as 1/a 2 with > 3, thus effectively diluting anisotropies. As such, the ekpyrotic scenario has been shown to be very robust with respect to dynamically producing an isotropic universe, as demonstrated by analytic and numerical studies [20][21][22][23][24][25]. In fact, formally, ekpyrosis implies a no-hair theorem stating that the future big crunch is a stable isotropic arXiv:2109.11701v3 [gr-qc] 31 Jan 2022 singularity [26]. The theorem does not hold, however, if instead of an ekpyrotic scalar field one has an imperfect fluid with anisotropic pressures that satisfy ekpyrotic equations of state [27] (deviations from perfect fluids will be further discussed below). It is also to be noted that single-field ekpyrosis predicts a blue spectrum of scalar perturbations. This is resolved in the case of two-field ekpyrosis at the cost of introducing an additional degree of freedom.
In the context of matter bounce cosmology, where a scale-invariant power spectrum of adiabatic curvature perturbations is generated during a matter-dominated contracting phase [28][29][30], the growth of anisotropies represents a serious problem [31] (to the opposite of ekpyrotic cosmology), which prevents the simplest models from being viable. Of the very few possible resolutions to this problem, we can mention the hypothetical possibility of promoting the graviton to a massive spin-2 field with mass larger than the Hubble scale in the contracting phase [32]. Therefore, in matter bounce cosmology and in a more general context of a bouncing or cyclic universe (not ekpyrotic), the question of how could the universe become isotropic enough for some structure formation scenario to successfully work and/or for a non-singular bounce to be achieved 1 remains mostly unsolved.
Most approaches to cosmology from an effective field theory point of view often assume the matter content to be represented by minimally coupled scalar fields or perfect fluids. However, non-viscous fluids are an approximation to more realistic fluid dynamics models. For instance, scalar fields non-minimally coupled to gravity [38,39], neutrinos that are free streaming [40][41][42][43][44][45][46], or any realistic interacting fluid all depict some form of viscosity. Therefore, the influence of viscosity on early-and late-time cosmologies has been studied in the context of Refs.  and many more. In particular, non-singular solutions have been found in the context of bulk viscosity (see, e.g., [54]), and the effect of shear viscosity has been studied as an isotropisation mechanism [4,[62][63][64]. However, the formulation of the shear viscosity term in analytic form -while accounting for relativistic effects -is challenging. Eckart [65] and Landau-Lifshitz [66] formulate a hydrodynamic relativistic theory of shear viscosity for models whose characteristic motion timescales are much larger than the relaxation time of the system to equilibrium. Close to a singularity, most characteristic motion should cease, so this approximation would not 1 Most realisations of a non-singular bounce usually simply rely on the presumption of isotropy. However, non-singular bounces with sizable anisotropies are possible (see, e.g., [33][34][35][36]), but anisotropies should be at least small enough after the bounce at the onset of radiation-dominated expansion to match later observational constraints from the cosmic microwave background [37]. Furthermore, if we include curvature in our bouncing model, significant anisotropies close to the bounce may not allow the universe to re-expand and create a singularity in the Weyl curvature tensor.
apply. This situation is applicable to the case of a contracting universe close to a bounce when the anisotropy energy density would grow the fastest in the absence of any other isotropising mechanism. Moreover, the theory formulated by Eckart and by Landau-Lifshitz allows for the superluminal propagation of viscous excitations. The Israel-Stewart [67] theory is able to rid the formalism of this problem. For the purposes of this work, we will be using a restricted version of the Israel-Stewart formalism to model the shear viscosity term. The restriction will apply in that we assume that the relaxation time to equilibrium for the fluid under consideration is very small. There exists no closed form for the viscosity term for non-zero relaxation times.
With these considerations, one can arrive at a phenomenological model for the coefficient of viscosity η as a power law of the energy density ρ, i.e., η ∝ ρ n . This was the approach of previous studies, e.g., [4,[62][63][64] in the context of approaches to singularities 2 , but this still remains a phenomenological model. There is no microscopic model -analogous to the kinetic theory picture of colliding hard spheres -of the origin of this viscosity for a cosmological model. It is our intention in this work to provide the beginnings of such a microscopic realisation for viscosity embedded in concrete cosmological scenarios.
Two main avenues will be explored: an interacting scalar field in a thermal bath and black holes. The former has been extensively studied in quantum field theory (QFT), with sophisticated techniques to compute the viscosity coefficient (see, e.g., [68][69][70]). Also, strongly interacting QFTs often have gravity duals (in a holographic description), from which computations have led to a viscosity bound conjecture (see, e.g., [71][72][73]). This conjecture implies that realistic, interacting fluids always have a minimal amount of viscosity, at least of the order of their entropy density. All of this motivates us to consider a simple QFT in a cosmological background as a first microphysical realisation of viscous cosmology.
The second avenue involves black holes, which are often ubiquitous in cosmological scenarios involving a phase of contraction prior to a bounce. Indeed, black holes could form from direct collapse of inhomogeneities [74][75][76] or already exist from preexisting structures (as in a cyclic universe). Such black holes are expected to potentially dominate the universe near a big crunch or bounce (except possibly in regions which could undergo ekpyrotic contraction [77,78]), and as such, a dense 'gas' of black holes has been proposed as a state of matter at very high energies, as studied in string theory (see, e.g., [79][80][81][82], as well as [83][84][85][86] for the related holographic scenario and [87,88] for string-size black holes). In the context of black holes forming in a contracting universe, there is a serious possibility that such black holes could persist through a bounce, thus transitioning into our expanding universe as primordial black holes [89][90][91][92] or remnants thereof [93][94][95]. The important novelty of this work is in realising that black holes, due to their gravitational attraction and intrinsic non-deformability [96][97][98], can be treated collectively as a non-perfect fluid with shear viscosity. Therefore, under certain approximations where the hydrodynamical approximation is valid, dissipative effects form, which tend to isotropise the cosmology.
Outline We shall begin in Sec. II by reviewing the concepts of stress, shear, viscosity, and their phenomenological implications for anisotropic cosmologies, with an emphasis on the models that are later studied in this paper. We then demonstrate in Sec. III some microphysical examples of shear viscosity: the case of an interacting scalar field theory at finite temperature and a gravitationally interacting gas of black holes, both in its dilute and dense limit. We study the effect of the viscosity coefficients derived in these scenarios on the small and the large anisotropy limits of the background universe in Sec. IV. We briefly comment on the implications for gravitational waves in Sec. V, and finally in Sec. VI, we present our conclusions.
Notation Throughout this paper, we use the mostly plus metric signature (−, +, +, +). Latin indices at the beginning of the alphabet run over spacetime coordinates (a, b, c, d, . . . ∈ {0, . . . , 3}), while Latin indices from roughly the third of the alphabet run over spatial coordinates only (i, j, k, . . . ∈ {1, 2, 3}). We also work with units where the speed of light, the reduced Planck constant, and the Boltzmann constant are set to unity (c = = k B = 1), and M 2 Pl := 1/(8πG N ) defines the reduced Planck mass in terms of the Newtonian constant of gravitation G N .

II. REVIEW OF STRESS, SHEAR AND VISCOSITY
A. The definition of the shear and stress-energy tensors and the meaning of viscosity Spatially homogeneous, anisotropic models can be investigated using the orthonormal frame formalism from dynamical systems analysis (see, e.g., [99,100]). The geometry is split into a fluid moving orthogonally to the homogeneous spatial hypersurface, with the timelike fluid 4-velocity u a being equal to the unit normal vector of the spatial hypersurface, hence g ab u a u b = −1. In the spirit of the 3 + 1 decomposition of the spacetime manifold, the fluid velocity vector that defines the foliation can be used to find a projection tensor, which represents the induced metric on the spatial hypersurface. The corresponding extrinsic curvature of the spatial hypersurface is then given by where the last equality defines the spatial covariant derivative, i.e., the spacetime covariant derivative projected on the spatial hypersurface. With simple tensorial algebra, the above can be used to show that the extrinsic curvature tensor can also be written as where the time derivative of the fluid velocityu a := u c ∇ c u a defines the acceleration of the fluid. The extrinsic curvature tensor can be decomposed into an expansion tensor and a vorticity tensor as K ab = Θ ab +ω ab , which are respectively symmetric (Θ ab = K (ab) ) and antisymmetric (ω ab = K [ab] ). We assume throughout that the spacetime has no vorticity, so we set ω ab ≡ 0. The expansion tensor can be further decomposed as where Θ := g ab Θ ab = D a u a = ∇ a u a is the trace part known as the expansion scalar, while the traceless part defines the shear tensor σ ab (so g ab σ ab = 0). Gathering the above, the shear tensor can be written fully in terms of the fluid velocity as or alternatively as σ ab = ∇ b u a + u bua − h ab Θ/3 under the no-vorticity assumption. Albeit the shear tensor is traceless, a useful scalar characterisation of the shear is defined as σ 2 := σ ab σ ab /2, and this is used throughout this work. The symmetric energy-momentum (or stress-energy) tensor for a generic fluid can be written as where ρ = u a u b T ab is the energy density, p = h ab T ab /3 is the pressure, π ab = h (a c h b) d T cd − ph ab is the anisotropic stress tensor (the components are also known as the anisotropic pressures), and q a = −h a b u c T bc is the heat conduction vector measured by an observer comoving with the fluid. For the purposes of this work, we ignore heat transfer (q a ≡ 0). The extra term for a nonperfect fluid represented by the anisotropic stress, which will soon be related to shear viscosity, has to satisfy the constraints π ab = π ba , g ab π ab = 0, and u a π ab = 0, by virtue of being the projected (symmetric) traceless part of the energy-momentum tensor. Upon introducing a dissipative term in the energy-momentum tensor such as π ab , one has to be aware that the fluid may deviate from its thermodynamic equilibrium, and the relaxation time τ to the equilibrium state (a.k.a. the collision time or Maxwell time) may generally be non-zero. In such a case, there is no closed analytic expression for these viscous anisotropic pressures. Instead, they are defined via a differential equation [4,101], where η is known as the viscosity coefficient. Then, the entropy density, considering only shear viscous terms and setting the bulk modulus to zero, can be expanded as [51] (see also [100]) where s 0 represents the entropy density at equilibrium and T is the fluid equilibrium temperature. In order for a fluid description of the system to still be valid, we have to assume that the system is fairly close to the equilibrium state. This can thus be quantified by the following approximation, where the approximate equality on the right-hand side precisely holds when τ is small. For our purposes, we neglect τ in comparison to the equilibrium entropy density, and so the differential equation (7) collapses to the simpler expression as can be seen in standard textbooks (e.g., [66,102]). Since π ab is the (projected) symmetric traceless part of T ab , it is natural for it to be proportional to the other projected symmetric traceless tensor defined above, namely the shear tensor. This theory with τ → 0, though standard, may be plagued by a superluminal propagation of shear excitations -the velocity of this propagation is given by [4,102] c s ∼ η ρτ .
This issue is discussed further below. Deviations from a perfect fluid are sometimes written as (ignoring heat transfer) wherep now denotes the perfect fluid pressure or average pressure. The deviation from a perfect fluid can then generally be written as a linear combination of the trace and traceless parts of the expansion tensor as which implies that the energy-momentum tensor becomes and hence the 'total pressure' is p =p − ζΘ. As shown in, e.g., Refs. [44,99,100], the proportionality coefficients η and ζ have the thermodynamical interpretation of shear viscosity and bulk viscosity, respectively. Bulk viscosity has the effect of modifying the pressure term.
In the case of a flat universe, the bulk viscosity term, which is proportional to the volume expansion rate Θ, can be expressed as a non-linear equation of state (EoS) p = p(ρ), with the pressure being a quadratic function of energy density. Quadratic equations of state have been shown to admit non-singular bouncing solutions [63,[103][104][105]. Given the relevance of anisotropic stress and shear anisotropies for this work, we are mostly concerned by the shear viscosity entering in (10), hence we assume no bulk viscosity throughout (ζ ≡ 0, so the 'total pressure' and 'perfect fluid pressure' have the same meaning, i.e.,p = p). The resulting energy-momentum tensor, T ab = ρu a u b + ph ab − 2ησ ab , is the same as motivated in the previous paragraph.
To gain some intuition about shear viscosity, let us consider a Minkowski background for the time being. We can do this without loss of generality to derive the viscosity coefficient since it is an intrinsic property of the fluid (just like an EoS). In other words, by the equivalence principle of general relativity, we are free to consider a locally Minkowski space to derive the properties of the fluid and later apply such properties to a curved spacetime. Equations (5) and (10) for a Minkowski metric tell us that where we are specialising ourselves to the spatial components. 3 Let us further consider a simplified setup in (2 + 1) dimensions, where one has a fluid in between two infinite-dimensional plates. Let the fluid move in the +x direction (in Cartesian coordinates), with a velocity that only depends on the y direction, i.e., u i (x, y, z) = (u x (y), 0, 0). This is the typical setup to derive the heuristic expression for a fluid's viscosity from kinetic theory first principles (see, e.g., [106][107][108]). In this setup, it is clear that the above relation between stress, viscosity and shear reduces to hence the viscosity coefficient is the proportionality factor that relates the net momentum flux through a constant-y surface to the velocity gradient of the fluid. In other words, viscosity is a measure of the rate of momentum diffusion in the fluid. The mean distance between interactions among the fluid's microscopic constituents is characterised by the mean free path mfp , and thus, the velocity difference of particles moving through a constant-y surface is proportional to − mfp ∂ y u x . This assumes that the mean free path is much smaller than the overall size of the system, i.e., mfp L, where L is the distance between the two plates in this simplified setup. The momentum flux is then proportional to multiplying − mfp ∂ y u x by the energy density ρ and the mean propagation speed of the particles or the root-mean-square speed for a given statistical distribution; for the purpose of our work, as an approximation, we will simply associate this speed with the sound speed c s . Combining the above, we arrive at the expression where α is a proportionality constant of O(1) whose precise value depends on the exact microphysics at play, the statistical distribution, etc.; moreover, it shall encapsulate the uncertainty in our choice of mean velocity. In the second equality above (up to a proportionality factor of order unity, hence the sign ∼), we used the fact that we can write the mean free path as mfp = 1/(βnσ cs ) in terms of the number density n, related to the energy density by the energy of the individual particles E via ρ = En, and the cross sectional area σ cs , which is a measure of the interaction probability among particles. The constant proportionality factor β of order unity again depends on the exact statistical distribution. From the above, we see that the smaller the interaction probability (i.e., the smaller the cross section), the farther a particle travels before interacting with another one (i.e., the larger the mean free path), the easier the momentum transfer, and therefore the larger the viscosity is. However, one needs to be careful since it would appear the limit σ cs → 0 implies infinite viscosity, when one would rather believe that a fluid with no interactions should be viscous-free. Indeed, the issue with the vanishing cross section limit is that it implies an infinite mean free path, hence the assumption mfp L is broken. In cosmology, the size of the system of interest can be associated with the Hubble radius, L ∼ |H| −1 . We shall thus be particularly careful with this assumption throughout this work in order to remain in the regime of validity for the expression (17) to hold. Nevertheless, if mfp ∼ L, it does not mean that viscosity goes away. In fact, the approximation can often be pushed to that limit within order 1 corrections that slightly reduce the viscosity (see, e.g., [106]). However, when mfp L, the above expression for viscosity definitely breaks down, and one generally expects η → 0 as mfp → ∞.
Let us mention that the mean free path and the relaxation time (the average time between collisions) are related by the average velocity: mfp ∼ c s τ . Hence, one can see that (11) and (17) are consistently related. This allows us to re-express the approximation (9) as an upper bound on the mean free path, mfp T s 0 ρσ 2 =: max .
Moreover, one could demand the speed of propagation not to surpass the speed of light, which from (11) amounts to a lower bound on the mean free path. Combining those, and from the discussion of the previous paragraph, we arrive at the following regime of validity: Therefore, given a model for which one can compute the viscosity thanks to Eq. (17), the above lower and upper bounds essentially tell us the regime of validity of that expression in terms of the size of the fluid's mean free path. This shall be the basis of our consistency checks throughout this work.

B. The effect of shear viscosity in anisotropic cosmology
In order to study the effect of shear viscosity of the form (10) in cosmology, let us write down the Einstein equations with no cosmological constant, G ab = M −2 Pl T ab , as follows when q a = ω ab =u a = 0 [100], where (3) R ab and (3) R are, respectively, the 3-dimensional Ricci curvature tensor and scalar on the spatial hypersurface. This system has within it cosmologies containing anisotropies in the expansion, i.e. different expansion rates in the 3 different spatial directions, as well as containing anisotropies in the 3-curvature.
In order to study the effects of anisotropic pressure on an anisotropic universe, let us specialise to a simple flat anisotropic universe. This is known as the Bianchi type-I universe. It represents the case of maximal anisotropy when it is empty, in which case it is called the Kasner solution. It also only has expansion anisotropy, instead of anisotropy in the 3-curvature as well. The metric can be represented as with the constraint denotes the anisotropy in direction x i and a(t) denotes the spatially averaged scale factor (in the sense that ln a = ln a (i) with a (i) = ae β (i) ). From this, H(t) :=ȧ/a defines the spatially averaged Hubble parameter, and the hypersurface geometry is characterized by (3) In particular, 2 (i) . The resulting equations of motion (EOMs) are where ρ σ := M 2 Pl σ 2 defines the energy density in shear anisotropies.
In the presence of a perfect fluid, the stress tensor vanishes, and we recover the shear equation 4 and so ρ σ ∝ 1/a 6 , according to which shear anisotropies essentially contribute to the Friedmann equations as a perfect fluid with stiff EoS p σ = ρ σ . In particular, one recovers the result that anisotropies typically dominate the energy budget of the universe near cosmological singularities since, as a → 0, ρ σ ∝ 1/a 6 grows faster than ρ ∝ 1/a −3(1+w) for a background perfect fluid with EoS w := p/ρ ∈ (−1, 1). As a result, the spacetime near the singularity is well approximated by the anisotropic Kasner metric, and the approach to the metric is of BKL type, as mentioned in the Introduction. An immediate loophole is if the matter EoS satisfies w > 1, which is known as an ultra-stiff ekpyrotic EoS, in which case the background energy density dominates as the scale factor goes to small values, hence isotropising the universe such that it becomes well approximated by a Friedmann-Lemaître-Robertson-Walker (FLRW) metric. In the same situation, the general existence of anisotropic stresses acts as a positive source, and the shear equation of motion is modified according to Eq. (23d). This causes the energy density in the anisotropies to grow faster than a −6 , and hence an ekpyrotic fluid can no longer be reliably expected to isotropise the universe. On doing an extension of this study to anisotropic cosmologies with anisotropic 3-curvature as well as expansion anisotropies (for example, in Bianchi type IX), one finds that anisotropic stresses even if they are ultra-stiff on average, fail to isotropise the cosmology. In fact, a bounce fails to occur as the geometry approaches an anisotropic singularity [27].
Let us now try to gain some intuition about how shear viscosity might change this picture. This was discussed initially in the context of neutrino viscosity and its effects on isotropisation [41]. The discussion was extended to derive a possible phenomenological form of such a shear viscous term in [4]. There, one postulates a shear viscosity coefficient in an anisotropic stress of the form (10), which is dependent on a power law of the energy density as η ∝ ρ 1/2 . The power is 1/2, which allows for isotropisation and an attractor behaviour to a Friedmann singularity [4]. This form of the viscous anisotropic stress, though, allows for the propagation of super-luminal excitations, which we have at the cost of the viscous stresses having a closed form.
If the shear viscosity enters the stress tensor as in Eq. (10), then the matter conservation equation and the shear EOM are generally modified as follows: Together with Eqs. (23a)-(23b), those are typically coupled, first-order ordinary differential equations (not necessarily linear), for which analytic solutions can be found only in special cases. Moreover, the viscosity coefficient is in general time dependent (i.e., it may depend on background quantities such as a, ρ, H, etc.). Let us first consider the simplest case where it is simply a constant, i.e., η = constant =: κ, with mass dimension 3. We use a different variable κ here to denote the constant viscosity since we will use such a positive, dimensionful 5 constant of proportionality for the viscosity coefficient throughout, i.e., it will serve as a reference scale in the time-dependent examples below. Accordingly, the shear EOM becomes whose general solution can be written in the form as was already found by Misner [40]. We notice that the 1/a 6 behaviour is modified due to the constant viscosity coefficient κ by an exponential factor in time. Applying this to our physical considerations of interest, we note that the BKL approach to a singularity 6 is probably not affected too much since, as a → 0 and t → 0, one gets very close to the situation where ρ σ ∼ 1/a 6 → ∞. Nevertheless, the exact growth rate of the shear anisotropies is modified in the approach to a singularity, but its exact value can only be recovered provided a solution for a(t) is also found, which requires additional input. While a constant viscosity coefficient may not isotropise a singularity, it might still dilute anisotropies over an intermediate timescale thanks to the above exponential suppression in ρ σ . Such an example will be explored in greater detail in the subsequent section. As a second example, let us explore the possibility that η = κ/a 3 , which we will motivate in the next section. In fact, this will appear as a possible scaling of the viscosity coefficient in the context of an interacting scalar field theory in a radiation bath. In such a context, the shear EOM can be rewritten as where a prime here denotes a derivative with respect to a.
As a result, one can solve the above differential equation, and the evolution of shear anisotropies is modified as Interestingly, if t 0 < 0 (contraction), one finds that ρ σ → 0 as a → 0, and so it appears that anisotropies have been fully washed out by the time of a big crunch. Alternatively, if t 0 > 0 (expansion), one finds that ρ σ badly blows up in the backward approach to the big bang, exponentially more severely than in the BKL case. Equivalently, it means the anisotropies very quickly decay under forward time evolution in an expanding universe. However, in both instances (contraction and expansion), the meaning of viscosity near the singularity might be lost, as will be discussed in greater detail in the next section. As a last example for this section, let us consider the possibility that η = κ|H| (in this case, κ has mass dimension 2), which will be further motivated in the next section. 7 For simplicity, let us rewrite this expression as η = εκH with ε = +1 for H > 0 (expansion) and ε = −1 for H < 0 (contraction). The shear EOM in this case reduces to whose solution is immediately read off to be Interestingly, the growth rate of the shear anisotropies is modified in such a case by adding a correction to the power of the 1/a 6 scaling; in fact, one can write it as ρ σ ∝ 1/a 6+δ with δ = 4εκ/M 2 Pl . One then notices that 7 We note that this corresponds to the case η ∼ √ ρ in a regime where ρ ρσ according to Eq. (23a). This is the scaling that was noticed to lead to perfect isotropisation [4,64].
for ε = −1 (contraction), the shear anisotropies grow less fast than the typical behaviour in the approach to a big crunch 8 , while for ε = +1 (expansion), they grow faster in the (backward) approach to the big bang; they also correspondingly decay faster under forward time evolution out of the big bang.
This analysis can be extended to cases where there is both expansion and curvature anisotropy. One such example is the Bianchi type-IX universe. This is the case that is taken to be the generic approach to the singularity according to the BKL analysis [3]. Due to the presence of the anisotropic 3-curvature, this cosmology on contraction shows infinite chaotic mixmaster oscillations on a finite time interval (when the lower limit of that time interval is 0), which is an attractor behaviour. An isotropisation mechanism that is successful would be able to resolve this attractor behaviour in the form of chaotic oscillations. In the absence of anisotropic pressures, numerical studies by [21] among others show that ekpyrosis is successful in doing this. In a separate work [63], a viscosity coefficient of the form η = κρ 1/2 for some constant κ of mass dimension 1 is shown to successfully isotropise a Bianchi-IX universe as well as mitigate the mixmaster chaos.

III. SOME EXAMPLES OF MICROPHYSICAL MANIFESTATIONS OF VISCOSITY
A. Interacting scalar field theory at finite temperature Let us consider a finite-temperature scalar field theory with action of the form where m > 0 and 0 < λ 1 are the mass and the selfinteraction coupling constant, respectively. The physical mass of the field is well approximated by m(1 + O(λ)) m at weak coupling (after renormalization, up to radiative corrections and at zero temperature). Denoting µ := m/M Pl , one could imagine having the following hierarchy of scales: 0 < µ µ/λ λ 1. This allows one to distinguish various regimes where the system behaves very differently as a function of the temperature T of the thermal bath [69]. Let us emphasize two such regimes: • When 0 < T /M Pl µ, the system is effectively composed of a non-relativistic, dust-like scalar field. Indeed, the potential is dominated by the zerotemperature mass, i.e., V (φ) m 2 φ 2 /2. In an FLRW background, as long as |H| m (or |H|/M Pl µ in dimensionless units), the field is coherently oscillating with vanishing time-averaged effective pressure, i.e., the EoS is that of dust. If one explores the limit of the universe getting smaller, the Hubble scale and the temperature of the thermal bath both grow as |H| ∼ ρ 1/2 ∼ a −3/2 and T ∝ a −1 , so radiation with ρ ∝ a −4 will quickly become dominant. Nevertheless, the mass term in the Lagrangian remains important in intermediate regimes within µ T /M Pl µ/λ.
• When T /M Pl µ/λ, the mass term becomes negligible, so the potential is dominated by the interaction term, i.e. the λφ 4 term. In this regime, the EoS is that of radiation, p = ρ/3, with energy density growing as ρ ∝ T 4 ∝ a −4 . What is crucial is that in this regime the interactions imply a scattering cross section already at the level of the 2 → 2 tree diagram, and consequently, the fluid should have shear viscosity. We will expand on this below.
In an anisotropic background, the anisotropies would quickly begin to dominate in the approach to a singularity for such a scalar field model, ignoring viscous effects. This is most simply seen in the case of Bianchi I, which contains only expansion anisotropies and in which the energy density in the anisotropies grows as a −6 . However, the thermal bath would remain, and thanks to the rising temperature, the regime T /M Pl µ/λ would be reached. From then on, the presence of the interaction term implies the appearance of viscosity, which may alter the evolution of anisotropies even if those are a priori dominant and large. How efficient viscosity may be at damping the anisotropies will be addressed later.
Let us now consider the scalar field above dominated by its self-interaction potential of the form λφ 4 in a thermal bath with temperature T . (This follows the example discussed in Ref. [73]). Then, this scalar field in a thermal bath behaves like radiation with background energy density, number density, and temperature scaling as ρ ∝ a −4 , n ∝ a −3 , and T ∝ a −1 , respectively; in particular, n ∝ T 3 . Additionally, QFT at finite temperature gives a cross section 9 σ cs ∼ λ 2 /T 2 . Putting everything together, the mean free path is 9 Some intuition for this goes as follows [68]: the typical cross section of a λφ 4 theory goes as σcs ∼ λ 2 /s, where s is the square of the center-of-mass energy here. In the limit of interest, in particular for T m, one can argue that the only relevant energy scale is the temperature, hence s ∼ T 2 and σcs ∼ λ 2 /T 2 .
As viscosity must be measured on scales much larger than the mean free path, it follows that one cannot take the decoupling limit, λ → 0, at which the mean free path goes to infinity. We recall that one should demand mfp |H| −1 on cosmological scales. If we are in a radiationdominated background, we have M 2 Pl H 2 ∼ ρ ∼ a −4 ∼ T 4 , and this would mean T /M Pl λ 2 . We would thus have to be in the regime µ/λ T /M Pl λ 2 λ 1, so one cannot take λ too small for the regime to exist in the first place. Of course, once viscosity is taken into account in the Einstein equations, the background evolution is expected to be modified, and one has to find the proper regime of validity then.
Another aspect that must be considered for the above to hold is thermalization. Indeed, the finite-temperature interaction cross section only holds if the scalar field is in thermal equilibrium with the thermal bath, which is the case as long as the interaction rate Γ = nσ cs v is greater than the Hubble rate |H|. In the high-temperature relativistic limit discussed above, the field is relativistic with unit average velocity v , hence the interaction rate is simply equal to the inverse mean free path, Γ = 1/ mfp . Correspondingly, the requirement for thermal equilibrium is the same as the one for the kinetic theory viscosity derivation to hold, i.e., mfp < |H| −1 , which we discussed above and which we will check explicitly when solving the full set of equations in the next section. Certainly, since Γ ∼ T and |H| ∼ T 2 in a radiation-dominated universe, one does not expect thermal equilibrium to hold up to arbitrarily high energy scales. 10 A caveat to keep in mind, however, is that this only applies as a toy model, where the scalar field φ does not couple to any other fields. In a more realistic context, the physics becomes more complicated regarding thermalization. Indeed, a gauge singlet scalar field can couple to other degrees of freedom, such as the standard-model Higgs, fermions, etc. If so, as φ acquires a large vacuum expectation value (VEV), the standard-model fields would obtain large, VEV-dependent masses, which would suppress the amplitude of the scattering processes, hence delaying chemical equilibrium and thermalization. Such discussions can be found in supersymmetry (e.g., [110] and references therein). In our context, this implies that further analysis is certainly needed.
Since the scalar field in (32) has a canonical kinetic term, its sound speed is unity, and correspondingly, the viscosity coefficient can be evaluated according to (17) as The exact coefficient of proportionality is difficult to estimate, but to leading order in small λ and small m/T , it is expected to be in the O(1 − 10 3 ) regime (see, e.g., [68][69][70]). With this expression in hand, one can then solve for the Einstein equations to determine how the shear viscosity arising from the self-interacting scalar field affects the evolution of the shear anisotropies. Equation (35) suggests η ∝ a −3 , which as we saw in Sec. II B, yields the solution (29) upon assuming a radiation-dominated background, according to which the universe isotropises to the future. In the next section, we will solve the equations in more generality by means of numerical methods. This will also allow us to comment more specifically on the regime of validity (19) over which Eq. (35) applies.
Since the entropy density goes as s ∝ a −3 ∝ T 3 for a radiation bath, we notice that (35) implies For λ at least 1, this means there is a constant lower bound on the ratio of the viscosity coefficient over the entropy density, η/s 1. This is reminiscent of the viscosity bound conjecture (e.g., [71][72][73]), claiming η/s ≥ 1/(4π), which comes from anti-de Sitter black hole solutions in various gravitational theories that are holographically dual to strongly interacting QFTs (nonperfect fluids with shear viscosity). It is thus an interesting observation that many theories suggest a lower bound on the viscosity coefficient, further motivating the investigation of this work.

B. Dilute gas of black holes
In a contracting universe that is not perfectly homogeneous, perturbative inhomogeneities grow under contraction in a similar manner to anisotropies. This growth of inhomogeneities could lead us to an endpoint where a pre-bounce early universe could be populated by a gas of black holes. It was shown in [75,76] that a perfect fluid with quantum vacuum initial conditions in the asymptotic past or thermal initial conditions at a finite time inevitably end up collapsing into Hubble-size black holes at a scale that is determined by the smallness of the fluid's sound speed. If structures already exist when the universe starts contracting (such as in a cyclic context), smaller black holes form first. Initially, these black holes are dilute -this has been modelled as each black hole being situated on a lattice in [90,92]. As the universe contracts, these black holes become denser with a contracting Hubble radius. There is then a dense limit where the Schwarzschild radius R is of the size of the Hubble radius, R ∼ |H| −1 . This case will be dealt with in the next subsection.
In this current subsection, we shall be interested in the possible dilute case where R |H| −1 , far in the contracting phase when the universe is still very large, i.e., very far away from the putative ultimate crunching singularity. We shall model a dilute gas of black holes as a set of hard balls of radius R that nevertheless attract one another gravitationally. To that level of approximation, these could in fact just be any astrophysical objects (small elliptical galaxies, stars, etc.), which might as well populate the universe in this scenario. Interactions between these hard spheres would give rise to a viscous drag, very similar to that derived in kinetic theory. Although the black holes in the dilute gas would have to coalesce to form a gas of larger black holes to ultimately enter into the dense limit R ∼ |H| −1 , effects of black holes (or other astrophysical objects) coalescing is not taken into account in this analysis in the dilute limit. In fact, we do not know exactly when the perturbations become too large as to not trust the approximations, hence we must add a word of caution. While the approximations might hold initially, it is unclear how long they may last, and this has to be taken into account when drawing conclusions.
We start by saying that the cross section for two black holes as described above to interact is given by (see, e.g., [111]) where the sound speed c s of the gas represents the average velocity of the distribution 11 of black holes. We note that in the relativistic limit where c s → 1 the expression for the cross section reduces to that for non-interacting hard spheres, σ cs ∼ R 2 , as expected. Alternatively, in the pressureless limit c s → 0, the cross section tends to infinity. This is understood from the fact that if all the black holes in the gas were perfectly static (say at some initial time), they would inevitably merge in some finite time due to the infinite-range gravitational attraction between them, hence the certain collision probability. Making use of the relation between the Schwarzschild radius and mass (upon specialising our attention to black holes), R ∼ M/M 2 Pl , the number density of the gas is related to its energy density via n ∼ ρ/(RM 2 Pl ), from which we can read the mean free path mfp ∼ (nσ cs ) −1 as The viscosity can then be evaluated as 11 One would generally expect a distribution of masses/radii and velocities for the gas of black holes. Here we are thus referring to R and cs as the mean radius and velocity, respectively. We are not making any assumption about the distribution since too many factors come into play in the formation of such a gas of black holes. Beyond idealised analytical estimates as in [75,76], this would potentially require numerical simulations, which would nevertheless be very dependent on the chosen initial conditions. Therefore, we remain agnostic about exact values for R and cs and treat them as free parameters.
which is just a constant since in the limit R |H| −1 one does not expect the Schwarzschild radius to be affected much by the cosmological background under the present approximations.
Let us comment on the regime of validity of the above expression for viscosity, recalling (19). The inequality η/ρ mfp , which followed from demanding a subluminal propagation speed of shear viscosity excitations, is satisfied provided c s 1. This is not a surprise as we expect the sound speed of the dilute gas of black holes to precisely be the propagation speed of viscosity excitations, and this sound speed is certainly expected to be subluminal.
The inequality on the right-hand side of (19) is less trivial though. To tackle it, let us first make the observation that for a gas of black holes to first form one generally has to be in a background that is relatively close to isotropy. This can certainly be envisioned in the context of a cyclic universe, where a prior expanding phase can efficiently isotropise the universe. Thus, we can assume here that shear is initially subdominant, or at most of the order of the fluid's energy density, i.e., σ 2 ρ/M 2 Pl ∼ H 2 . How the shear subsequently evolves, given the viscosity (39), will be addressed in the following section, but we already saw that a constant viscosity coefficient can lead to temporary exponential suppression of the shear [recall (27)]. Under the assumption that shear is subdominant, the mean free path (38) can be written as mfp ∼ c 4 s /(RH 2 ). The requirement that this is smaller than the Hubble radius thus reads where the right-hand side is expected to be much smaller than unity since we are considering R |H| −1 . Therefore, Eq. (39) for viscosity is expected to apply only if the sound speed is very small. While c s remains at the level of a free parameter given the uncertainties stipulated earlier, one certainly does not expect black holes to have large peculiar velocities upon formation from gravitational collapse, so the above inequality does not appear unreasonable.
The other requirement from the right-hand side inequality of (19), which comes from the small Maxwell time assumption, is generally found to be less stringent than (40) as long as the shear remains subdominant. To see this, let us express the temperature of the dilute black hole gas assuming a Maxwell-Boltzmann distribution of velocities, such that T ∼ c 2 s M . The entropy density can also be read from the sum of the black hole's individual entropies, s ∼ n(RM Pl ) 2 ∼ ρR, which dominates over the 'ideal gas' entropy in this context. Putting those together, we arrive at T s/(ρσ 2 ) ∼ c s RM Pl /σ, and thus the mean free path (38) is smaller than max as long as c 3 s R 2 ρ/(M Pl σ). If σ 2 ρ, this is not a very severe constraint. Even if the shear is of the order of the energy density, then the constraint reduces to c 3 s R 2 |H|M Pl , which is generally no more restrictive than (40) since we expect to be in a deeply sub-Planckian cosmological regime (|H| M Pl ).

C. Dense gas of black holes
A contracting universe that isotropically evolves with a dilute gas of black holes will arrive at a phase where the separation distance is of the order of their Schwarzschild radius. As the black holes are pushed closer together, we arrive at the dense black hole gas picture (e.g., [74]). In this picture, the EoS resembles a stiff fluid p = ρ, a result that is derived from thermodynamic considerations in this current section. In the dense black hole gas picture, every Hubble patch can be thought to be filled by a Hubble-size black hole, i.e., R = |H| −1 . As the universe keeps contracting, one expects some form of quantum instability that allows black holes to 'continuously' bifurcate into smaller black holes such that the relation R = |H| −1 holds as a function of time (this is forbidden classically [112]). Though such a phase is highly hypothetical, it is not violating the second law of thermodynamics as we will see below, and it might well occur if black holes at high densities are to be replaced by stringy counterparts (see, e.g., [79-82, 87, 88]). At the level of semi-classical gravity, such a phase would inevitably still ultimately lead to a collapse of the whole universe into a singularity, but again, this is poorly studied and new physics might well come into play.
Let us consider a region of physical volume V containing N black holes of Schwarzschild radius R, so N ∼ V /R 3 . The total energy in the volume is then given by E ∼ N M ∼ V M 2 Pl /R 2 , where M ∼ M 2 Pl R is the Schwarzschild mass of the black holes. One must keep in mind the following: we assume that we can describe the black holes by their usual Schwarzschild mass and radius coming from the Schwarzschild metric of a single black hole embedded in Minkowski space (i.e., asymptotically flat). This might not hold in a universe that has a possibly infinite number of black holes and that could be dynamical, but we have no good prescription in that situation, so we will stick with the usual Schwarzschild description -more comments are to be given in the discussion section. Then, if the entropy of each black hole is given by the Bekenstein-Hawking entropy, the total entropy in the volume is given by S ∼ N M 2 Pl R 2 ∼ V M 2 Pl /R. These relations can be combined to yield S ∼ M Pl √ EV , or in terms of densities, Using standard thermodynamic relations such as 1/T = ∂ E S and p = T ∂ V S, one finds a temperature T ∼ √ ρ/M Pl ∼ M 2 Pl /M and a pressure p = ρ. It is in that sense that the dense black hole gas picture is akin to a stiff fluid.
To then get the viscosity (which has never been considered before for a dense black hole gas), let us estimate the interaction cross section by σ cs ∼ R 2 , where in analogy with the dilute gas of the previous subsection the propagation speed of fluctuations is essentially taken to be unity for a stiff fluid. The mean free path follows as mfp ∼ 1/(nσ cs ) ∼ R since n = N/V ∼ R −3 . Already, we see that the mean free path is of the order of the Hubble radius by construction. Indeed, if black holes are expected to fill each Hubble patch, then it takes a distance R ∼ |H| −1 before black holes interact with one another. As such, we expect the naïve kinetic expression (17) for viscosity to be only a rough order of magnitude estimate. Nevertheless, it should convey the right scaling as a function of energy density. The above mean free path implies η ∼ ρR, but the energy density is actually related to the black holes' radius as ρ ∼ (M Pl /R) 2 as we saw above, hence we finally obtain It is interesting to notice that, from the results above, the ratio of viscosity to entropy density is constant (and of order unity), as was the case for the finite-temperature interacting scalar field. As already mentioned, we might indeed expect the conjectured bound η/s ≥ 1/(4π) to hold. In fact, by assuming the conjecture, the above result has already been guessed and consequences thereof explored in [81]. Another interesting observation is that (42) implies M 2 Pl H 2 ∼ ρ, which is the Friedmann constraint equation with no anisotropies. This is perhaps not a surprise since the derivation essentially assumes the universe to be isotropic enough for the dense black hole gas to form in the first place. However, it seems to suggest already that no anisotropy is allowed to form when the universe is dominated by such matter. As we saw from Sec. II B, a viscosity coefficient of the form of (42) does indeed lead to isotropisation, i.e., the energy density in anisotropies always remains subdominant compared to the energy density of a stiff fluid (the dense black hole gas in this case). Therefore, this picture of a dense black hole gas represents the only microphysical origin known to the authors of a stiff fluid with viscosity given by the scaling η ∝ ρ 1/2 , which was previously phenomenologically understood to perfectly isotropise the universe [4,62,64], i.e., leading to a Friedmann singularity if taken all the way to a big crunch.

IV. THE EVOLUTION OF ANISOTROPIES IN VARIOUS SCENARIOS
In the previous section, we presented three fluids for which one can derive a viscosity coefficient from a microphysical perspective. The dilute and dense black hole gases can in fact be viewed as a single fluid in opposite limits, while the finite-temperature interacting scalar field is unambiguously different in nature. In deriving the properties of the black hole gas (in both the dilute and dense limits), we had to resort to the assumption that the background cosmology was isotropic to a good approximation in the first place. The question of how anisotropies (even if small initially) can evolve subsequently remains well posed. The goal of this section is thus to explore the evolution of anisotropies for the fluids described in the previous section, first under the assumption of small anisotropies initially (which can apply to both black hole gases and the scalar field example), and then conversely, in the limit of large initial anisotropies (which can only be applied to the scalar field model).

A. Small anisotropy limit
For a black hole gas in the dense limit (where η ∝ ρ 1/2 ), the evolution of anisotropies is already known from the analytical solution (31) and previous works [4,62,64] as already discussed, which confirms the isotropising power of such a fluid. In the dilute limit (where η is constant), we also already obtained an analytical solution in (27), but a solution for the background scalar factor a(t) is needed to fully quantify the evolution of anisotropies in such a case. This is where assuming small anisotropies initially (so approximately FLRW) can be useful analytically.
Under the assumption of a FLRW metric initially, we can solve for the evolution of a Bianchi-I metric for a wide class of viscous fluids. For the sake of generality, let us consider a phenomenological parametrisation of the viscosity coefficient as where the constant n determines how viscosity scales as a function of ρ (e.g., n = 0 for a dilute black hole gas, while n = 1/2 for a dense black hole gas), and κ ≥ 0 is the proportionality factor, whose exact value can be derived from the microphysics of the fluid. In the above, we set up the dimensions such that κ is a dimensionless constant this time. With this parameterisation of viscosity, let us rewrite the background EOMs (23) as follows, where w := p/ρ defines the matter EoS and where we defined to be the ratio of the shear energy density to the matter energy density, which we dub the shear-to-matter ratio. The logic to solve the above analytically shall be the following: consider a contracting universe in which the energy density in anisotropies is initially contributing at most as much as the matter content, i.e., the ratio Ω σ = ρ σ /ρ is at most order 1 initially. Then, one can say that, initially, 3H 2 ρ/M 2 Pl andρ + 3H(ρ + p) 0 (provided κ is also not too large) as a rough approximation. In other words, one assumes that the spacetime is approximately FLRW at the onset of the analysis and check whether that approximation may remain valid under time evolution (i.e., whether it improves or worsens). Practically speaking, this means checking whether or not Ω σ remains ≤ 1. With no viscosity, we already saw that ρ σ grows as a −6 , so even if we start out with small anisotropies, the contraction will cause these anisotropies to grow to such an extent that the universe quickly becomes anisotropy dominated, certainly more than allowed for a successful bounce to occur or for a structure formation scenario to be realised. We are now asking the question whether the inclusion of shear viscosity from a fluid that can reasonably be expected to be present can mitigate the growth of these anisotropies. Our question is thus whether ρ σ may remain subdominant, and in fact, how much it may decay as the universe contracts under the influence of shear viscosity.
The solution to the EOM for σ i j , Eq. (44d), reads where t i is the time at which the initial conditions are set. Given the approximate FLRW background, we have where we shall be looking at the regime where t < 0 for a period of contraction. Performing the integral in (46), it follows that as long as n = 1/2 (one has to treat the n = 1/2 case separately) and where we used the fact that we have A first thing to notice from (48) is that with no viscosity (κ = 0), one is left with Ω σ ∝ |t| −2(1−w)/(1+w) , and therefore, one recovers the usual result that anisotropies grow, are constant, or decay compared to the background energy density as t → 0 − if w < 1, w = 1, or w > 1, respectively (assuming the fluid's EoS parameter is always at least greater than −1). Then, reinserting viscosity with κ > 0, we note that if n > 1/2, the term in the exponential becomes dominated by −(2n − 1)(t i /t) 2n−1 , which goes to −∞ as t → 0 − . Consequently, anisotropies are (exponentially) infinitely suppressed, and the BKL instability is resolved in this regime. However, we do not know of a realistic fluid, which would have a well-defined viscosity all the way to high energy scales with n > 1/2. 12 The case n = 1/2 is treated separately later, so let us focus on the cases when n < 1/2. One can see from (48) that as t → 0 − , the factor in the exponential only goes to a finite negative constant, so while the anisotropies are exponentially suppressed compared to the non-viscous solution, the approach to the singularity remains highly anisotropic for w < 1. The exponential suppression remains interesting though, especially in a context where the universe might not reach a singularity or even Planckian scales. Indeed, it might be interesting to see if there could be significant isotropisation before a bounce occurs. To explore this question, let us first observe that demanding the time derivative of (48) to be negative at the initial time t i , we find that the shear-to-matter ratio Ω σ is initially decaying (demandingΩ σ (t i ) < 0, so the universe is initially isotropising) as long as assuming w > −1, and where H i = 2/(3(1 + w)t i ) is the initial value of the Hubble parameter. What this shows is that, for w = 1, any positive non-zero viscosity coefficient suffices to start isotropisation, i.e., for the shearto-matter ratio Ω σ to start decreasing. When w < 1, however, the viscosity coefficient cannot be arbitrarily small for that matter; it needs to be larger than a minimal value dubbed κ min . The smaller the EoS parameter, the larger κ needs to be. Also, the higher the initial energy scale, the larger the viscosity coefficient must be to be able to begin isotropisation. Vice versa, when the universe is initially very large (small initial energy scale), the viscosity coefficient can be smaller. This dependence on the initial Hubble scale is most important the closer n is to 0, but it becomes less important for n closer to 1/2. Provided isotropisation starts, we can then check at what point there is a turnaround, i.e., a point where the shear-to-matter ratio Ω σ starts growing again. By solv-ingΩ σ = 0, we find that this occurs at an energy scale This expression only applies for w < 1; for w ≥ 1, Ω σ always decreases until a big crunch or a bounce is reached. In fact, the larger w is, the higher the energy scale at which isotropisation stops can be. The same applies for the viscosity coefficient, as expected. Additionally, the closer n is to the value 1/2, the more efficient isotropisation is.
In the cases where isotropisation stops, Ω σ starts growing again, and one can approximate its subsequent evolution by the usual power-law scaling without viscosity, where t = 2/(3(1 + w)H ) is the end-of-isotropisation time following from (50), and Ω σ, := Ω σ (t ) is the corresponding value of the shear-to-matter ratio at that time. The time at which the ratio reaches unity, , represents the moment when anisotropies start dominating again, and so the moment when the initial assumption breaks down and the above solutions do not apply anymore. Past that point, we essentially expect the universe to reach its chaotic mixmaster behavior toward the big crunch.
Let us explore the above timescales in a specific model of interest. Let us consider the case of a constant viscosity coefficient corresponding to n = 0, which was already solved in (27). If we now assume the background to be approximately FLRW with matter having the EoS of dust (w = 0), we can write this sub-case of (48) as This is thus the solution for the evolution of anisotropies in the example of a contracting universe containing a dilute gas of black holes with effective EoS w = 0, which is initially isotropic to a good approximation. The evolution of Ω σ in this case is shown in the top plot of Fig. 1 as a function of the e-folding number defined according to The bottom plot of Fig. 1 shows similar computations, but applying (48) for some phenomenological 13 case with w = 1/12 and n = 1/6. Curves of different color show different values of the viscosity coefficient of proportionality κ, whose value as a fraction of the minimal isotropising coefficient κ min can be read off from the color bar. There, we can see that for κ close to κ min (the curves with lighter color), the universe does start by isotropising, but this is not very efficient, and Ω σ quickly turns over and grows as a power law beyond Ω σ = 1, indicating a future shear-dominated universe. It is only for values of κ that are about 1 or 2 orders of magnitude larger than κ min that we start seeing long-lasting isotropisation (of the order of tens of e-folds). In those cases (darker 13 Such arbitrary values could potentially correspond to some intermediate regime, in between a dilute and a dense black hole gas for instance, or in the case of the scalar field example, in between the matter-and radiation-dominated regimes.  (49). The initial conditions are set at a time ti = −10 80 t Pl , and the initial shear-to-matter ratio is set to the threshold value Ωσ,i = 1. This value is highlighted by the horizontal dotted grey line. curves), isotropisation is extremely efficient (exponential as expected) for the first few e-folds before turnaround and power-law growth. However, since Ω σ shrinks to exponentially small values at first, it takes several tens of e-folds before shear becomes dominant again.
The problem with the above description in the case of the physically motivated dilute black hole gas is that large viscosity coefficients κ/κ min ∼ O(10 2 ) are not expected to respect previously discussed approximations. To see this, let the constant viscosity coefficient η = κM 3 Pl be given according to (39) for a dilute black hole gas. Together with (49) when w = n = 0, we thus find which needs to be at the very least greater than 1 for isotropisation to work, i.e., one needs c 5 s > R|H i |. However, this is clearly incompatible with the requirement that the mean free path has to be smaller than the Hubble radius, cf. (40). Therefore, we conclude that a dilute black hole gas is not viscous enough for isotropisation to start, even less so for an isotropic background to be sustained.
The evolution of anisotropies in the case of an interacting scalar field at high temperature in the small anisotropy limit was already found in Sec. II B. Indeed, assuming the background to be FLRW and radiation dominated and taking the viscosity coefficient to be η ∝ a −3 in accordance with (35), one finds the solution (29), which depicts isotropisation as a → 0. This solution is equivalent to (48) with w = 1/3 and n = 3/4. As mentioned earlier, n > 1/2 immediately implies isotropisation in this limit, but approximations most likely break down before reaching a singularity in this case. For this reason, this requires greater scrutiny, and so we defer the analysis of this scenario to the next subsection, where we look at the large anisotropy limit numerically, making no approximation about the background.
To end this subsection, we come back to the special case of n = 1/2. In such a case, the integral (46) yields Isotropisation thus occurs as t → 0 only if In the case of a stiff fluid with w = 1, we see that any non-zero positive viscosity coefficient of proportionality leads to isotropisation, in accordance with the expectations previously mentioned. We note that for a general EoS such a lower bound on the viscosity coefficient of proportionality has already been derived in [64]. In fact, there it is found that for n = 1/2, whenever the future crunching singularity is a stable Friedmann singularity (i.e., the universe fully isotropises by then). This has been derived for all Bianchi classes, and thus, it may explain the more stringent proportionality factor of 3 compared to √ 3/2 found in our simplified Bianchi-I analysis under the assumption of small anisotropies.
B. Large anisotropy limit

Bianchi I
As mentioned in the previous subsection, the case of an interacting scalar field at finite temperature has strong potential isotropising power, although this remained at the level of assuming small anisotropies initially. The strength of this model, though, lies in the fact that one does not have to make any assumption about the 'formation' of the fluid or its previous history. In other words, even if the universe is highly anisotropic to start with, we would still reasonably expect a λφ 4 scalar field in a thermal bath to exhibit viscosity, and consequently, given its contribution to the coupled Einstein equations, affect the subsequent evolution of the anisotropies. This would even be true if one started in a maximally anisotropic homogeneous flat universe, also known as a Kasner universe. We shall be interested in this initial limit in this subsection, i.e., when anisotropies are dominant over everything else.
Let us first restrict ourselves to the case of flat spatial sections, i.e., to a Bianchi type-I metric (the case with curvature anisotropy, Bianchi IX, is treated separately later). We now seek to solve the corresponding equations (23) numerically, where in the case of the field theory model introduced in Sec. III A, we can use (35) for the viscosity coefficient together with the usual scaling of temperature T ∝ 1/a. In this case, the matter and shear EOMs (25) reduce tȯ assuming the matter EoS w = 1/3 for the radiation bath, and where we denote the (expected order 1) constant of proportionality in the viscosity coefficient (35) by α.
For the purpose of the numerical analysis, we simply set α = 1, and the scalar field self-interaction coupling constant is taken to be λ = 10 −3 . Other numerical values have been explored, but we focus our attention here on the free parameter T 0 , which sets the initial temperature of the thermal bath at the initial scale factor value a 0 . Exploring a range of values for T 0 shall encapsulate different choices for the combination of parameters αT 0 /λ 2 .
We are now in position to numerically solve the set of coupled ordinary differential equations (23b) and (58), which respect the constraint (23a). Solutions are shown in Fig. 2 for the shear-to-matter ratio Ω σ (in the top plot) as a function of the e-folding number N . Initial conditions are picked at the Hubble scale H 0 = −10 −50 M Pl such that the initial shear-to-matter ratio is Ω σ,0 = 10 15 , i.e., we want the anisotropies to be dominant over the matter at the initial time and see how this changes under time evolution. The curves of different color show different values of the initial thermal bath temperature T 0 , ranging from colder (10 −25 M Pl , blue) to warmer (10 −17 M Pl , red).
Starting from the colder temperatures in blue in Fig. 2, we see that Ω σ first grows as the usual power law in a Kasner universe, before starting to saturate. In fact, after about 10 e-folds, Ω σ reaches a constant, already showing that viscosity has started becoming effective in mitigating the otherwise unbounded growth of anisotropies. For the lighter shades of blue and the green/yellow curves , we can see that Ω σ starts by decreasing, demonstrating isotropisation, but the exponential damping does not last. Rather, Ω σ saturates at some constant value greater than 1, meaning that the universe remains anisotropy dominated (though with bounded shear).
The situation changes once we consider initial temperatures warmer than about 10 −19.5 M Pl . For the darker orange curves, we see that there is initial isotropisation followed by saturation, but then there is a second phase of exponential isotropisation, which brings Ω σ to exponentially small values, well below unity, such that the universe is isotropic to a very good approximation. For the red curves, this isotropisation occurs all at once, with no intermediate saturation phase, and the universe becomes isotropic within a few e-folds (or even a fraction of an e-fold for the temperatures closer to 10 −17 M Pl and above).
In all of these cases, it is important to consider whether the viscosity approximation is valid though (and whether thermal equilibrium holds). In its simplest iteration, this can be stated as the situation when the mean free path mfp remains less than the characteristic length scale of the system. In our case, the characteristic length scale is given by the size of the horizon |H| −1 , where H is the average expansion rate as before. For this reason, we show the ratio mfp /|H| −1 in the middle plot of Fig. 2. For our purposes, we use the expression for the mean free path in this field theory, derived in Eq. (34), with the constant of proportionality set to 1. We see that in all cases considered the approximation remains valid for at least 40 e-folds. For the higher initial temperatures that successfully lead to an isotropic universe, we see that the approximation remains valid even longer, up to at least 60 e-folds. Once the mean free path becomes of the order of the Hubble radius and even surpasses it, the expression used for viscosity does not apply anymore. In fact, one would rather expect viscosity to go to zero as thermal equilibrium is lost in the limit where the averaged volume shrinks to zero. Therefore, one cannot realistically expect isotropisation to remain effective all the way to a crunching singularity. Rather, a Kasner singularity is anticipated. Yet, considering the efficiency of the exponential damping of shear within the regime of validity of the theory initially when T 0 is high enough, even if one were to turn off viscosity altogether once mfp ∼ |H| −1 , it would take several hundreds of e-folds (if not more) before the power-law growth in anisotropies would bring Ω σ back to values greater than unity. Therefore, it is expected that in any realistic scenario where the universe would undergo a non-singular bounce at some highcurvature scale such that a singularity is never reached, the universe would still be highly isotropic at the onset of the transition from contraction to expansion.
At last, let us point out that according to (19) another approximation should be satisfied, namely mfp / max 1, where max is defined in Eq. (18). This can be computed using the usual radiation entropy relation s ∝ T 3 . The result is shown in the bottom plot of Fig. 2, where it can be seen that the ratio mfp / max remains well below unity throughout the evolution and for any initial temperature in the given range. In fact, the approximation improves under time evolution and for warmer initial temperatures. Therefore, the requirement that the Maxwell relaxation time be small enough does not represent a threat to the validity of the viscosity approximations in this context.

Bianchi IX
The above discussion only concerns expansion anisotropies, where the underlying geometry is a flat anisotropic universe. However, the generic approach to a singularity and, in our case, the endpoint of contraction is the closed anisotropic universe described by the Bianchi type-IX metric (see, e.g., [113] and references therein). In this case, the anisotropy energy density is not just stored in the expansion tensor, but there is also an anisotropy 'potential'. This potential is nothing but the anisotropic 3-curvature terms that arise in the closed Bianchi type-IX universe, and which are responsible for the chaotic mixmaster oscillations on approach to a singularity.
The Bianchi-IX metric takes the general form of a homogeneous spacetime as follows, Here, h ij is the spatial metric, and the σ i 's are one-forms, which take the simple Cartesian form dx i in the case of a flat anisotropic universe [Bianchi I, cf. (21)]. In the case of homogeneous spacetimes with non-trivial curvature, they can always be chosen so that h ij always remains strictly a function of time. In the case of the Bianchi-IX universe, these one-forms take the following shape, which are the differential forms on a 3-sphere with coordinate ranges 0 ≤ θ ≤ π, 0 ≤ ϕ ≤ 2π, and 0 ≤ ψ ≤ 4π. In the frame in which the metric h ij is diagonal and strictly a function of time, it takes the form The volume averaged expansion is given by the scale factor a(t). The variables β ± (t) are the Misner variables that are used to parameterise the anisotropies. The shear anisotropy has only two independent components as the anisotropic shear is traceless. In this formalism, the three-dimensional curvature (3) R on spatial hypersurfaces of constant coordinate time is given by where the curvature potential U (β + , β − ) is given by The Einstein equations (20) in this formalism become (we set M Pl = 1 for the rest of this subsection) where the shear energy density is σ 2 = 3(β 2 + +β 2 − ). If one ignores the presence of viscosity and simply set η ≡ 0, then one recovers the usual chaotic mixmaster behaviour in the approach to a singularity. To see this, let us numerically solve the above set of ordinary differential equations when the matter content is radiation-like (p = ρ/3). The initial conditions are set in a contracting phase with H 0 = −10 −40 , β +,0 = 10 −10 , β −,0 = −10 −1 , β +,0 = −1, β −,0 = 0, and a 0 ≈ 9.1 × 10 39 , where a prime here denotes a derivative with respect to the efolding number N := − ln(a/a 0 ), which turns out to be an easier time variable 14 to work with in Bianchi IX, numerically speaking. Such values are chosen such that, initially, ρ 0 /(3H 2 0 ) = 1/10, σ 2 0 /(3H 2 0 ) = 1, and − (3) R 0 /(6H 2 0 ) = −1/10. Physically, this means that we choose shear anisotropies to be dominant over matter (radiation) initially since σ 2 0 /ρ 0 = 10, but we want curvature anisotropies to be small. Taking |β ± | 1 initially, the anisotropy potential is negative (the potential minimum is −3/4), indicating positive spatial curvature ( (3) R > 0), and the curvature radius is set to be small by taking a 0 large. In other words, we want to start in a large universe relatively close to a flat Bianchi-I spacetime in this example.
The result of the evolution is shown in the left plot of Fig. 3. There, we see that, without viscosity, the radiation contribution (orange curve) rapidly goes to 0, while anisotropies dominate. In fact, there is a chaotically oscillatory exchange between shear anisotropies (purple curve) and curvature anisotropies (olive curve), representative of the mixmaster dynamics as the universe approaches a BKL singularity.
When viscosity is introduced, the situation changes, as can be seen in the right plot of Fig. 3. There, we numerically solve the same previous set of equations with the same initial conditions, except now the viscosity coefficient is taken to be η = αT 3 0 (a 0 /a) 3 /λ 2 as in the previous subsection. Numerical values for this example are taken to be α = 1, λ = 10 −3 , and T 0 = 10 −16 .
For the first e-fold or so, the evolution with and without viscosity is very similar. However, as the scale factor decreases, the temperature rises and so does the viscosity coefficient. Accordingly, the radiation component is not diluted with respect to the anisotropies; rather, it remains more or less constant and starts growing after a few e-folds. Counterbalancing, the contribution from shear starts decreasing already after about 2 e-folds. By ≈ 8.98 e-folds, radiation becomes dominant over shear and curvature anisotropies (Ω σ becomes smaller than unity). From then on, the spacetime becomes more and more isotropic, with anisotropies decaying to exponentially small values, and with the chaotic oscillations in the anisotropies stopping. Through this evolution, the mean free path mfp is found to remain smaller than the Hubble radius, up to approximately 35.66 e-folds. By then, log 10 Ω σ ≈ −67.71. As discussed in the previous subsection, beyond this point one cannot fully trust the approximations leading to the viscosity coefficient, which should in fact start decreasing. Nevertheless, even if viscosity were to suddenly become negligible again, it would take more than about 78 e-folds before anisotropies would become dominant again. Therefore, we can say that the model is isotropic to a very good approximation for more than 100 e-folds in total, and the warmer the initial temperature T 0 of the thermal bath, the longer the isotropic phase, in the same spirit as seen in Fig. 2 for Bianchi I. In the end, it seems that isotropisation due to viscosity in a finite-temperature field theory is robust against curvature anisotropies, i.e., the same qualitative results hold whether the spacetime is of Bianchi type I or IX.

V. IMPLICATIONS FOR GRAVITATIONAL WAVES
So far, we have discussed the process of isotropisation in a contracting universe. We have shown that the addition of shear viscous anisotropic stress leads to a reduction of the fractional contribution of the shear anisotropies in many instances. The shear anisotropies that we have studied so far have been in spatially homogeneous cosmological settings -they are nonperturbative by definition. Then, the perturbative limit of this represents a homogeneous and isotropic universe, but containing gravitational wave perturbations. The concept of an isotropic universe sourced by gravitational waves being equivalent to an anisotropic universe is not a new one. For example, in an open or flat isotropic Friedmann model, gravitational waves superimposed upon the background only leave homogeneity untouched, and hence reproduce the corresponding spatially homogeneous, anisotropic cosmology when the wavelength is infinitely long [114]. This fulfills the assumption of homogeneity as the periodicity of a propagating wave with finite wavelength would actively violate it. There are some exceptions to this rule, such as in the case of circularly polarised gravitational waves [114], where the average quantities coincide with the scenario of a gravitational tensor representing a homogeneous isotropic cosmological model being sourced by gravitational wave anisotropies. The approach to a singularity also becomes quasi-isotropic and resembles the Friedmann solution. In general, one can assume that whatever physical effect modifies the propagation of shear will also correspondingly affect the propagation of gravitational waves (and vice versa). Examples include massive gravity [32], neutrinos and more (e.g., [46,55,57,58,[115][116][117][118][119][120]).
As there appears to be an inexorable link between the shear anisotropies we have been studying and gravitational waves, it would be interesting to see how a characteristic spectrum of gravitational waves would be affected at the end of shear viscosity driven contraction. For the purposes of this computation, we shall restrict ourselves to a flat background. This indicates the case of Bianchi I. In fact, we can even assume the background to be flat FLRW as the only anisotropies present are in expansion and can be written as part of the energy density. This is a very simple example of the general idea that perturbative shear anisotropies on a homogeneous background can be represented as gravitational waves on an isotropic background.
The general equations of motion in a homogeneous background given in (20), written in terms of the electric and magnetic parts of the Weyl curvature tensor denoted by E ab and H ab , are given by [100] These equations are written assuming linear perturbations around a flat background and follow the covariant and gauge-invariant approach to perturbation theory outlined in [100]. They are easily generalisable to the fully non-linear case as for example in [100,121]. This is different from the metric perturbation approach where we linearise the metric around a background and then trace the time evolution of the metric perturbations through the perturbed Einstein equations. The disadvantage of this approach is of course that it is hard to generalise to non-linear perturbations. The relative advantage of the latter approach is that we start out from the full nonlinear equations (20) and then linearise around a given background, in our case it would be the FLRW background. Gravitational wave perturbations in the metricperturbation approach are the tensor modes born out of perturbations to the ij components of the metric tensor and then traced through the Einstein equations. In contrast, in the covariant, gauge-invariant approach, gravitational wave perturbations are expressed as curvature perturbations that propagate and manifest themselves in the evolution of the electric and magnetic parts of the Weyl tensor. Pure tensor modes must be transverse and tracefree and therefore cause the divergence of the electric and magnetic parts of the Weyl tensor to disappear, i.e., D b E a b = 0 and D b H a b = 0, as well as the divergences of the shear anisotropy tensor and the anisotropic stress to disappear, D b σ a b = 0 and D b π a b = 0. In the linearised limit around FLRW and the presence of anisotropic stress of the form of (10), one can take a time derivative of (65a) and use (65b) to derive a wave equation for the shear anisotropies σ ij [100], reminiscent of the wave equation obeyed for gravitational waves, In fact, perturbing the spatial metric as where γ ij is the transverse and traceless tensor perturbation corresponding to the gravitational wave perturbation, the shear anisotropy tensor is related to the metric tensor perturbation as follows (e.g., [122]), This is because the shear tensor ultimately is the traceless part of the expansion tensor defined by (4), which is related to the time derivative of the metric variables in a homoegeneous spacetime. In drawing the equivalence between the metric perturbation approach to perturbation theory and the covariant gauge invariant approach, this relation would allow us to recover the familiar evolution equation for the metric tensor modes γ ij through Eq. (65a). 15 The corresponding equation is of the form agreeing with, e.g., [57,123]. In the infrared limit, i.e., on large super-Hubble scales where ∂ 2 /a 2 → 0, the equation becomes This can also be found by substituting (68) into the previously derived equation (25b), which makes the connection between anisotropies and long-wavelength gravitational waves explicit.
A key aspect of the above, either viewed through (66) or (69), is that shear and equivalently gravitational waves receive a damping factor (in the form of a friction term) 15 We can also see this by noting that the electric part of the Weyl tensor E ab is related to the traceless part of the 3-Ricci tensor denoted by (3) R ab as R ab = E ab + 1 2M 2 Pl π ab − Hσ ab + σ c a σ b c This relation is taken to be in the absence of vorticity, as in all of this work. The full equations are found in the Appendix of [100].
due to the presence of viscosity with η > 0. The negativity of the Hubble parameter in a contracting universe typically implies the growth of shear and of gravitational waves (most easily seen on super-Hubble scales). 16 The viscosity coefficient can counterbalance this effect though, such that anisotropies are damped (resulting in isotropisation) and so are gravitational waves. In fact, in the FLRW limit, one can solve (70) for the longwavelength ∂ t γ i j in the same way we solved for σ i j in (46), from which we can translate the results. For a constant viscosity coefficient and a pressureless EoS, one finds an exponential damping initially [in the form of (27), where we should think of ρ σ being replaced by A similar result was derived in [47] for a constant coefficient of viscosity. In the context of matter bounce cosmology, this damping would not realistically resolve the large tensor-to-scalar ratio problem if the viscosity is coming from a dilute gas of black holes (for the same reason it could not realistically lead to isotropisation within the regime of validity of the approximations). For an interacting field theory at finite temperature with η ∝ 1/a 3 and a radiation EoS, one recovers exponential damping in the form of (29). For a dense black hole gas with η = κ|H| (when H < 0) and a stiff EoS, one finds in a similar way to (31) that ρ GW ∝ 1/a 2(3−2κ/M 2 Pl ) , and hence gravitational waves are completely damped out by the time a → 0 if κ > 3M 2 Pl /2.

VI. DISCUSSION AND CONCLUSIONS
Bouncing cosmologies present an alternative to traditional expanding cosmologies by avoiding an initial singularity. The expense occurs by hypothesising some possible new physics at the bounce, which causes the universe to re-expand after an initial phase of contraction. However, there are a few problems regarding the growth of anisotropies and inhomogeneities in the contracting phase itself. Traditionally, a phase of ekpyrosis, where a fast-rolling scalar field mediates a slow contraction, exhibits an effective EoS p ρ and is able to dominate over the anisotropies and inhomogeneities.
Other dissipative mechanisms, such as particle creation and other quantum effects (e.g., [109,[126][127][128][129]), a nonlinear EoS (e.g., [63,103]), and the introduction of shear viscosity have been studied in the context of anisotropy reduction. In this work, we have studied possible microphysical realisations of such a dissipative model of shear viscosity. We have studied this in the context of a gas of black holes, both in the dilute and the dense limit. We find that the coefficient of viscosity remains constant and is temporarily effective in suppressing anisotropies in the dilute limit. However, the viscosity approximation is violated unless the viscosity coefficient is small enough, in which case isotropisation cannot occur. In the dense black hole gas case (which is considerably more speculative), we have the beginnings of a microphysical picture of understanding how a coefficient of viscosity that scales with energy density as η ∝ ρ 1/2 can be realised and, as has been seen in the literature, give rise to successful isotropisation and lead to a Friedmann singularity (if allowed to evolve to a crunch) even in the most general of anisotropic spatially homogeneous universes.
Another microphysical example that we have studied is the case of a λφ 4 interacting scalar field theory at finite temperature. The effective evolution of the background is that of a radiation-dominated universe. We studied the evolution of anisotropies in the case of a flat Bianchi type-I universe containing only expansion anisotropies, as well as in the case of a spatially curved closed anisotropic Bianchi type-IX universe. We found that in both cases the viscous damping dissipates the energy density in the anisotropy into radiation. The viscosity approximation itself remains valid in both cases, at least for enough efolds for the exponential suppression of anisotropies to be effective, under assumptions of high initial temperature and a universe that does not start out curvature dominated deep in the contracting phase for the case of Bianchi IX. Similar results have been found in the same context, but using different analyses and in the context of particle creation and semi-classical gravity [129]. Finally, as the anisotropy tensor itself is related to the time derivative of the tensor modes, the effect of the shear dissipation is equivalent to a damping of the amplitude of long-wavelength gravitational waves (see, e.g., [111,130] for additional implications of this principle).
While the λφ 4 model is an interesting toy model, which successfully manifests isotropisation, it does not constitute a complete theory of the very early universe. In particular, it cannot explain the formation of structures, i.e., it does not generate a nearly scale-invariant spectrum of curvature perturbations on large scales by itself. The addition of a spectator field (e.g.,à la curvaton [131]) could potentially resolve this issue, but this would require further investigation, especially with regard to the competition between quantum and thermal fluctuations in such a model. Alternatively, a contracting λφ 4 model could be part of a larger scenario that includes a period of inflation (e.g., [132][133][134]), which takes care of generating the right perturbations.
For the matter bounce scenario, where scale-invariant curvature perturbations are generated during a phase of matter-dominated contraction, it appears viscosity can serve as an isotropising mechanism to keep the model close enough to FLRW. However, this remains phenomenological since viscosity is actually hard to generate in a fluid that weakly interacts by definition. For example, we showed in this paper that a dilute gas of black holes could not realistically provide sufficient vis-cosity to keep the universe isotropic. Thus, unless one modifies the gravitational theory, e.g., with a graviton mass [32], which suppresses both anisotropies and gravitational waves, or with a specific non-minimal coupling to gravity (e.g., [135][136][137], but see also [138]), the matter bounce scenario remains unviable.
In any more realistic bouncing scenario hoping to explain the origin of the cosmic microwave background, one has to be aware that requiring isotropy with Ω σ < 1 for a certain number of e-folds might not be sufficient. Indeed, Ω σ might have to be several orders of magnitude below unity for the bounce to be achievable and for cosmological perturbations not to receive significant contributions from the shear. This is due to the fact that shear enters as a source term in the scalar, vector, and tensor perturbations of an anisotropic universe such as Bianchi I (see, e.g., [139]). Therefore, one expects an upper bound on the size that σ 2 may be allowed to reach in any given scenario [140].
Another aspect that needs to be taken into consideration in a more realistic scenario is the presence of shear due to quantum fluctuations in addition to the classical anisotropies discussed in this work. For instance, stochastic fluctuations of a scalar field could produce an anisotropic stress sourcing shear. However, when the background EoS satisfies w ≥ 0 as studied in this work, the resulting quantum shear only becomes dominant near the Planck scale [141]. Therefore, any 'lowenergy' bounce could evade this issue, though it remains an important contribution to shear that needs to be considered seriously in light of the previous paragraph.
Let us end by commenting on the dense black hole gas. As we mentioned, this remains the only known model resulting in η ∝ √ ρ and thus in full isotropisation within the approximations. Such a gas remains a fairly exotic toy model though. To start, the possible formation channels of such a gas remain hand-wavy; dealing with large inhomogeneities and their collapse into black holes would certainly have to be tackled numerically as in, e.g., [90,142]. Also, there is a great lack of understanding of the evolution of black holes embedded in cosmological backgrounds (apart from approximately Schwarzschild-de Sitter and McVittie spacetimes -see, e.g., [143][144][145][146][147]), and refining the corresponding approximations made on that front would definitely improve the description of the dense black hole gas. Nevertheless, if such a gas could really exist in nature in the very early universe (near a crunching singularity for instance), it remains interesting to ask the question of what could be the possible subsequent evolution of the gas. Could the black holes pass through a bounce and become primordial black holes as suggested in [89][90][91][92] or evaporate into remnants accounting for dark matter [93][94][95]? Could the black holes become stringy in nature at high energies and be part of a greater string-cosmology scenario [87,88]? Or could the black holes evaporate and emit specific electromagnetic signals or merge and emit specific gravitational-wave signals [148,149]? All those questions deserve closer scrutiny and could open up the path to a new understanding of the physics near the highest cosmological energy scales.