Viscous attenuation of gravitational waves propagating through an inhomogeneous background

We consider the propagation of gravitational waves in the late-time Universe in the presence of matter distribution inhomogeneities, and we also consider the cosmic fluid to be viscous. In this work, we investigate the cumulative effect of inhomogeneities and viscosity of the cosmic-fluid on the observables associated with the sources of the gravitational waves. Employing Buchert’s averaging procedure in the backreaction framework, we consider a model of spacetime in which matter is distributed in-homogeneously across space. Using the modified redshift versus distance relation, through the averaging process in the context of the model, we study the variation of the redshift-dependent part of the observed gravitational wave amplitude for different combinations of our model parameters while simultaneously considering damping of the gravitational wave amplitude due to viscosity of the cosmic-fluid. Then, we investigate the differences occurring in the variation of the redshift-dependent part of the observed gravitational wave amplitude due to consideration of viscous attenuation. We show that there are significant deviations after the inclusion of viscous attenuation in our analysis, depending on the chosen value of the coefficient of viscosity. Our result signifies the importance of the effect of viscosity, within the model of an inhomogeneous Universe, on precision measurements of parameters of compact-binary sources of gravitational waves.


Introduction
The century-old prediction of the possible existence of Gravitational waves (GWs) in Einstein's theory of General Relativity [1,2] has recently found confirmation from the Laser Interferometer Gravitational-Wave Observatory (LIGO), and Virgo scientific collaborations [3][4][5][6][7][8], which has opened a new window to decipher the mysteries of Universe.With more and more GW data pouring in, one expects that the GWs will provide more insight into diverse phenomena, such as the origins of black holes, the extreme conditions inside neutron stars, the chronicle of how the Universe structured itself into galaxies, the physics of the first few moments in the aftermath of the Big-bang and the standard picture of Universe itself.
Since gravitational wave observations are used to infer various fundamental properties related to the source of emission, it is important to have a complete understanding of the physics of their propagation from the source to us through the intervening background which in the present Universe is dominated by the dark components, viz., dark matter and dark energy.Properties of the cosmic fluid could be cause for significant attenuation of the amplitude of GWs propagating through it.It may be noted that though electromagnetic (EM) waves are not affected by the viscosity of matter, GWs are indeed affected by viscosity [9][10][11][12][13][14][15].GWs have to work against the viscous matter while passing through it, resulting in loss of its energy, which is manifested in attenuation of amplitude or damping of the GWs. 1 Studies [16] have been done analyzing the interaction between a viscous fluid (the primordial plasma in this case) and the primordial gravitational waves using a relativistic hydrodynamic theory.There also exist studies [17,18] describing damping of GWs due to non-collisional media during the propagation in alternate gravitational theories.For example, it is possible for the longitudinal scalar modes of the GWs from Horndeski theories, to be damped by the non-collisional ensemble of massive particles.
Cosmologists have used viscosity in wideranging studies over the years.The initial singularity at the big bang can be avoided by invoking shear, and bulk viscosity [19,20].Viscosity has been used to explain dark energy [21][22][23], and it has been shown that the Universe's accelerated expansion can be due to the effect of viscosity [19,[22][23][24][25][26][27].In certain other schemes, the neutrino mass [28] and the 21-cm emission temperature [29] can also be constrained using viscous cosmology.The viscous matter in the path of propagation of GWs could be most likely in the form of certain types of dark matter [24,[30][31][32][33], though even some of the visible stellar matter may also have some amount of viscosity.Such dark matter with dissipative dynamics induced by viscosity can possibly settle the strain between Planck Cosmic Microwave Background (CMB) and Large Scale Structure (LSS) observations [34].GWs in the presence of viscosity have been suggested as probes of such viscous cosmological models [14].
Cosmological observations have also revealed that though the Universe is smooth and uniform on the very largest of scales, given the standard picture of the Big Bang and the known age of the Universe, this is not true for smaller scales.The transition from homogeneity to inhomogeneity at smaller scales has been indicated by various tests from cosmological observations like the WiggleZ Dark Energy Survey [35].Analysis based on Sloan Digital Sky Survey [36] have revealed absence of homogeneity in the large-scale galaxy distribution.Current estimates to analyze large-scale fluctuations in the luminous red galaxy samples based on higher-order correlations have found significant (more than 3 σ) deviations from the Λ cold dark matter (ΛCDM) mock catalogues on samples as large as 500 h −1 Mpc [37].Thus, inhomogeneities due to structures may have important effects on length scales even as large as 500 h −1 Mpc.
It has been argued that local inhomogeneities may impact the overall evolution of the Universe, through the backreaction arising from the process of averaging [38][39][40][41][42][43][44].The backreaction effect quantifies the non-linear process of structure formation on the mean global evolution of the Universe.Several works have been done to explain the accelerated expansion of the late-time Universe through backreaction [45][46][47][48], though there exists a debate in the literature [49] whether inhomogeneities could account for the accelerated expansion.Recently, the Hubble tension [50,51] arising from a discrepancy in the inferred value of the Hubble parameter from local measurements compared to that from early Universe physics, has attracted a lot of attention.It may be noted that backreaction induced curvature may possibly explain the larger values of the Hubble parameter obtained locally [52].
Generally, GW analysis is done by deeming that GW propagates through a homogeneous and isotropic spacetime, described by an FLRWmetric.However, GW sources which are the subject of the present observations [3][4][5][6][7][8], lie well within the scale of 500 h −1 Mpc.The analysis of the effect of inhomogeneities on the propagation of GWs may be of significance for precision measurements in the emerging field of GW astronomy.
The motivation of the present analysis is to study the propagation of GWs through the background containing viscous dark matter in the presence of inhomogeneities.It has been shown earlier that the inclusion of the effect of local inhomogeneities leads to a non-trivial impact on the propagation of EM waves in averaged spacetime [53][54][55][56][57][58][59].The cosmological-distance versus redshift relations gets modified due to averaging over inhomogeneities [60][61][62], leading to interesting prospects for the detection of signatures of inhomogeneities through observations of Hubble expansion [63].GWs act as fellow carriers of information to EM waves, and analysis pertaining to the former opens up a new avenue to the physics of the evolution of the Universe starting from early times, offering insight into the nature of gravity itself, and ranging to the current cosmic acceleration.
GWs have particular relevance for those sources which do not emit any EM signals.Here we consider compact objects in binary formations, e.g., black hole -black hole (BH-BH) binaries, from which emitted GWs can be detected after traversing through the background viscous fluid.We consider the background dynamics arising from the backreaction of inhomogeneities due to structure formation.Specifically, we employ the Buchert formalism [43,44] to quantify the effect of backreaction.Buchert's approach on backreaction has been analysed earlier in various efforts to obtain concurrence with cosmological observations related to the current acceleration without resorting to dark energy [45-48, 53, 61, 64-67].It has been shown using the Buchert framework [68] that the amplitude of GWs produced from binaries could deviate substantially from that in the case of a homogeneous spacetime described by the ΛCDM model.
In the present work, we show that the local viscous-inhomogeneities in the path of propagation of GWs may have a considerable impact on the GW observables.The inclusion of viscosity affects the GW observables in ways different from the case of its absence [68].In the context of a simplified two-partitioned model of inhomogeneities within the context of the Buchert framework, here we evaluate the attenuation of GWs resulting from our model in comparison with the standard analysis of the ΛCDM model.Our approach clearly brings out the quantitative differences in the GW signal due to the inclusion of effects of viscosity and inhomogeneities, in combined as well as separate ways.
The paper is organized as follows.A brief description of the background dynamics is provided in Sec. 2.Here we first discuss a viscous ΛCDM model and next describe our inhomogeneous two-partitioned model (for both viscous and non-viscous cases) within Buchert's averaging formalism.In Sec. 3, the modification of the redshift-distance relation due to the averaging procedure is presented.In Sec. 4, the effect of local viscous inhomogeneities on the redshift-dependent part of GW amplitude is demonstrated.Finally, we present concluding remarks in Sec. 5.

Background dynamics
For ΛCDM model (without viscosity), the Hubble parameter is given by the Friedmann equation where a is the scale factor, H 0 is the present value for the Hubble parameter, Ω m0 denotes the fractional matter density components (assumed pressure less) of the Universe, Ω r0 is the fractional radiation density term, Ω k is the term related to the curvature, and Ω Λ denotes the cosmological constant component.In practice, the contribution of the radiation at late times (i.e., at the time of structure formation) is negligible compared to the matter and cosmological constant terms.Also, observations indicate that the geometry of the Universe is almost flat, viz., Ω k ≈ 0.

Viscous ΛCDM model (vΛCDM)
GW may propagate through dark matter in its path from its source to the observer.There are various theoretical models of dark matter.One of these is the Self Interacting Dark Matter (SIDM) model [30][31][32][33].In this model, self-interaction is introduced between the dark matter particles, which results in dissipation in the dark matter fluid.The outcome of this dissipation is the introduction of coefficients of shear, and bulk viscosities [24].In our approach, the dark matter behaves as a viscous/dissipative component.The general structure of this model is given by the field equation where R µν represents the Ricci tensor, R represents the Ricci scalar, g µν is the metric tensor, Λ is the cosmological constant and T µν stands for the energy-momentum tensor of the viscous matter.This tensor possesses both the perfect fluid structure as well as the possible dissipative effects such that [69,70] where ρ is the density, p is the pressure and the component ∆T µν is the viscous contribution to the fluid, Here ξ is the bulk viscosity, η is the shear viscosity, Θ = u µ ;µ is the expansion, u µ is the 4 -velocity and ";" represents the covariant derivative.For simplicity, we set the pressure, p = 0.Then, our dark matter possesses only the viscous pressure given by [71] p In the FLRW metric, the bulk viscous pressure reduces to Dark matter physics can incorporate some possible dissipative mechanisms [24,[30][31][32][33][34].Only the bulk viscosity remains compatible with the assumption of large-scale homogeneity and isotropy.The other processes, like shear and heat conduction, are directional mechanisms that decay as the Universe expands.Shear viscosity has mostly been neglected in these studies on the grounds of not contributing to a homogeneous and isotropic universe, which is undoubtedly true at the large-scale background level [71,72].Hence, for our purpose, for this viscous ΛCDM model, shear viscosity does not contribute to the dynamics of an isotropic and homogeneous background.However, shear viscosity does indeed play a role in the attenuation of gravitational waves, as we will see later in Sec. 4. Now, for a viscous ΛCDM model, using the FLRW metric, the Friedmann equation reads, Here, ρ v stands for the density of viscous matter and denotes all the matter components.We have assumed that all the matter components are endowed with viscous properties.For our purpose here, a proper separation between baryons and dark matter is unnecessary.One should also recall that baryons contribute about 1/6 th of the present total matter distribution; thus, this is not expected to lead to appreciable changes in our analysis, as baryonic matter is a subdominant component in comparison to dark matter.By defining the fractional densities , where H 0 is the present value for the Hubble parameter, the Friedmann equation (Eq.7) becomes Using now the fluid equation for ρ v , one gets Using (Eq.6), (Eq.9) can be recast as an equation for the fractional density Ω v as where we have defined the fluid equation of state parameter for the viscous dark matter fluid, ω v , as Using this formalism, Ω v as a function of the redshift z is calculated.The corresponding quantity in the ΛCDM case is Ω m0 (1 + z) 3 (Eq.1).

Buchert's formalism in a two-partitioned model
A popular approach for studying the effect of inhomogeneities is based on an averaging framework, and several averaging techniques have been proposed [38][39][40][41][42]60].Since the Einstein equations are non-linear, the solutions for an overall homogeneous matter distribution differ from the averaged solution for a general locally inhomogeneous matter distribution.In other words, the evolution of the homogeneous Universe at large scales may be slightly different from that of an inhomogeneous Universe, even though inhomogeneities, when averaged over a sufficiently large scale, might be negligible.The difference between the evolution of these models of the Universe gives the backreaction effect.It quantifies the non-linear effect of structure formation on the mean global evolution of the Universe.
In the averaging framework, [43], the problem is simplified and restricted to scalar quantities only.Einstein equations are decomposed into a set of dynamical equations for scalar quantities.Averages on flow-orthogonal spatial hypersurfaces (vorticity is assumed zero) are defined as proper volume averages.Under certain assumptions, this leads to Buchert's equations with a kinematical backreaction term [43].Such an approach has sparked considerable interest as it has been shown that backreaction could lead to an agreement with cosmological observations without resorting to dark energy [45-48, 53, 61, 64-67].For details of Buchert's averaging procedure, one may refer to Ref. [43,44].Here we provide a brief overview of Buchert's formalism required in context of the present analysis.In Buchert's averaging scheme for scalars, averages of scalar quantities on flow-orthogonal spatial hypersurfaces are defined as where D is a spatial domain.One fundamental quantity characterizing this domain is its volume, which is given by, The normalized dimensionless effective volume scale factor a D is defined by which is normalized by the volume V D0 of the domain D at some reference time t 0 which we can take as the present time.
Spatially averaging the Raychaudhuri equation, the Hamiltonian constraint and the continuity equation, one obtains the equations for effective scale factor in Buchert's formalism, which are respectively, where local averaged matter density ρ D , averaged spatial Ricci scalar R D and the Hubble parameter H D are domain dependent and are functions of time.Λ is the cosmological constant (Λ = 0 for our model and for our purpose here).Q D is called the backreaction term which quantifies the averaged effect of the inhomogeneities in the domain D and is defined as where θ is the local expansion rate and σ 2 := 1 2 σ i j σ j i is the shear-scalar.Q D is zero for a homogeneous domain.The departure from homogeneity is ingrained in this term.Q D and R D are inter-related by the equation: (Eq. 19) couples the time evolution of the backreaction term with the time evolution of averaged intrinsic curvature and signifies the departure from FLRW-cosmology, where there is no such coupling.
In this framework, the domain D is partitioned into non-interacting subregions F l composed of elementary space entities ∅ for all α = β and l = m.The average of any scalar function f on the domain D is given by, where λ l = V F l /V D is the volume fraction of the subregion F l such that l λ l = 1 and f F l is the average of f on the subregion F l .The above equation governs the averages of scalar quantities ρ, R and H D .But Q D , due to the presence of θ 2 D term, does not follow the above equation.Instead, the equation for where Q l and H l are defined in the subregion F l in the same way as Q D and H D are defined in the domain D [44].
We can also define scale factor a l for the individual subregions in the same way as a D has been prescribed for the domain D. Since, the domain D comprises the different subregions F l and all these subregions are disjoint, therefore V D = l V F l , which results in a 3 D = l a 3 l .Twice differentiating this relation with respect to foliation time gives us, Within the context of the above framework, we consider a two-partitioned model.In our model, the domain D, through which GWs propagate from the source to the observer, is an ensemble of two types of disjoint FLRW regions [73].These are -(i) the overdense region described by an FLRW region which can be assumed to be spatially flat, and (ii) the underdense region, which is a nearly empty FLRW region which is assumed to have a small density (compared to the overdense region) and negative spatial curvature.We consider the overdense region to have viscous nature in our present work to study the effect of viscosity on the GW amplitude.Therefore for our two-partitioned model, (Eq.21) effectively becomes, where λ o denotes the volume fraction of the overdense region.Now, (Eq.19) couples the time evolution of the backreaction term Q D with the time evolution of the averaged 3-Ricci scalar curvature.Therefore, we can choose the curvatures of our individual sub-regions in such a way that the Q l term for these sub-regions becomes effectively zero [44], i.e., Q o = 0 and Q u = 0.This has been done by taking the curvature of our overdense region to be zero, i.e., our overdense region is flat.On the other hand, we have assumed our underdense region to have Friedmann-like a −2 u constant curvature term.These assumptions along with (Eq.19) results in, Q o = 0 and Q u = 0.This stipulation to FLRW is an approximate assumption governing our model (in the more general case, the sub-domains may not necessarily be FLRW regions).From (Eq. 23), it can be seen that control over global backreaction can be achieved only if the individual backreaction terms are not set to zero.
In this work, we investigate the simultaneous impact of inhomogeneities and viscosity of the matter contained in the overdense region, on the amplitude of GWs for the same type of sources.For this, we consider two cases -(a) when the GW is assumed to propagate through a homogeneous and isotropic spacetime (with and without viscosity of associated matter), described by the FLRW metric, in the ΛCDM model, and (b) when the GW propagates through an inhomogeneous spacetime (both viscous and non-viscous cases), described by our model specified above.The physical interpretation of the underdense region in our model is that they represent the cosmic voids in the path of propagation of the GWs, and the overdense region represents all the matter content in that path.In this work, we have considered both scenarios (viscous and non-viscous) for the matter present in the overdense region of our model.A point to note here is that we are not considering viscous matter in addition to non-viscous matter here.The total matter content in both the non-viscous case and the viscous case is the same.It's just that in viscous case the matter content has viscous properties too.The overdense viscous region in our model portrays a viscous dark matter fluid.

Non-viscous case (no bulk and shear viscosity)
Expressions for the scale factors for the two types of regions, overdense (non-viscous in this case) and underdense, for this model, are given as [65,66], Here o represents the overdense region, and u represents the underdense region.In this model, the scale factors of the two regions are proportional to cosmic time raised to some powers α and β for over and under-dense regions, respectively.α varies from 1/2 to 2/3 since the evolution of a u is expected to be faster than that for the radiationdominated case (1/2) and limited by the maximum value for the matter-dominated case (2/3).β varies from 2/3 to 1 to denote any behaviour ranging from a matter-dominated region (β = 2/3) up to dark energy-dominated region (β > 1).c u and c o are constants of proportionality and are given respectively as: where a D0 is the scale factor at the present time for the domain D. For present time t 0 ≈ 13.8 Gy, a D0 = 1 and H 0 = 100 h km sec −1 M pc −1 .Therefore, Using (Eq.22) for our model we get, Here λ o = V o /V D is the volume fraction for the overdense region, V o is the volume of the overdense region and λ u is the volume fraction of underdense region such that λ u + λ o = 1.Now, using (Eq.14), V o can be written in terms of scale factor and initial volume fraction, , where V o0 is the volume of the overdense region at some reference time t 0 .This in turn gives us, where is a constant and λ o0 , a D0 , t 0 are the present volume fraction of overdense region, the present global scale factor and the present time, respectively.(Eq.29) shows that the volume fraction of the overdense region is a function of α and β (through a D as a D is a function of both α and β).Similarly, it can be shown that the volume fraction of the underdense region is a function of α and β.This implies that (α, β) governs the volume fractions of the 2 subregions in our 2-domain model.The present values of volume fractions of the 2 region are taken as (λ o0 , λ u0 ) = (0.09, 0.91) [44].Solving (Eq.28) gives us an expression for a D , and using that expression in (Eq.15) gives us an expression of ρ D for our model.

Viscous case
The sub-regions in our backreaction model are essentially FLRW regions.Hence, their dynamics are also governed by the standard Friedmann equations (Eq. 1, 7).Only our overdense region has viscous matter contained in it.So, for our overdense region, the Friedmann equation reads (we have taken Λ = 0), Here, a v o is the scale factor for our viscous overdense region, H v o is the Hubble parameter for our viscous overdense region and ρ v o stands for density of viscous matter in our overdense region and denotes all the matter components.
We have assumed that all the matter components of our overdense region are endowed with viscous properties similar to the ΛCDM model with viscosity.Also, since our overdense region is an FLRW region (homogeneously and isotropically overdense), therefore here too the shear viscosity with coefficient η, doesn't have any effect on the dynamics of the background.The fluid equation for ρ v o is given by, where p v o is the bulk viscous pressure given by By defining the fractional densities Ω v o = 8πGρ v o /(3H 2 0 ), equation for the fractional density Ω v o is given as, where we have defined the fluid equation of state parameter for the viscous dark matter fluid for our overdense region, ω v o , as Since (Eq.33) is a first-order differential equation, we require one boundary condition to solve it.We take this boundary condition from our non-viscous backreaction model, given by where Ω o is the fractional density of the overdense region from the non viscous case.The scale factor for our overdense region in our viscous backreaction model is calculated using (Eq.34, 33 and 30).The boundary condition used for solving (Eq.33) is (Eq.35).The a v o that we get from the above-mentioned analysis is a function of (α, β).The scale factor for the underdense region, a u is still given by (Eq.25).The scale factor for domain D in this case a v D is given by, The equation for the backreaction term in this case Solving (Eq.36) gives us an expression for a v D .(Eq. 37) leads to an expression for Q v D .Using these expressions in (Eq.15) we get an expression of ρ v D for our viscous backreaction model.The present value of the bulk viscosity parameter has been estimated in the literature based on various theoretical considerations and observations [15,27,74,75].By solving the energy conservation equation for bulk viscous flat Friedmann universes, the coefficient of bulk viscosity at present time ξ has been estimated to be ∼ 10 6 Pa sec in [75].Other analyses based on comparing the theoretical curves for H = H(z) with observations [76], and by studying the asymptotic behaviour in the equivalent phase space in a Friedmann model of the Universe with bulk viscous matter [77] have also been performed.In [78], a value of ∼ 10 7 Pa sec for ξ has been suggested.Considering all the above studies, a reasonable range for ξ may be taken as 10 5 Pa sec < ξ 0 < 10 7 Pa sec.For our present work, we have used the mid-range value of 10 6 Pa sec in the red-shift range 0 ≤ z ≤ 5.
In Fig. 1 we provide a plot of the Hubble parameter versus the redshift for ΛCDM, viscous ΛCDM and for our viscous backreaction model for various values of the backreaction model parameters.The underdense region in our model represents the cosmic voids in the path of propagation of the GWs, and the overdense region represents all the matter content in that path where the effect of viscosity is directly evident.In Fig. 1, for the plots of our viscous backreaction model, we fix the value of β at 0.8 (representing a mid-range value in the range of values of β -(2/3, 1)), and we vary the values of α.As expected, for a fixed β, curves with larger value of α lead to larger values of Hubble parameter.It can be seen that backreaction from inhomogeneities may lead to departure in the background Hubble evolution that may get accentuated for higher redshifts.

Redshift and distance relation
Buchert's averaging scheme provides us with a method of spatially averaging scalar quantities in the backreaction framework.Such quantities need to be related to cosmological observables.A possible approach relies on the study of distance-redshift relation in an inhomogeneous universe using an approximate metric [60].In our present investigation, we consider a more consistent scheme based on the procedure of averaging.The covariant scheme proposed by S. Räsänen [54,79] provides us with a way of doing this.This scheme gives the relation between effective redshift and angular diameter distance D A in the following way, (Eq. 38) provides an expression of effective redshift z in terms of the scale factor a D of the domain D. Certain conditions must be satisfied to apply the covariant scheme.These are: (i) spatial averages are determined on hypersurfaces of statistical homogeneity and isotropy, and (ii) structure evolution is slow compared to the travel time of GW from source to the observer.For our model, domain D is the region of spacetime through which the GW travels while propagating from the source to the observer.Domain D could have any combination of fractions of underdense and overdense regions, i.e. any combination of (λ u , λ o ) as long as λ u + λ o = 1 is satisfied.(λ u , λ o ) in turn are governed by (α, β).In this work, we take various combinations of allowed values of (α, β) in our analysis.Using the expressions for H D and ρ D that we calculated from our model for the two cases -viscous (with only ξ, but no η since η doesn't affect the background) and non-viscous, and using the covariant scheme (Eq.38, Eq. 39), we can calculate D A for our model.
In (Fig. 2), we plot the ratio of angular diameter distance for our model to the present Hubble-length (D H = cH −1 0 , value of H 0 (present value of H) used is 100 h km s −1 M pc −1 ≃ 0.07 Gyr −1 , with h = 0.7) as a function of effective redshift with different combinations of (α, β).We have studied earlier the effect of varying α on the dynamics (Fig. 1).Here in Fig. 2 we explore the effect of varying β on the angular diameter distance choosing a fixed value of α = 0.5.The value of ξ used is 10 6 Pa sec in the range 0 ≤ z ≤ 5 [15].From this figure, one can see that for low redshifts (z < 0.5), curves for our model overlap with each other and with the ΛCDM curve, but as we increase the redshift, curves start deviating from the ΛCDM curve.Inclusion of viscosity in the analysis results in deviation in the plots for both the ΛCDM model and our model.For the ΛCDM case, plot of viscous case (dashed) has higher magnitude than the non-viscous case (solid).It is observed that for lower values of z (z 2), plots of viscous cases (dashed curves) for our backreaction model have a higher magnitude of D A for the same value of z in comparison to corresponding non-viscous cases (solid curves) and for z 2, the magnitude of D A for viscous cases (dashed curves) are smaller than those for the non-viscous cases (solid curves) for our model.
In (Fig. 3), the ratio of luminosity distance D L to the present Hubble length D H is plotted with respect to redshift z for both viscous and non-viscous cases of ΛCDM and our backreaction model.D L is calculated from D A (Fig. 2) using the relation D L = (1+z) 2 D A .There are two insets in the figure.It can be seen here too that for low redshifts (z < 0.5), the curves for our model overlap with each other and with the ΛCDM curve, but as we increase the redshift, the curves start deviating from the ΛCDM curve.Inset (a) shows the magnified portion of the plots at lower values of redshift, particularly around the value of z for which the plots of D A /D H turn around in (Fig. 2).As in (Fig. 2), plots for the viscous backreaction model (dashed curves) have a larger magnitude than the non-viscous case (solid curves) for about z 2. For z 2.5, plots for the viscous case have a smaller magnitude than the non-viscous case similar to plots in (Fig. 2).This behaviour is shown in inset (b).A point to note is that this behaviour is not observed for the ΛCDM case, where viscous case plots have a larger magnitude throughout the range of z of our interest.In inset (b), purple (backreaction (β = 0.8)) and red (backreaction (β = 0.7)) plot lines overlap, which is similar to the case of (Fig. 2) at this value of z.Another feature of the plots for our backreaction model can be observed at higher redshift z.It can be seen that the difference between the solid and dashed lines (representing the non-viscous and viscous cases, respectively) increases with larger values of β.Specifically, the difference between the red lines is smallest for the lowest values of β, while the difference between the brown lines is largest for the highest values of β.

Gravitational wave amplitude
The amplitude of GW from a binary of compact objects of masses m 1 and m 2 in the early inspiral stage, where Keplerian approximations are well valid, is given by (for the cross(×)-polarization) [80] where ω is the observed angular frequency of the binary of compact objects and D L is the luminosity distance of the binary from the observer.For the plus (+)-polarization, the peak of the amplitude remains identical.For a constant observed frequency, the redshift-dependent part in the GW amplitude is (1 + z) 5/3 /D L .
The amplitude of GW given in (Eq.40) is derived without considering the effect of viscosity on the propagation of the GW.In Ref. [14], the authors have studied the effect of the viscosity of the cosmic fluid, particularly dark matter, on the GW amplitude and have estimated its effect on the GW amplitude.Assessing the impact of viscosity, after travelling a proper distance L = ar, the GW gets attenuated by the factor (see, Appendix A) [14] A where L * is the proper source distance for zero shear viscosity and γ = 16πGη, where η is the coefficient of shear viscosity for the region through which GW is propagating.The attenuation factor in terms of luminosity distance D L is given by : where D L = (1 + z) 2 L. D L * is the luminosity distance of the source, for zero shear viscosity.Therefore, the total redshift-dependent part of the attenuated GW amplitude (let's represent this quantity by F (z)) is given by For non viscous case (η = 0), The luminosity distance D L is calculated from angular diameter distance, D A as D L = (1 + z) 2 D A .Since, D L = (1 + z) 2 L, for the flat FLRW spacetime, D A gives a good measure of the proper distance L. D A for flat spacetime is given by where D C is the comoving distance which is given as where D H is the Hubble distance (= c H0 , where H 0 = 100 h km s −1 M pc −1 ≃ 0.07 Gyr −1 , with h = 0.7.) and E(z) ΛCDM ≡ Ω M (1 + z) 3 + Ω Λ for the ΛCDM model (nonviscous case) (throughout this work, parameters used for ΛCDM are Ω m = 0.31, Ω Λ = 0.69. ) and E(z) vΛCDM ≡ √ Ω v + Ω Λ for the viscous ΛCDM model, where Ω v is the fractional density for viscous matter which has a contribution from the bulk viscosity ξ (subsection 2.1).Hence, now, for viscous ΛCDM case (vΛCDM ), where D Lv is the luminosity distance in the presence of viscosity, i.e. it is calculated using E(z) vΛCDM .For non viscous ΛCDM case, where D L * is the source luminosity distance in the absence of viscosity, i.e. it is calculated using E(z) ΛCDM .
The value of the shear viscosity parameter η has been estimated in earlier works to lie within the range ∼ 10 6 Pa sec and ∼ 10 9 Pa sec [14,15,81].Gravitational wave observation data from LIGO has been used to put constraints on the value of η [14].Ref. [81] used the relativistic Boltzmann equation to examine the theory of viscosities at the time of neutrino decoupling.The analysis of [81] to calculate the present-day value of η was subsequently employed in [15], to relate the value of η at the time of neutrino decoupling with the present value using a scaling relation, leading to the value of 10 8 Pa sec.Our present analysis uses the above value of 10 8 Pa sec as the most favored one.
In (Fig. 4), we plot the overall redshiftdependent part, F (z) (both F (z) vΛCDM & F (z) ΛCDM ), given in (Eq.46 & Eq.47) w.r.t.redshift z, for the range of values of η between ∼ 10 6 Pa sec and ∼ 10 9 Pa sec.The plots for viscous cases in (Fig. 4) have contributions from both bulk viscosity ξ (through D Lv ) and shear viscosity η.It was argued in Ref. [14] that since ξ only couples to scalar perturbations, it doesn't play a role in the attenuation of GWs, and only η affects the GW amplitude.However, from our analysis, it is clear that both ξ and η affect the GW amplitude, with the role of the former entering through the modified background dynamics due to bulk viscosity.The expressions given in (Eq.46 & Eq.47) are valid for a space with homogeneous mass distribution.To examine the variation of the redshift-dependent part of GW amplitude for our model, we have to modify the expressions accordingly.For our backreaction model, the relation D L = (1 + z) 2 D A is valid too, but in this case D A is calculated using the covariant scheme (Eq.38, Eq. 39).Our model represents a space with inhomogeneous mass distribution, where there are two types of regions, viz.over-dense and under-dense, and viscosity is only associated with the overdense region, as the under-dense region is assumed to be empty.
Hence, in the exponential term of the F(z) in (Eq.43), e − γ 2(1+z) 2 DL , for our model, D Lo would replace D L , where D Lo is the luminosity distance traversed by the GW through the overdense region only.So, the exponential term of the F(z) in (Eq.43) for our model is now given by e − γ 2(1+z) 2 DL o .In the ΛCDM model, incorporating viscosity results in an attenuation factor with D L in the exponential factor, where D L is the total luminosity distance traversed by the GW.Viscosity is not distributed through the entire path of the GW but is concentrated only in some regions.Using D L in the ΛCDM model results in a larger deviation between the attenuated and unattenuated cases in the ΛCDM model, compared to our model.
An important consideration for the propagation of EM waves for a model like ours (inhomogeneous 2-domain model) is that the ratio of distances travelled by EM waves through the two regions is equal to the ratio of proper volumes of the two regions [61].It is clear that this also holds for GWs.This condition gives, where D Lo and D Lu are the luminosity distances traversed by the GW through the over-dense and under-dense regions, respectively.If the total luminosity distance travelled by the GW in our model is D L v2d (v2d stands for viscous 2-domain inhomogeneous model, this D L v2d is calculated for the viscous case of our model (subsubsection 2.2.2) using D A calculated from the covariant scheme (Eq.38, Eq. 39)), then D L v2d = D Lo +D Lu , which gives, or, Therefore, the exponential term in the total redshift dependent part of attenuated GW amplitude for our model now becomes, and the total redshift-dependent part of the attenuated GW amplitude for our viscous model is given by, The redshift-dependent part of the GW amplitude for our non-viscous model (subsubsection 2.2.1) is given by, where 2d stands for non-viscous 2-domain inhomogeneous model and D L * 2d is the source luminosity distance for the case of our non-viscous 2-domain inhomogeneous model.Therefore, for our model, the deviation of the redshift-dependent part of GW amplitude due to viscous-attenuation, can be written as, To summarize, the redshift dependent part of GW for a homogeneous and non-viscous spacetime (ΛCDM model) is given by (Eq.47).Now, if we introduce matter distribution inhomogeneities, then F(z) gets modified due to modification of the redshift-distance relations, and is now given by (Eq.52).It can be seen that, Therefore, the effect of inclusion of matter distribution inhomogeneities on redshift dependent part of GW amplitude is equivalent to multiplication by the factor . On further introduction of viscosity in the analysis, GW in our backreaction model gets attenuated by the attenuation factor (using (Eq.50)), where D L * 2d is the luminosity distance for the case of our non-viscous 2-domain inhomogeneous model, calculated using the covariant scheme.Thus, the total redshift dependent part of GW in the presence of viscous inhomogeneities becomes (using (Eq.51),(Eq.52), (Eq.54) & (Eq.55)), In (Fig. 5), we plot the redshift dependent part of the GW amplitude, F (z) vs z for the ΛCDM model and for our model for model parameter (α, β) = (0.67, 1).There are 3 curves for each model.The solid curve represents the non-viscous case; the dotted curve represents the case in which only the bulk viscosity has been included in the analysis, and the dashed curve is the case with both bulk viscosity and shear viscosity in the analysis.As can be seen from (Fig. 5), even if we just consider bulk viscosity then also there is deviation of the redshift dependent part of GW amplitude (dotted curve) with respect to the non viscous case (solid curve).This is because redshift dependent part of GW amplitude consists of D L and ξ affects this D L via the quantity E(z) (Eq.45).On further inclusion of η in the analysis, redshift dependent part of GW amplitude gets attenuated (dashed curve).In this case, F (z) for viscous ΛCDM model is given by (Eq.46), for non viscous ΛCDM model by (Eq.47), for our non viscous model with inhomogeneities by (Eq.52) and for our viscous model with inhomogeneities, it is given by (Eq.51).In (Fig. 6), we plot the redshift dependent part of the GW amplitude, F (z) vs z for the ΛCDM model and for our model for different combinations of model parameter (α, β).Here, we have kept β constant = 0.8, and we vary the value of α.The solid curves represent the non-viscous cases, and the dashed curves represent the viscous case with contributions from both ξ and η.The redshift dependent part of GW amplitude, i.e., the quantity (1 + z) 5/3 /D L in the case of ΛCDM, has a minimum at z min ≃ 2.63 [82].In (Fig. 6) we display the minima points of the various curves.Circular points represent the minima for the solid (non-viscous) curves while triangular points represent the minima positions for the dashed (viscous) curves.We observe that the minima for the vΛCDM curve deviates significantly from the ΛCDM case.For our backreaction model (both viscous and non viscous), minima points deviate significantly from the ΛCDM case.It can be seen from the figure that for the non-viscous cases of our model, all three plot lines and their minima points for different values of α overlap.F (z) for the non-viscous cases is given by (Eq.52).Since the value of (1 + z) 5/3 is the same for all (α, β), the only quantity varying with changing (α, β) is D L .The volume fraction of the overdense region at the present time is taken as 0.09, and the range of variation of α, which governs the evolution of the overdense region, is from 0.5 − 0.67.Since the overdense volume fraction is so small; therefore variation of α over this small range (0.5 -0.67) doesn't have much effect, and hence, the plot lines for the non-viscous cases overlap in (Fig. 6).
From (Fig. 6), it can be seen that there is substantial attenuation of the redshift-dependent part of the GW amplitude due to viscosity.F (z) for viscous cases is given by (Eq.51).As compared to the non-viscous case (Eq.52), there are now additional terms.From (Eq. 24) and (Eq.26), the scale factor of the overdense region is given by, where t 0 = 13.8Gyr.For a given value of t (t < t 0 ), as we increase α, a o decreases.It can be seen in (Fig. 6), as we increase the value of α keeping β constant, the redshift dependent part of GW, F (z), for the viscous case, gets smaller in magnitude.Physically, this can be explained from the fact that from (Eq. 14) for an overdense region, a o can be defined as where V o0 is the volume of the overdense region at the present time which we have fixed for all values of (α, β) to be 0.09 by fraction of the total volume.As a o decreases, the volume of the overdense domain decreases.Hence, the distance travelled by the GW through the overdense region will be less, leading to less attenuation suffered by the GW.GW which travels the largest distance through the overdense region will suffer the most attenuation.This distance travelled by the GW is directly proportional to the scale factor of the overdense region.In (Fig. 6), the green dashed curve has the largest value of α (hence the smallest value of a o ), corresponding to the largest amplitude among the viscous cases since it has suffered the least attenuation.In contrast, the black dashed curve has the smallest value of α (hence the largest value of a o ), corresponding to the smallest amplitude among the viscous cases since it has suffered the most attenuation.The above feature is not observed for the non-viscous cases (solid curves) as the overdense region is non-viscous; hence, there is no viscous-attenuation.

Conclusions
In this work, we have studied the propagation of GWs from compact binary sources through a viscous inhomogeneous Universe governed by a model based on the averaging procedure for scalars in Buchert's backreaction framework [43,44].Dynamics under this model lead to a modification of the redshift-versus-distance relation from that for the ΛCDM model.The extent of variation depends on the combination of the model parameters.In the present work, we have considered viscosity to be present in the matter content within the over-dense regions of inhomogeneous spacetime described by our model, which causes the attenuation of GW amplitude when the GW passes through those regions of spacetime.We have incorporated the viscous attenuation of GW amplitude within our model of inhomogeneous spacetime and have derived an expression for the resultant redshift-dependent part of the GW amplitude.
In the ΛCDM model, incorporating viscosity results in a GW attenuation factor with D L in the exponent, where D L is the total luminosity distance traversed by the GW.It is worth noting that in the real Universe, viscosity resulting from dark matter interactions is not distributed uniformly through the entire path of the GW but is concentrated only in some regions.Therefore, in our model, we have used D Lo (luminosity distance of the overdense region) instead of D L .Using D L results in a more significant deviation between the attenuated and unattenuated cases for the ΛCDM model, compared to using D Lo for our model.It has been argued earlier [14] that since bulk viscosity couples only to scalar perturbations, it doesn't play a role in the attenuation of GWs.However, as shown here, bulk viscosity indirectly impacts the GW amplitude through its effect on the luminosity distance.Moreover, the effect of shear viscosity on GW attenuation is clearly demonstrated in the backreaction model due to inhomogeneities.
Our analysis demonstrates a substantial deviation in the redshift-dependent part of the GW amplitude due to the inclusion of viscous attenuation, compared to the case when viscosity is considered negligible or absent within the model of inhomogeneous spacetime.We have further shown that the rate of expansion of the overdense region (characterized by the parameter α, which governs the time evolution of the scale factor of the overdense region) plays a vital role in the magnitude of attenuation.It may be emphasized that consideration of the effect of viscosity on GW observables for compact binary sources is significant in the context of local inhomogeneities in the Universe.
To summarize, we would again like to highlight the importance of considering realistic background effects in the study of GW propagation, since using incorrect background dynamics for analyzing GW data from detectors could result in incorrect inferences about GW sources.Towards this end we have taken into account these two aspects of inhomogeneities and attenuation of GW due to viscosity together in our present work.Our analysis leads to several interesting features in the red-shift dependent part of the gravitational wave amplitude.Additionally, the role of bulk viscosity is a new feature that is also brought out in the form on its indirect contribution towards GW attenuation through modification of the background dynamics.Our analysis paves the way for obtaining more precise estimations of GW observables and bounds on the viscosity parameters of dark matter in future work involving more realistic backreaction models and data-analysis techniques in GW astronomy.
8πGT ij gives us the wave equation for GWs in a viscous fluid, ḧij + (3H + 16πGη) ḣij − ∇ 2 a 2 h ij = 0 (A5) ξ doesn't come into (Eq.A5) as it only couples to scalar perturbations.Next, performing the Fourier transform of (Eq.A5) and defining h ij as µ ij /a one gets, μij + (H + 16πGη) μij Defining conformal time τ as dt = adτ and using it in (Eq.A6) leads to where ′ denotes derivatives with respect to τ .On sub-horizon scales k 2 >> a ′ a , (Eq.A7) reduces to, The strain h ij of the GW in cosmic time t now becomes, where L = ar is the source distance, L 0 is the source distance for zero shear viscosity and ω p = ω a is the physical angular frequency.Therefore, the attenuation factor is given by L0e − β 2 L

Fig. 1 6
Fig. 1 Plot of Hubble parameter vs redshift z for the ΛCDM model, viscous ΛCDM model and our viscous backreaction model.Plots for our viscous backreaction model are for varying values of α with fixed value of β = 0.8.The value of ξ used is 10 6 Pa sec.The value of H 0 used is 100 h km s −1 M pc −1 ≃ 0.07 Gyr −1 , with h = 0.7.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.The inset shows the magnified portion of the plots at higher redshift.

Fig. 2
Fig. 2 Plot of the ratio of angular diameter distance D A to the present Hubble length D H w.r.t.redshift, for the ΛCDM case and for our backreaction model with different combinations of β and fixed value of α = 0.5 for both viscous and nonviscous cases.Value of ξ used is 10 6 Pa sec.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.The inset shows a magnified portion of the plots for our backreaction model.

Fig. 3
Fig. 3 Plot of the ratio of luminosity distance D L to the present Hubble length D H w.r.t.redshift, for the ΛCDM case and for our backreaction model with different combinations of β and fixed value of α = 0.5 for both viscous and non-viscous cases.Value of ξ used is 10 6 Pa sec.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.There are two insets in the figure.Inset (a) shows the magnified portion of the plots at lower values of redshift, from z = 1.4 to z = 1.6.Inset (b) shows the magnified portion of the plots at higher values of redshift, from z = 4.5 to z = 5.

Fig. 4
Fig. 4 Plot of F(z) for the ΛCDM model for both viscous and non-viscous cases.Value of ξ used is 10 6 Pa sec.Plots for viscous cases are for the applicable range of values of η between ∼ 10 6 Pa sec and ∼ 10 9 Pa sec.Dashed orange plot represents the favored value of η.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.The inset shows the magnified portion of the plot for the ΛCDM model and for our backreaction model (η = 10 6 Pa sec).

Fig. 5
Fig. 5 Plot of F (z) vs z for ΛCDM model and for our model for (α, β) = (0.67,1).The solid curves are for non-viscous cases.The dotted curves are for those cases in which only ξ has been included.The dashed curves have contributions from both ξ and η.Value of ξ used is 10 6 Pa sec.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.Inset shows the magnified portion of the solid curve and the dotted curve for the two models to illustrate the difference between these two curves.

10 8Fig. 6
Fig. 6 Plot of F (z) vs z for ΛCDM model and for our model for different combinations of (α, β) while keeping β = 0.8 as constant.The solid curves are for non-viscous cases.The viscous cases represented by dashed curves have contributions from both ξ and η.Value of ξ used is 10 6 Pa sec.Parameters used for ΛCDM are Ωm = 0.31, Ω Λ = 0.69.Circular points represent the minima for the solid (non-viscous)curves while triangular points represent the minima positions for the dashed (viscous) curves.Inset shows the magnified portion of the three dashed curves of our backreaction model to illustrate the difference between these curves.