Evolution of linear perturbations in spherically symmetric dust spacetimes
Abstract
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 gigaparsec void, the evolution of the density contrast is well approximated by the previously studied decoupled evolution only for very largescale 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.
pacs:
98.80k,98.80.Jk,98.65.Dx,04.25.D,04.25.Nxams:
35Q75,8308,83F05,83C251 Introduction
We present results from a numerical code implementing a new method that solves the first order gaugeinvariant linear perturbation equations in a Lemaî treTolmanBondi (LTB) background. LTB models are spherically symmetric but inhomogeneous dust solutions of the Einstein field equations. Compared to a FriedmanLemaîtreRobertsonWalker (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 nonlinear models of Newtonian systems, or relativistically but only in the linear regime. This leaves an important area unexplored: nonlinear relativistic aspects of structure formation – see, e.g., [2] for a review. Certain aspects of this have started to be taken into account in Nbody 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 nonlinear background — either an overdensity 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 gigaparsecs 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 secondorder 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 TypeIa 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 Gpcvoid 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 largescale cosmological void that is asymptotically Einsteinde Sitter (EdS) and that fits the distanceredshift 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 subpercent level when the largescale 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 percentlevel inaccuracies to the 2point 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,
(1) 
where [5]
(2) 
with a prime being shorthand for . The angular and radial scale factors, and , respectively, are associated with individual expansion rates
(3) 
where an overdot denotes . While we only focus on a dust energydensity 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:
(4) 
where
(5) 
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
(6)  
(7) 
respectively. The total energy density (including dust and a cosmological constant) may be expressed as
(8) 
We require the solution (i.e. ) to Equation (4), from which everything else follows. Integrating Equation (4) we find
(9) 
where is the bang time function.
Cosmic time (Gyr) 



Radial distance (Gpc) 
In this work we consider a test case to demonstrate the method, and so focus on an open (), dustonly background model () that is asymptotically Einsteinde 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:
(10) 
where is the lengthscale 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
(11)  
(12) 
where
(13)  
(14)  
(15) 
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) 
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 scalarvectortensor 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 ReggeWheeler (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 4velocity field in the polar sector is given by (indices which are capitals run over , and indices run over )
(17) 
where and , and , both functions of , are the radial and angular (peculiar) velocities, respectively, is the gaugeinvariant metric perturbation, and . The energymomentum tensor
(18) 
defines the density contrast . Note that all our perturbation variables, i.e. both metric and fluid perturbations, are automatically gaugeinvariant. 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 1storder perturbed Einstein equations for the case reduce to:
(19)  
(20)  
(21)  
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
(22)  
(23)  
(24) 
The conservation of the perturbed energymomentum, i.e. , implies that our solutions must satisfy:
(25)  
(26)  
(27) 
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 gaugeinvariant 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 nonzero background parts for are:
(28)  
(29) 
The nonzero perturbed parts for the electric and magnetic Weyl tensors are:
(30)  
(31) 
(32)  
(33) 
and
(34)  
(35) 
respectively. Here the trace () and trace free () parts arise from a decomposition on the 2sphere. 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 4thorder RungeKutta solver. The system can be recast in terms of dimensionless variables through the following transformations:
(36)  
(37) 
where the last equality, motivated by the standard form, defines the conformal time for central observers. Then, using
(38)  
(39) 
we have
(40)  
(41)  
(42)  
(43) 
We also introduce the dimensionless angular peculiar velocity,
(44) 
and relate the cosmic time to the dimensionless conformal time by
(45) 
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:
(46) 
4.1 Discretisation
We introduce discretized time and space coordinates and given by
(47)  
(48) 
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 CourantFriedricsLewy (CFL) condition requires that
(49) 
for numerical stability, where . For our simulations we use
(50) 
Spatial derivatives are calculated using 2ndorder finite difference operators. For some quantity evaluated at time and position at the baseline resolution,
(51)  
(52) 
We evaluate the RHS of Eqs. Equation (19)Equation (21) on a slice and evolve forward in time using a standard 4thorder RungeKutta integrator. The overall scheme is 2ndorder accurate due to the spatial finite differencing used.
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:
(53)  
(54) 
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 ):
(55) 
where the hatted variables are all polynomials of even power in . Using
(56) 
we find that
(57)  
(58)  
(59)  
(60) 
which vanish at . Fixing the value of all variables to zero at the origin is sufficient for regularity^{1}^{1}1Note 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
(61) 
where our characteristics approach 45 at large distances in our coordinates. An appropriate value for which is sufficiently removed from the measurement domain is
(62) 
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) 
5 Convergence Tests
To establish the accuracy of the discretization, we carry out a standard convergence test. We check the 2ndorder 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
(63) 
where is the norm on the fixedresolution grid,
(64) 
with () the number of spatial gridpoints stored for analysis within the range . We use the following dimensionless measure to quantify how well the constraints are satisfied:
(65) 
where , is one of Equation (25)Equation (27), and we estimate via a centered difference, i.e.
(66) 
For all of our evolution variables and constraints, we observe the expected 2ndorder 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 1storder PDE). Here, cases 1, 3 and 5 correspond to , and cases 2, 4 and 6 correspond to .
Initialising 




CASE 1  CASE 2  

Initialising 



CASE 3  CASE 4  
Initialising 



CASE 5  CASE 6  
Radial distance 
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.
CASE 1 CASE 2 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 Einsteinde Sitter. Since and are relatively subpercent in amplitude, on each radial shell more or less behaves as expected in an open, dustdominated 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 nonpropagating decay with some radiation can be seen in the figures. It is proportional to the tracefree 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, quasiFLRW 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 subpercentlevel, 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 tracefree 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.
CASE 3 CASE 4 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 quasiFLRW regions, but decreases more quickly deep inside the void due to the faster expansion rate there. Compared to the initial amplitude of , and remains subpercent in magnitude. CASE 5 CASE 6 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 subpercent, 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:
(67)  
(68)  
Percentage errors on and when neglecting the full coupling  
1st peak 3rd peak 5th peak  1st peak 3rd peak 5th peak  






Multipole moment 