Evolution of linear perturbations in spherically symmetric dust spacetimes

Evolution of linear perturbations in spherically symmetric dust spacetimes

S February, J Larena, C Clarkson and D Pollney Astrophysics, Cosmology and Gravity Centre, and, Department of Mathematics and Applied Mathematics, University of Cape Town, Rondebosch 7701, South Africa. Department of Mathematics, Rhodes University, Grahamstown, 6140 South Africa sean.february@uct.ac.za

We present results from a numerical code implementing a new method to solve the master equations describing the evolution of linear perturbations in a spherically symmetric but inhomogeneous background. This method can be used to simulate several configurations of physical interest, such as relativistic corrections to structure formation, the lensing of gravitational waves and the evolution of perturbations in a cosmological void model. This paper focuses on the latter problem, i.e. structure formation in a Hubble scale void in the linear regime. This is considerably more complicated than linear perturbations of a homogeneous and isotropic background because the inhomogeneous background leads to coupling between density perturbations and rotational modes of the spacetime geometry, as well as gravitational waves. Previous analyses of this problem ignored this coupling in the hope that the approximation does not affect the overall dynamics of structure formation in such models. We show that for a giga-parsec void, the evolution of the density contrast is well approximated by the previously studied decoupled evolution only for very large-scale modes. However, the evolution of the gravitational potentials within the void is inaccurate at more than the 10% level, and is even worse on small scales.


1 Introduction

We present results from a numerical code implementing a new method that solves the first order gauge-invariant linear perturbation equations in a Lemaî tre-Tolman-Bondi (LTB) background. LTB models are spherically symmetric but inhomogeneous dust solutions of the Einstein field equations. Compared to a Friedman-Lemaître-Robertson-Walker (FLRW) background, perturbations around a LTB background are complicated by the fact that they cannot be decomposed into the standard scalar vector and tensor modes, as the reduction of symmetry (through radial inhomogeneity) causes these modes to couple [1].

The method presented here can be used to model a variety of different astrophysical/cosmological scenarios. For example:

Relativistic corrections for structure formation.

Currently structure formation in cosmology is modelled either using non-linear models of Newtonian systems, or relativistically but only in the linear regime. This leaves an important area unexplored: non-linear relativistic aspects of structure formation – see, e.g., [2] for a review. Certain aspects of this have started to be taken into account in N-body methods which make use of relativistic corrections to the potentials [3, 4]. Our model can be used to analyse the growth of structure on top of a strongly non-linear background — either an over-density such as a cluster, or a large void, both of which generate large curvature and shear. The coupling of density perturbations to vector and tensor degrees of freedom can be explored and the errors induced by neglecting this coupling quantified.

Evolution of perturbations in void models.

If we were to live in a large, underdense void of a few giga-parsecs in diameter, distant supernovae would appear fainter than expected in a FLRW model, and the dark energy phenomenon could be explained on purely relativistic terms without invoking any new physics (see, e.g., [5, 6, 7, 8] and [9] for a comprehensive review). Clearly, structure formation places a key constraint on such models by probing their difference from the concordance model, and so serves as a test of the Copernican principle [10]. Structure formation in LTB models has only been quantified for the special case which neglects the coupling of the scalar gravitational potential to vector and tensor degrees of freedom [11, 10, 12, 13, 14, 15]. While this seems reasonable, the accuracy has not been quantified. A recent alternative approach to this problem, based on second-order perturbation theory in FLRW, can be found here: [16].

Weak lensing of gravitational waves

Gravitational waves (GWs) from supermassive black hole mergers act as precise “standard sirens,” promising to significantly improve upon standard(-isable) candles such as Type-Ia Supernovae (SN1a) (see e.g., [17]). However, weak lensing of the GWs by the intervening dark matter distribution distorts the signal, degrading their use as cosmological distance estimators [18]. A particular problem is that the GW wavelength is comparable to the size of the dark matter halos which produce a portion of the lensing effect. Thus the geometric optics approximation cannot be used to model the expected lensing. By modelling a dark matter halo using a LTB model, and scattering gravitational waves off it using our method, we can hope to quantify the lensing of gravitational waves more accurately.

In this paper, we focus on the linear evolution of perturbations in Gpc-void cosmological models, and present the results for this case (other scenarios shall be analysed elsewhere). To illustrate the performance of the method and the physics of the evolution of perturbations, we concentrate on polar perturbations in a large-scale cosmological void that is asymptotically Einstein-de Sitter (EdS) and that fits the distance-redshift relations given by SN1a observations as well as age data. We show results for several different spherical harmonic frequencies, considering a variety of initial conditions at different locations throughout the void. We find that couplings can only be neglected if one is interested in the matter density growth function on very large scales. Indeed, for this quantity, the evolution without couplings is accurate at the sub-percent level when the large-scale quadrupole is considered. However, the metric perturbation that is the closest to the standard “Newtonian” gravitational potential gets corrections of up to 10% even in this case (which implies percent-level inaccuracies to the 2-point correlation function given in [10]). Furthermore, for higher frequency perturbations even the density contrast suffers more than 10% inaccuracies. This has significant impact on constraints on LTB models which use such approximations for analyses of structure formation.

2 Background evolution

The general unperturbed LTB line element may be written as,


where [5]


with a prime being shorthand for . The angular and radial scale factors, and , respectively, are associated with individual expansion rates


where an overdot denotes . While we only focus on a dust energy-density component in our analysis, we include the cosmological constant in our equations for completeness. The angular component of the expansion then obeys the following Friedmann equation with different curvature, i.e. , on each radial shell:




In these expressions, is the angular Hubble parameter today, a boundary condition related to the matter content within a comoving shell of radius , and we set by convention. The expansion and shear scalars are given by


respectively. The total energy density (including dust and a cosmological constant) may be expressed as


We require the solution (i.e. ) to Equation (4), from which everything else follows. Integrating Equation (4) we find


where is the bang time function.


Cosmic time (Gyr)

Radial distance (Gpc)
Figure 1: The spacetime evolution of selected contrasts in the background dynamics, illustrating the growth of the background void over time. Left: Contrast in the energy density . Centre: Contrast in . Right: Contrast in . Note that scales of the vertical and horizontal axes apply to all such 2D plots in this paper. The values of and (respectively the maximum and minimum of the color scale) can be read at the top of each 2D plot.

In this work we consider a test case to demonstrate the method, and so focus on an open (), dust-only background model () that is asymptotically Einstein-de Sitter, and which is known to comfortably accommodate distances to SN1a (see e.g. [7]). Here we model the total (dimensionless) matter density profile today according to:


where is the length-scale of the inhomogeneity, and = 70 km/s/Mpc 0.23 Gpc is the Hubble constant. The superscript “in” denotes evaluation at . For consistency with an inflationary paradigm, we ignore decaying modes by forcing the bang time function to be uniform throughout space; setting is sufficient. In this case, we may write the solution to Equation (9) in the following parametric form




Figures 1 and 2 show the behaviour of various illustrative quantities in this model for the values of the constants quoted above.

Shear scalar ()

Cosmic time (Gyr)

Radial distance (Gpc)
Figure 2: The spacetime evolution of the shear scalar, . In the vicinity of the origin and far outside the void, the spacetime is effectively FLRW (depicted in green). The presence of propagating modes is more apparent in all the variables here. The maximum and minimum values of the colour scale are given in brackets above.

3 Polar perturbations

Perturbations on a spherically symmetric background are decomposed into spherical harmonics, where tensorial quantities of degree 1 or 2 split into two parities – polar and axial (this is analogous, but not identical, to the usual scalar-vector-tensor split in a FLRW background). In this paper, we restrict our attention to the polar (even parity) sector since this is where the density perturbation is defined. Furthermore, we only consider modes with spherical harmonic index . The modes corresponding to and obey different equations (see [1]) and should be treated separately, although in a similar fashion as far as the numerical integration is concerned. A similar treatment could be straightforwardly applied to the axial (odd parity) sector. Details of the derivation of the perturbation equations we present below can be found in [1]. Let us emphasize here that we are not trying to develop a full analysis of realistic structure formation in a LTB Universe. Rather, we would like to demonstrate that the method we develop can integrate the perturbation equations, allowing us to study a few remarkable features of the evolution of perturbations in a large cosmic void. Note that we have adapted the equations to include the cosmological constant. We have also written them in terms of partial radial derivatives rather than frame derivatives, in readiness for numerical integration.

3.1 Formalism

The general form of polar perturbations to the background LTB metric, in the Regge-Wheeler (RW) gauge, is expanded in spherical harmonics as:

where , and are functions of (independent of due to the spherical symmetry of the background), are the scalar spherical harmonic functions. For notational convenience, an implicit sum over is implied whenever a quantity is multiplied by . The 4-velocity field in the polar sector is given by (indices which are capitals run over , and indices run over )


where and , and , both functions of , are the radial and angular (peculiar) velocities, respectively, is the gauge-invariant metric perturbation, and . The energy-momentum tensor


defines the density contrast . Note that all our perturbation variables, i.e. both metric and fluid perturbations, are automatically gauge-invariant. This is due to the fact that the all perturbations conveniently reduce, in the RW gauge, to the corresponding variables arising from a general gauge (coordinate) transformation (see [20] and references therein).

The 1st-order perturbed Einstein equations for the case reduce to:


The three equations (19), (20) and (21) represent the master equations of our problem, and are the evolution equations for the metric perturbations , and . Knowing these master variables, one can then obtain , and , i.e. the behaviour of the matter perturbations. These are given by


The conservation of the perturbed energy-momentum, i.e. , implies that our solutions must satisfy:


These three equations can thus be seen as constraints that a solution to the previous system, i.e. Eq.’s Equation (19)Equation (21) along with Equation (22)Equation (24), must satisfy. These constraints can be used to test the accuracy of the numerical implementation.

3.2 Weyl curvature

In studying perturbations of a LTB background, it is interesting to consider the evolution of the Weyl curvature tensor, as it describes genuine relativistic effects. In particular the magnetic part of the Weyl tensor is zero in the background and is a gauge-invariant tensor sourced by purely relativistic effects — frame dragging (vector modes in FLRW perturbation theory, for example) and gravitational waves.

In the background is zero, and the only non-zero background parts for are:


The non-zero perturbed parts for the electric and magnetic Weyl tensors are:




respectively. Here the trace () and trace free () parts arise from a decomposition on the 2-sphere. The magnetic part of the Weyl tensor has a parity opposite to the electric part, and carries a bar to denote that the axial part of it appears in the polar equations (and vice versa).

4 Numerical Implementation

We use a method of lines approach [19] to solving the system Equation (19)Equation (21), whereby the spatial domain is discretized by standard finite differences and integrated pointwise in time using a 4th-order Runge-Kutta solver. The system can be recast in terms of dimensionless variables through the following transformations:


where the last equality, motivated by the standard form, defines the conformal time for central observers. Then, using


we have


We also introduce the dimensionless angular peculiar velocity,


and relate the cosmic time to the dimensionless conformal time by


Thus, to evolve the system from some initial time () until today (), we compute the corresponding initial and final values for , i.e. and , respectively, using:


4.1 Discretisation

We introduce discretized time and space coordinates and given by


where and are grid spacings in and , respectively, and and . The factor determines the grid resolution relative to an baseline (with the number of grid points and increased proportionally to cover the same domain). For our system of equations, the Courant-Friedrics-Lewy (CFL) condition requires that


for numerical stability, where . For our simulations we use


Spatial derivatives are calculated using 2nd-order finite difference operators. For some quantity evaluated at time and position at the baseline resolution,


We evaluate the RHS of Eqs. Equation (19)Equation (21) on a slice and evolve forward in time using a standard 4th-order Runge-Kutta integrator. The overall scheme is 2nd-order accurate due to the spatial finite differencing used.

Figure 3: An illustration of the spacetime domain of the model (left), along with the generic form of the initial conditions used (right), according to Eq. Equation (53).

4.2 Initial and boundary conditions

For this work, we have used a generic set of initial conditions for each of the three master variables:


where Gpc is an array of five equally spaced points between and our desired region of interest , and Gpc sets the width of each pulse see the right panel of Fig. 3. The conditions are set at (corresponding to a time Gyr, or a redshift in a fiducial cosmology). These initial data are not meant to represent a realistic physical state of perturbations in the early Universe; they simply allow us to test the method and to extract the physical behavior of perturbations.

Regularity conditions determine the variables in the neighbourhood of the origin according to the prescription of [20]. Near , we require (for all ):


where the hatted variables are all polynomials of even power in . Using


we find that


which vanish at . Fixing the value of all variables to zero at the origin is sufficient for regularity111Note that while our choice of initial conditions are not exactly zero at the origin, they are well below the machine precision.

We require an additional boundary condition at the outer edge of the computational domain, . This boundary condition is necessarily artificial since we do not compactify the spatial coordinate. We place the boundary at a distance from the origin such that it is causally disconnected from the region where the initial perturbations (53) and (54) are set for the duration of the evolution. This prevents artificially reflected signals from influencing the evolution of these perturbations. We can estimate an appropriate distance by tracing null geodesics inward from the outer boundary. Using the background LTB line element, radial null geodesics are given by


where our characteristics approach 45 at large distances in our coordinates. An appropriate value for which is sufficiently removed from the measurement domain is


where is the outer boundary of the domain in which we would like to analyse the behaviour of perturbations between times and . Since we are working in a single spatial dimension, this grid extension to remove outer boundary effects is not overly costly in terms of memory or computation time, though in the future it may be useful to consider a logarithmic radial coordinate. Given that the spacetime in our model is effectively homogeneous above , we choose a conservative region of interest of .

Cosmic time (Gyr)
Figure 4: Amplitude of the measure of error over time for each constraint equation, in the case of an initialised , for , 10, 200, 1000. Starting from a reference resolution of , we include curves of double () and four times () the resolution. We multiply the higher resolution errors by factors of 4 and 16, respectively, indicating 2nd-order convergence where the curves line up. We see clear convergence in each case, though for , the resolution must be pushed to at least for the curves to align satisfactorily. The sudden drop in error seen in some plots around 4 Gyr is associated with the exiting of the initial pulses from the measurement domain.

5 Convergence Tests

To establish the accuracy of the discretization, we carry out a standard convergence test. We check the 2nd-order convergence empirically by carrying out a series of runs of the same initial data at successively doubled resolution, corresponding to , , and in Eqs. Equation (47) and Equation (48). The rate of convergence for a variable is given by


where is the -norm on the fixed-resolution grid,


with () the number of spatial grid-points stored for analysis within the range . We use the following dimensionless measure to quantify how well the constraints are satisfied:


where , is one of Equation (25)Equation (27), and we estimate via a centered difference, i.e.


For all of our evolution variables and constraints, we observe the expected 2nd-order convergence rate (). Examples are plotted in Fig. 4, which shows how well the constraint equations perform for various multipole moments in the case of an initial . Using a reference resolution of , we include curves of double () and four times () the resolution, multiplying each by a factor of 4 and 16, respectively. When the curves line up, then , i.e. we obtain second order convergence as expected. It is clear that, in the case , one must go to a resolution of at least to obtain the correct convergence. We have used a base resolution of Gpc throughout, typically with .

6 Results

6.1 Evolution of perturbations

To study the evolution of perturbations in a LTB cosmological void, we concentrate, for the sake of illustration, on the evolution of the perturbation variables for the spherical modes and . For each , we consider three distinct cases: we initialise the profile of any one of , and according to Equation (53) while setting the others to zero, and apply Equation (54) to all variables (except since it satisfies a 1st-order PDE). Here, cases 1, 3 and 5 correspond to , and cases 2, 4 and 6 correspond to .




Cosmic time




Radial distance
Figure 5: Spacetime evolution of each of the master variables in the various cases considered. In cases 1 and 2, excites both and to about the sub-percent level. The propagating modes resulting from is visible in , or alternatively the trace-free part of the magnetic Weyl tensor Equation (35), thus clearly showing relativistic degrees of freedom at work. In cases 3 and 4, while an initial decays away quickly due to the Hubble friction, it still manages to excite the other two variables, albeit to a low level. The final two cases, 5 and 6, are rather interesting: an initial generated inside a void excites a relatively significant amount of and . The presence of propagating modes is more apparent in all the variables here. The maximum and minimum values of the colour scale are indicated respectively in brackets above and below each 2D plot.

The evolution of each of the variables is presented in Fig. 5 for each of the corresponding cases. The resolutions used in all of the plots are typically in the range , .

Cases 1 and 2:

In these cases, we initialise , and set initially. We clearly see the “bleeding” of the modes due to the coupling.

 Time evolution             Profile today  Time evolution             Profile today
Cosmic time (Gyr)         Distance (Gpc) Cosmic time (Gyr)         Distance (Gpc)
Figure 6: Temporal and spatial slices through the spacetime evolution observed in the top panel of Fig. 6. The variable is clearly unaffected on the outskirts of the void in which the spacetime is almost Einstein-de Sitter. Since and are relatively sub-percent in amplitude, on each radial shell more or less behaves as expected in an open, dust-dominated FLRW model, i.e. decays with time.

On the 2D plot, Fig. 5, behaves like a propagating degree of freedom, evolving along the characteristics of the spacetime, and radiating energy away from each pulse. The behaviour of is more difficult to qualitatively describe because it is a mixture of frame dragging and gravitational wave degrees of freedom – the combination of non-propagating decay with some radiation can be seen in the figures. It is proportional to the trace-free part of the magnetic Weyl curvature Equation (35), and thus represents a true relativistic degree of freedom. Then, follows a standard evolution throughout the spacetime: staying constant around the EdS region, while decaying faster deep inside the void.

The top panel of Fig. 6 presents the profile of today for these cases, as well as its time evolution along selected radii. As expected, remains constant in the outer, quasi-FLRW regions of the void, given that it essentially satisfies the Bardeen equation there. Deep inside the void, decreases for the most part as the usual Bardeen potential would in an open FLRW dust model, but there is evidence of influence from and at least at the sub-percent-level, as can be seen by the amplitudes of the latter in the middle and bottom panels of Fig. 6; see also section 6.2 for a discussion on the importance of the couplings.

We show the spacetime configuration of in Fig. LABEL:Tab:DerivedQuantities: its growth appears to follow the peaks were is concentrated, suggesting that the tiny and generated by the dynamics have little impact on the profile of density perturbations. Also included in Fig. LABEL:Tab:DerivedQuantities is the resulting radial peculiar velocity and the trace-free part of the electric Weyl curvature.

Cases 3 and 4:

Here, we initialise , and set initially.

From the middle panels of Fig. 5 we see that decays very quickly in fact, roughly proportional to along the peaks from where it is initially located, while sourcing and . As expected, is very well described as a propagating degree of freedom, but one also sees that the sourced has a propagating component that follows the characteristics of the background spacetime and escapes the void.

 Time evolution             Profile today  Time evolution             Profile today
Cosmic time (Gyr)         Distance (Gpc) Cosmic time (Gyr)         Distance (Gpc)
Figure 7: Temporal and spatial slices through the spacetime evolution observed in the central panel of Fig. 6. The variable decays roughly in the quasi-FLRW regions, but decreases more quickly deep inside the void due to the faster expansion rate there. Compared to the initial amplitude of , and remains sub-percent in magnitude.
 Time evolution             Profile today  Time evolution             Profile today
Cosmic time (Gyr)         Distance (Gpc) Cosmic time (Gyr)         Distance (Gpc)
Figure 8: Temporal and spatial slices through the spacetime evolution observed in the bottom panel of Fig. 6. Here it is clear that the evolution of is dominated by a propagating mode. It is interesting to note that the level of generated is of a similar order of magnitude than the initial , and again showing a slower ( 25%) growth rate inside the void compared to the outskirts. The variable is also produced to a significant proportion early on, but nevertheless decays quickly with time.

The middle panel of Fig. 7 presents the profile of today for these cases, including the time evolution along selected radii. It’s clear that decays, for the most part, approximately as . (In the FLRW limit this would be a pure vector mode with this exact decay rate.) The greater decay in seen in the central regions of the void can be attributed to the faster expansion rate there. The top and bottom panels of Fig. 7 show the profiles for the other variables, and .

We also show the spacetime configuration of in Fig. LABEL:Tab:DerivedQuantities. Remarkably, the density contrast generated initially by the presence of the perturbation also decays very rapidly at the peak locations, except deep inside the void (first few peaks) where it starts to grow at later times (much more at small angular scales than at large ones). This is associated with the potential deepening in this region at the same time: the decay of into is associated with a growth of structure deep within the void.

=2 =10
            % diff.                     Profile today            % diff.                    Profile today

          Time (Gyr)                Distance (Gpc)          Time (Gyr)               Distance (Gpc)
 Initial profile             % diff.            Profile today  Initial profile             % diff.            Profile today

    Distance                Time                Distance     Distance                Time                Distance
Figure 10: Comparison of coupled to uncoupled runs, for cases 1 and 2. On the largest scales (), and deep within the void (), (top panel) is enhanced by 10% when the coupling is present, while (bottom panel; normalised to its maximum value in space at the initial time, ) is enhanced by 1%. As we approach the outskirts of the void, the differences are sub-percent, as expected.
Cases 5 and 6:

Here, we initialise , and set initially.

According to the bottom panels of Fig. 5, and the generated propagate to the outskirts of the void along the characteristics of the background, resulting in the localised generation of the potential , and the associated growth of density perturbations, as is shown in Fig. LABEL:Tab:DerivedQuantities.

The bottom panel of Fig. 8 presents the profile of today for these cases, as well as the time evolution along selected radii, while the profiles for the other variables, and , are shown in the top and middle panels.

All of these cases demonstrate that , and are much more difficult to interpret than on a FLRW background. As emphasised in [1], they are mixtures of scalar, vector and tensor modes and therefore their coupling is an essential ingredient of first order perturbation theory around a LTB background: in principle, they cannot be treated as separate, independent modes that describe different physical aspects of perturbations. In the next subsection, we compare the behaviours of the fully coupled perturbation system to cases where they are decoupled “by hand” as has been done before in various ways to simplify the analysis [12, 13, 10]. In particular, we would like to analyse the errors in and when and are neglected, not only for initial data but also all during the evolution of the system.

6.2 Coupled vs uncoupled dynamics

In this section, we quantify the errors induced when assuming that the coupling of to and is negligible by considering models which initialise only. We compare these to cases where Eqs. Equation (19), Equation (20) and Equation (21) are solved retaining terms with no coupling between and , that is, by solving the reduced system:

Percentage errors on and when neglecting the full coupling
   1st peak            3rd peak          5th peak    1st peak            3rd peak          5th peak
Multipole moment