Incorporating curved geometry in cosmological simulations

TL;DR

This paper introduces a fully relativistic framework for cosmological simulations that incorporate curved spatial geometry, enabling more accurate modeling of large-scale structures and observables in a universe with spatial curvature.

Contribution

The authors develop a novel relativistic simulation framework that embeds a curved spacetime region into a flat exterior, allowing for consistent modeling of curvature effects in cosmological simulations.

Findings

Framework successfully models curved spacetime in simulations

Boundary conditions are effectively managed via embedding in flat space

Applicable to Newtonian codes with minimal modifications

Abstract

Spatial curvature is one of the fundamental cosmological parameters that is routinely constrained from observations. The forward modelling of observations, in particular of large-scale structure, often relies on large cosmological simulations. While the so-called separate universe approach allows one to account for the effect of curvature on the expansion rate in small sub-volumes, the non-Euclidean geometry is harder to accommodate. It becomes important when observables are computed over large distances, e.g. when photons travel to us from high redshift. Here we present a fully relativistic framework to run cosmological simulations for curved spatial geometry. The issue of consistent boundary conditions is solved by embedding a spherical cap of the curved spacetime into a hole within a flat exterior, where it can undergo free expansion. The geometric nature of gravity is made explicit in our framework, allowing for a consistent forward modelling of observables inside the curved patch. Our methodology would also work with any Newtonian code to a good approximation, requiring changes only to the initial conditions and post-processing.

Equations

$$d{s}^2=-e^{2\psi }d{t}^2+a^2(t)e^{-2\varphi }\left[d{r}^2+r^{2}d{\Omega }^2\right],$$

$$4{\nabla }^2\varphi +2{\left(\nabla \varphi \right)}^2+6{\left({\partial }_t\varphi -H\right)}^2{e}^{-2\psi }=16\pi G\rho ,$$

$$d{s}^2=-d{t}^2+\frac{{\left({\partial }_{r}R\right)}^2}{1+2E(r)}d{r}^2+R^{2}d{\Omega }^2,$$

$$R(t,r)$$

$${\left(\eta -\sin \eta \right)}^{2/3}$$

$$GM(r)=\frac{3}{2}{H}^2{a}^3{\int }_0^r{\tilde{r}}^2\left[1+\delta (\tilde{r})\right]d\tilde{r},$$

$$\delta (r)=\{\begin{array}{ll}{\delta }_1& r<r_1=r_2{\left(1+{\delta }_1\right)}^{-1/3}\\ -1& r_1<r<r_2\\ 0& r>r_2.\end{array}$$

$$f(r)=\{\begin{array}{ll}{\delta }_1& r<r_1\\ \frac{r_2^3}{r^3}-1& r_1<r<r_2\\ 0& r>r_2.\end{array}$$

$$E(r)=-\frac{5}{6}{H}_{\mathrm{in}}^2{a}_{\mathrm{in}}^{2}f(r)\left[1-\frac{1}{7}f(r)\right]r^2,$$

$$R(t,r)=a(t)\left[1+\frac{1}{3}\left(1-\frac{a(t)}{a_{\mathrm{in}}}\right)f(r)\left(1-\left(\frac{1}{3}-\frac{1}{7}\frac{a(t)}{a_{\mathrm{in}}}\right)f(r)\right)\right]r.$$

$${\tilde{H}}^2+\frac{K}{a_{\mathrm{in}}^2}=H_{\mathrm{in}}^2\left(1+{\delta }_1\right),$$

$$t^{(\mathrm{syn})}$$

$$r^{(\mathrm{syn})}$$

$$u_{\mu }^{\text{(P)}}=u_{\nu }^{\text{(syn)}}\frac{d{x}_{\text{(syn)}}^{\nu }}{d{x}_{\text{(P)}}^{\mu }},$$

$${\tilde{\Omega }}_K$$

$${\tilde{\Omega }}_{\mathrm{m}}$$

$${\tilde{\Omega }}_{\Lambda }$$

$$\tilde{H}$$

$$t_0^{\text{(syn)}}-t_{\mathrm{in}}^{\text{(syn)}}={\tilde{H}}_0^{-1}{\int }_0^{{\tilde{z}}_{\mathrm{in}}}\frac{d{z}^{\prime }}{1+z^{\prime }}{\left[{\tilde{\Omega }}_{\mathrm{m}}^0(1+z^{\prime }{)}^3+{\tilde{\Omega }}_K^0(1+z^{\prime }{)}^2+{\tilde{\Omega }}_{\Lambda }^0\right]}^{-1/2}.$$

$$H_{\mathrm{in}}\simeq \tilde{H}\left(1+\frac{1}{3}{\delta }_1+\frac{1}{21}{\delta }_1^2\right)\simeq \tilde{H}\left(1-\frac{1}{5}{\tilde{\Omega }}_K-\frac{8}{175}{\tilde{\Omega }}_K^2\right).$$

$${\Omega }_{\mathrm{m}}$$

$${\Omega }_{\Lambda }$$

$$t_0^{\text{(syn)}}-t_{\mathrm{in}}^{\text{(syn)}}\underset{\text{correction term}}{\underbrace{-T_1(t_{\mathrm{in}}^{\text{(syn)}},0)-T_2(t_{\mathrm{in}}^{\text{(syn)}},0)}}={\int }_{t_{\mathrm{in}}^{\text{(P)}}=t_{\mathrm{in}}^{\text{(syn)}}}^{t_0^{\text{(P)}}}{e}^{\psi }d{t}^{\text{(P)}}.$$

$${\partial }_r^2{\varphi }_1+\frac{2}{r}{\partial }_r{\varphi }_1-3H_{\mathrm{in}}^2{a}_{\mathrm{in}}^2{\varphi }_1=\frac{3}{2}{H}_{\mathrm{in}}^2{a}_{\mathrm{in}}^2\left[3{\varphi }_1-{\partial }_r\delta r_1-\frac{2}{r}\delta r_1\right],$$

$$\delta r_1=-\frac{2}{3H_{\mathrm{in}}^2{a}_{\mathrm{in}}^2}{\partial }_r{\varphi }_1+\frac{5}{r^2}{\int }_0^r{\tilde{r}}^2{\varphi }_1(\tilde{r})d\tilde{r}.$$

$${\partial }_r^2{\varphi }_1+\frac{2}{r}{\partial }_r{\varphi }_1-3H_{\mathrm{in}}^2{a}_{\mathrm{in}}^2{\varphi }_1=\frac{3}{2}{H}_{\mathrm{in}}^2{a}_{\mathrm{in}}^2\left[3{\varphi }_1-{\partial }_r\delta r_1-\frac{2}{r}\delta r_1+\frac{\delta m_1}{m}\right],$$

$$\delta r_1=-\frac{2}{3H_{\mathrm{in}}^2{a}_{\mathrm{in}}^2}{\partial }_r{\varphi }_1+\frac{r}{3}\frac{\delta m_1}{m}+\frac{5}{r^2}{\int }_0^r{\tilde{r}}^2{\varphi }_1(\tilde{r})d\tilde{r}.$$

$$\frac{\delta m_1}{m}=-\frac{15}{r_2^3}{\int }_0^{r_2}{r}^2{\varphi }_1(r)dr.$$

$$\begin{gathered} {\partial }_r^2{\varphi }_2+\frac{2}{r}{\partial }_r{\varphi }_2-3H_{\mathrm{in}}{a}_{\mathrm{in}}^2{\partial }_t{\varphi }_2-3H_{\mathrm{in}}^2{a}_{\mathrm{in}}^2\left({\psi }_2-{\varphi }_1^2\right)+2{\varphi }_1\left({\partial }_r^2{\varphi }_1+\frac{2}{r}{\partial }_r{\varphi }_1\right)-\frac{1}{2}{\left({\partial }_r{\varphi }_1\right)}^2 \\ =\frac{3}{2}{H}_{\mathrm{in}}^2{a}_{\mathrm{in}}^2[-{\partial }_r\delta r_2-\frac{2}{r}\delta r_2+3{\varphi }_2+\frac{\delta m_2}{m}-\left(3{\varphi }_1+\frac{\delta m_1}{m}\right)\left({\partial }_r\delta r_1+\frac{2}{r}\delta r_1\right) \\ +3{\varphi }_1\frac{\delta m_1}{m}+\frac{3}{r^2}{\left(\delta r_1\right)}^2+\frac{2}{r}\delta r_1{\partial }_r\delta r_1+{\left({\partial }_r\delta r_1\right)}^2+\frac{1}{2}{v}_1^2], \end{gathered}$$

$$\psi$$

$$\varphi$$

$$\frac{{\left({\partial }_{r}R\right)}^2}{1+2E(r)}$$

$$R^2$$

$$b_1$$

$$y_1$$

Full text

1 Introduction

Spatial flatness is a critical prediction of inflation [1] that should be tested and not assumed from the outset. Observations of anisotropies in the cosmic microwave background (see, e.g. [2] for early results and [3, 4] for recent data) put strong constraints on curvature which are, however, somewhat model dependent. Some analyses even suggest a preference for spatial curvature; see [5] for an example. Largely independent constraints can be obtained from distance measurements in the late Universe and certain probes of large-scale structure. For example, recent measurements of baryon acoustic oscillations (BAO) by the Dark Energy Spectroscopic Instrument (DESI) constrain the curvature parameter, $\Omega_{K}$, to $0.025\pm 0.041$ at $68\%$ confidence [6]. It has been argued that even a curvature as small as $0.2\%$, corresponding to a radius of curvature of $21H_{0}^{-1}$, would be important for cosmological constraints on the neutrino mass [7]. A conclusive detection of spatial curvature would have profound implications for the origin of our Universe of course, shattering the prevailing paradigm of inflation and offering a rare glimpse at the conditions from which it emerged.

With the wealth of observational data from surveys such as DESI, Euclid [8, 9], or the Vera C. Rubin observatory [10], obtaining direct constraints that are robust below the percent level will require careful modelling. In the context of large-scale structure, the forward modelling typically relies on sophisticated simulations of the gravitational dynamics and the kinematical projection effects on the past light cone of our observations. Such simulations are employed, for example, to estimate covariances [11], to train emulators that can be used for parameter inference [12, 13], or for simulation-based inference [14, 15].

In the late Universe, curvature has two main effects on large-scale structure probes. First, it affects the expansion history and therefore the growth rate. This effect is relatively straightforward and has been studied with the so-called separate universe approach [16, 17]. Since matter only moves a few megaparsecs in one Hubble time, distances can be considered approximately Euclidean on the scales of interest for structure formation. Simulations of small Euclidean sub-volumes that account for the correct expansion history on a curved background therefore capture the evolution of matter well.

The second effect concerns the non-Euclidean distance and volume when considering large scales. These are the scales probed by the photons as they travel cosmological distances. Intuitively, given a standard ruler such as the BAO scale, a curved geometry implies that the circumference of a circle around the observer is either larger (for an open Universe) or smaller (for a closed Universe) than $\pi$ times the geometric diameter when measured with the same ruler. Similar arguments hold for the surface of a sphere or for the physical volume of a redshift bin. This poses the question of how the sub-volumes comprising the separate universe simulations could be stitched together consistently to build the full light cone of a curved Universe.

Further interesting aspects of curvature concern the early Universe. Since standard inflation would always lead to a flat spatial geometry, it cannot be invoked to predict the primordial power spectrum, and thus we need to modify our prescription for setting initial conditions. Even if some inflation occurred, perturbation modes close to the curvature scale would always remain sensitive to pre-inflationary physics. In addition, the basis of mode functions that describe the perturbations also depends on the spatial geometry, with the spectrum becoming discrete in a Universe that is spatially closed. These aspects become important at exactly the scales we are interested in, and using them to gain insights into the origin of the Universe is certainly very compelling. However, from the point of view of simulations, they are relevant only for setting initial conditions of matter perturbations, which is not the main focus of this work.

In most simulation frameworks, in particular $N$-body simulations, accounting for a curved spatial geometry is not straightforward. The reason is that boundary conditions need to be imposed, and almost all frameworks have settled on periodic boundaries. Beyond the scale of periodicity, the volume therefore necessarily behaves like in a flat geometry. Although changing boundary conditions could be done in principle (see, e.g. [18] for an interesting example), it is often not desired because too much of the established code infrastructure depends on this assumption. We will therefore work within the constraints of periodic boundaries, meaning that we must incorporate curvature somehow below the scale of periodicity.

The remainder of this paper is structured as follows. In Sec. 2 we present our methodology, which is essentially based on translating a specific Lemaître–Tolman–Bondi model into the coordinates used in simulations. In Sec. 3 we show some first numerical results to validate our framework. We conclude in Sec. 4. Some details of our calculations are relegated to App. A. As a side product of our work, we present in App. B the Schwarzschild metric in a new coordinate system motivated by our setup.

2 Methodology

The basic idea of our approach is to embed a patch of a curved Friedmann–Lemaître–Robertson–Walker (FLRW) spacetime, which shall contain any observations we wish to model, into an exterior solution that is spatially flat, so that periodic boundaries can be employed. The junction conditions for joining different solutions of Einstein’s equations have been worked out in full generality in a work by Israel [19]. Using these prescriptions, one finds that any spherical patch of a FLRW spacetime with pressureless matter and arbitrary curvature, truncated at some fixed comoving radius, can be embedded into the Schwarzschild solution. The resulting spacetime is an exact solution to Einstein’s equation describing a uniform ball of matter undergoing free expansion in vacuum.

This solution is not yet what we want since Schwarzschild only becomes flat asymptotically at spatial infinity. However, we can repeat the process once more and embed the Schwarzschild solution (which embeds the curved FLRW model) into an outer flat FLRW model. We can then fit the entire system into a cubic simulation box with periodic boundary conditions. The setup is illustrated in Fig. 1. A similar idea was originally proposed by Einstein & Straus [20] to model an expanding universe filled with stars, and we emphasise once again that this construction is exact for dust-like matter. For the second matching radius to be larger than the radius of the curved interior FLRW patch, the curvature of the latter has to be positive, modelling a closed Universe. For small values of the curvature parameter, the two matching radii approach the same value, merging into one as the curvature goes to zero and the solution simply becomes a continuous flat FLRW spacetime. Open FLRW patches can also be embedded, but one needs to add some compensating mass shell inside the vacuum region to obtain a plausible sequence of matching radii. In this work, we avoid this minor complication and focus on embedding closed FLRW patches, but we keep our derivations general in most places so that the open case can be treated with minimal adjustments.

In Fig. 1, we also indicate two possible choices to place an observer inside the curved patch. Observer A is at the center of the patch, observing a curved Universe in all directions up to the inner matching radius $r_{1}$. While this might be the most obvious choice, we have to keep in mind that simulations always need to trade volume against resolution (for a fixed number of particles). If the observations we want to model go to some maximum redshift $z_{\mathrm{max}}$, this determines the minimum value of $r_{1}$ and thus puts a corresponding limit on the resolution. However, if the observations cover only a part of the sky, we may place the observer off centre to allow a larger redshift coverage or, equivalently, a smaller $r_{1}$ for the same $z_{\mathrm{max}}$. This situation is illustrated with observer B in the figure, who sits close to the perimeter of the curved patch and observes a cone-shaped field of view towards the far side of the patch.

The remainder of this section deals with the question of how to set up this system precisely. We start by noting that for pressureless matter, often called pure dust, the entire spacetime solution is known analytically, and the metric in matter-comoving synchronous coordinates is given by the Lemaître–Tolman–Bondi (LTB) metric (see [21] for details). We shall set initial conditions deep in the matter era, so this solution is a valid starting point. However, simulations are not performed in comoving synchronous coordinates, so we need to work out a suitable coordinate transformation. We choose to work in the Poisson gauge adapted to the exterior flat FLRW model, where the metric takes the form

[TABLE]

where $\psi$ and $\phi$ are the two Bardeen potentials and $a(t)$ is the scale factor of the exterior FLRW that, from the point of view of the simulation code, is the relevant cosmology. We use units where $c=1$. For most of the discussion presented here we do not consider matter perturbations inside the curved domain yet, although we shall get to that at the end. For simplicity, we therefore ignore vector and tensor perturbations and initially consider $\phi$ and $\psi$ as functions of $t$ and $r$ only.

To understand how spatial curvature is incorporated in this picture, it is instructive to consider the Hamiltonian constraint. In any foliation of spacetime, it takes the form ${}^{3}\mathcal{R}-K_{ij}K^{ij}+K^{2}=16\pi G\rho$, where ${}^{3}\mathcal{R}$ is the curvature of the three-dimensional constant-time hypersurface, $K_{ij}$ is the extrinsic curvature tensor describing the embedding into the four-dimensional manifold, $K$ is its trace, and $\rho$ is the matter density measured along the timelike hypersurface normal (see, e.g. [22] for details). For the Poisson gauge, the Hamiltonian constraint becomes

[TABLE]

where $H=\partial_{t}\ln a$ and $\nabla$ is the covariant derivative on the three-dimensional hypersurface, i.e. $\nabla^{2}\phi=a^{-2}e^{2\phi}\left[\partial_{r}^{2}\phi+\frac{2}{r}\partial_{r}\phi-\left(\partial_{r}\phi\right)^{2}\right]$ and $\left(\nabla\phi\right)^{2}=a^{-2}e^{2\phi}\left(\partial_{r}\phi\right)^{2}$. The first two terms of Eq. (2.2) represent the three-curvature ${}^{3}\mathcal{R}$, whereas the third term comes from the extrinsic curvature. In the flat exterior FLRW solution, the two potentials are set to zero and the Hamiltonian constraint reduces to the Friedmann equation. In the curved patch, the potentials are non-zero and $\phi$ gives rise to spatial curvature. As a minor point, note that spacetime is sliced in the rest frame of the exterior FLRW and hence the constant-time hypersurfaces are slightly tilted with respect to the homogeneous slices of the closed patch. This subtlety leads to a slightly non-uniform curvature on constant-time slices (beyond linear order), as well as to some gravitational time delays, which we will have to deal with later.

2.1 The Lemaître–Tolman–Bondi solution

The LTB metric,

[TABLE]

describes an isotropic, but not homogeneous, spacetime in comoving synchronous coordinates for a pressureless matter source. Einstein’s equations yield the parametric solutions

[TABLE]

for the two functions $E(r)$ and $R(t,r)$. Following [23], we set the mass profile to

[TABLE]

where

[TABLE]

Here, $\delta_{1}$ is a uniform initial density contrast inside the region $r<r_{1}$ whose relation with curvature is made explicit later. We therefore have $GM(r)=\frac{1}{2}H^{2}a^{3}r^{3}\left[1+f(r)\right]$ with

[TABLE]

Note that the scale factor $a$ and the associated Hubble function $H$ still refer to the flat outer FLRW region, where the synchronous and Poisson gauges coincide.

At first order in perturbation theory, the gauge transformation between the comoving synchronous gauge and the Poisson gauge is well known and has been applied to the same problem in [23]. Here we wish to extend this calculation to second order. The reason for this is two-fold. First, it allows us to set initial conditions more accurately, making our method more robust when curvature is large. Second, we are ultimately interested in the coupling between perturbations and curvature. If we count curvature itself as one order in perturbation theory, then this coupling appears at second order, and is therefore potentially as important as second-order corrections from curvature itself.

Proceeding as in [23], we rescale the coordinate $r$ so that at some initial time $t_{\mathrm{in}}$ we have $R(t_{\mathrm{in}},r)=a(t_{\mathrm{in}})r\equiv a_{\mathrm{in}}r$. Eliminating $\eta$ perturbatively to second order in $f(r)$, the function $E(r)$ that appears in the LTB metric then reads

[TABLE]

The corresponding solution for $R(t,r)$ reads

[TABLE]

Given the form of the LTB metric, Eq. (2.3), and observing the fact that $f(r)$ is constant inside $r_{1}$, we recover the metric of a curved FLRW spacetime for $r<r_{1}$ and can immediately read off the curvature parameter $K=\frac{5}{3}H_{\mathrm{in}}^{2}a_{\mathrm{in}}^{2}\delta_{1}\left(1-\frac{1}{7}\delta_{1}\right)$. On the initial hypersurface, the curved Friedmann solution has a Hubble rate $\tilde{H}$ that satisfies

[TABLE]

and hence $\tilde{\Omega}_{\mathrm{m}}\simeq 1+\frac{5}{3}\delta_{1}\left(1+\frac{11}{21}\delta_{1}\right)$, $\tilde{\Omega}_{K}\simeq-\frac{5}{3}\delta_{1}\left(1+\frac{11}{21}\delta_{1}\right)$. Here and in the following, we use a tilde to indicate cosmological parameters in the curved patch. The expansion rate is somewhat reduced at $\tilde{H}\simeq H_{\mathrm{in}}\left(1-\frac{1}{3}\delta_{1}+\frac{4}{63}\delta_{1}^{2}\right)$.

2.2 Coordinate transformation

Up to second order in perturbation theory, the relation between synchronous gauge and Poisson gauge is written as

[TABLE]

To avoid clutter, labels (syn) and (P) will only be used when necessary, that is, when two different coordinate types appear in the same equation. Explicit expressions for the coordinate transformations $T_{1}$, $T_{2}$, $L_{1}$, $L_{2}$, and the resulting expressions for $\phi$ and $\psi$ are derived in App. A. To obtain the initial density, we can either work out the gauge transformation of the stress-energy tensor or simply evaluate the Hamiltonian constraint (2.2) given the potentials. The peculiar velocity can be easily inferred by noting that

[TABLE]

with $u_{\nu}^{\text{(syn)}}=-(1,0,0,0)$. Our simulation code uses the canonical momentum as the phase-space coordinate for the particles, given by $q_{i}=mu_{i}^{\text{(P)}}$.

2.3 Cosmological parameters in the simulation

With the relation between the Bardeen potentials and the curved background cosmology in hands, we can initialise the simulation by providing the cosmological parameters and initial redshift as seen by the observer in the curved Universe. These need to be translated into a corresponding set of parameters for the exterior flat FLRW spacetime, as this is the one that is evolved from the point of view of the simulation code. The first step is to obtain the cosmological parameters at the initial time where the coordinate transformation is applied. Working in the curved patch, we have

[TABLE]

where the label [math] indicates parameters evaluated at present time and we assume that the density of radiation can be neglected. Furthermore, we choose $\tilde{z}_{\mathrm{in}}$ so that $\tilde{\Omega}_{\Lambda}$ is sufficiently small and we can use the pure dust solution to compute the coordinate transformations on the initial hypersurface. We also note that the proper look-back time in matter-synchronous time is given by

[TABLE]

We can now fix the initial cosmological parameters of the flat model. First, we determine the auxiliary variable $\delta_{1}$ by inverting the relation with $\tilde{\Omega}_{K}$, finding $\delta_{1}\simeq-\frac{3}{5}\tilde{\Omega}_{K}\left(1+\frac{11}{35}\tilde{\Omega}_{K}\right)$. Using Eq. (2.11) we then have

[TABLE]

Relating cosmological parameters to proper densities we have $\tilde{H}^{2}\tilde{\Omega}_{\mathrm{m}}=H_{\mathrm{in}}^{2}\Omega_{\mathrm{m}}\left(1+\delta_{1}\right)$ and $\tilde{H}^{2}\tilde{\Omega}_{\Lambda}=H_{\mathrm{in}}^{2}\Omega_{\Lambda}$, accounting for the fact that the matter density in the exterior is lower but the vacuum energy density is the same everywhere. Therefore,

[TABLE]

The last thing we have to do is to scale all parameters back to today’s values in the flat exterior spacetime. Or, more precisely, we need to tell the code for how long (in terms of exterior coordinate time) the simulation needs to run until the observer reaches $\tilde{z}=0$. This point turns out to be a bit subtle, because the exterior coordinate time is not synchronised with the proper time in the curved patch. Adding to the complication, the matter-synchronous clock runs differently in different places when viewed in the coordinates of the Poisson gauge.

The latter issue is fixed by choosing the location of the observer within the patch. Consider first an observer at the center (indicated as observer A in Fig. 1). This observer remains at rest in both coordinate systems used here, and the relation between matter-synchronous time (eigentime) and coordinate time is given by the lapse function evaluated at $r=0$. Specifically, the matter-synchronous clock runs slower, with $dt^{\text{(syn)}}=e^{\psi}dt^{\text{(P)}}$. This can also be seen from the generators $T_{1}$, $T_{2}$ of the time coordinate transformation.

The situation is illustrated in Fig. 2, showing the constant-time slices of the Poisson gauge in matter-synchronous coordinates. Choosing an initial redshift in the curved patch corresponds to setting the proper look-back time to the initial conditions. Since the constant-time slices are slightly tilted between the two coordinate frames, the simulation is initialised on a slightly different hypersurface. We choose the clocks to coincide in the exterior flat FLRW where the two coordinate systems match identically. This means that at the center of the curved patch, the proper look-back time to the initial hypersurface of the simulation needs to be corrected by a small amount that is given by the time coordinate transformation we have already computed. Finally, the simulation time is computed by taking into account the integrated effect of time dilation,

[TABLE]

In this equation, we solve for the unknown upper bound $t_{0}^{\text{(P)}}$ of the integral, which can be done with a simple root finding algorithm. One might also try to use the time coordinate transformation at the final time, but this is less robust for cases with large curvature. The expressions we derived for the coordinate transformation are valid in a small region around the initial time.

To evaluate the lapse perturbation $\psi$ in Eq. (2.23), we use the expressions derived in App. A. In the presence of a cosmological constant, we also account for the modified linear growth of the first-order solution in the exterior flat $\Lambda$ cold dark matter (CDM) model. If the exact behaviour of $\psi$ is difficult to estimate, one can also simply run a small low-resolution simulation and evaluate the integral on the right-hand side of Eq. (2.23) on the fly.

For an observer close to the perimeter of the curved patch, such as observer B in Fig. 1, the difference between the coordinate time in the Poisson gauge and eigentime becomes very small. In a closed Universe, it therefore reaches redshift $\tilde{z}=0$ sooner (in terms of simulation time) than the observer at the center.

Once the timeline of the simulation is fixed, we can infer relevant simulation parameters such as $\Omega_{\mathrm{m}}^{0}$, $\Omega_{\Lambda}^{0}$, $H_{0}$, $z_{\mathrm{in}}$, and so on, by evolving the exterior flat FLRW cosmology from the initial time.

2.4 Particle initial data and relativistic mass defect

In $N$-body simulations, the matter density is represented through a particle ensemble. For performance considerations, one typically chooses identical masses for all particles. To set initial conditions, the most common procedure is to lay down a homogeneous particle distribution, e.g. a regular crystal-like arrangement. Then the particle positions are perturbed to obtain the desired initial density field; see [24] for a review or [25] for the specific method we use which is suitable for a relativistic setting.

For the purpose of this discussion, let us work in the continuum limit, taking the number of particles to infinity. Following [23], we can obtain the particle displacement order by order in perturbation theory by demanding that the Hamiltonian constraint is satisfied on the initial hypersurface. For example, at first order we find

[TABLE]

where we use that $\psi_{1}=\phi_{1}$ and $\partial_{t}\phi_{1}=0$ in matter domination. Here, $\delta r_{1}$ is the displacement from the initial homogeneous particle distribution at first order. The solution for $\delta r_{1}$ in this case is

[TABLE]

The boundary condition was chosen $\delta r_{1}\rightarrow 0$ as $r\rightarrow 0$, required by regularity. A peculiar feature of this solution is that it does not vanish outside of $r_{2}$ due to the second term. This term implies that some additional mass is transferred into the curved region from infinity without perturbing the density in the flat region. This can be traced back to an effect called ‘relativistic mass defect’, see [21] for a discussion. In short, the active gravitational mass $M(r)$ is not simply the sum of the masses of all particles inside the radius $r$ because gravitational binding energy contributes to $M(r)$. Another way of looking at it is that the initial homogeneous template was prepared for a perfectly flat spatial geometry and, therefore, does not take into account the volume excess inside the curved patch: for a closed-universe simulation, the total proper volume of the initial hypersurface is larger than the proper volume of a flat FLRW model with the same boundary dimensions. If one computes the proper density before any displacement, one finds that it is reduced inside the curved patch compared to the exterior because of the extra volume.

As we need to ensure that the displacement field $\delta r_{1}$ has a solution compatible with periodic boundary conditions, we add the missing mass to the mass of the particles. Perturbing the particle mass inside the curved region by a constant $\delta m_{1}/m$, the Hamiltonian constraint at first order becomes

[TABLE]

leading to the solution

[TABLE]

The boundary condition $\delta r_{1}\rightarrow 0$ for $r\rightarrow r_{2}$ then implies

[TABLE]

At second order, the Hamiltonian constraint gives

[TABLE]

where $v_{1}$ is the first-order velocity given by $v_{1}=-\frac{2}{3}H_{\mathrm{in}}^{-1}a_{\mathrm{in}}^{-1}\partial_{r}\phi_{1}$. This equation can be solved for $\delta m_{2}/m$ in the same way as before.

In practice, instead of adding the extra mass only to the particles inside the curved domain, we evenly distribute it on all particles in the simulation. This implies that when initial conditions are set, indeed some mass needs to be transferred from the exterior to the curved patch, but we have ensured that the displacement field has a solution compatible with boundary conditions. We do not need to derive this solution explicitly, as it will be solved numerically by the iterative method we use to set the initial particle positions.

2.5 Matter perturbations in the curved patch

Including matter perturbations in the curved patch can be separated into two problems. First, one has to give a prescription for how to generate initial perturbations in a curved FLRW Universe. For convenience, we will assume that these can be provided in the Poisson gauge on a curved FLRW background. The curved background solution itself is the same in all gauges, so our coordinate transformations (2.12) and (2.13) can be used to transform these perturbations into Poisson gauge on a flat FLRW background. We already have everything in our hands for this second step, but the first one is a bit challenging in detail.

As noted in the Introduction, there are several interesting aspects that we want to briefly mention. In a curved Universe one cannot use standard inflation to predict the primordial power spectrum of perturbations, as inflation would always lead to a flat geometry. To the degree to which curvature survives inflation, the initial power spectrum will be sensitive to pre-inflationary physics. We want to emphasise again that this happens at exactly the scale we are interested in: the curvature scale of the Universe also defines the scale at which the power spectrum is modified. Interesting scenarios in which the primordial power spectrum can still be computed self-consistently include open inflation [26] or the no-boundary state [27] originally proposed by Hartle & Hawking [28]. Having observational evidence for any of these scenarios would be extremely exciting.

Another point to note is that the spectrum of perturbations is usually analysed in terms of eigenfunctions of the Laplace–Beltrami operator on the given background. For a flat background, these would be Fourier modes in a Cartesian basis or the corresponding spherical waves when using spherical polar coordinates. The mode functions in a closed Universe are hyperspherical harmonics, i.e. a generalisation of spherical harmonics to the three-sphere, which have a discrete spectrum. Mode functions on the three-hyperboloid can also be constructed and, like in the flat case, they have a continuous spectrum. We refer the interested reader to [29] where the implementation of such mode functions in a modern Einstein–Boltzmann solver is discussed in detail. What this means is that the initial perturbations are not easily expressed in terms of Fourier modes in our simulation box. An approximate correspondence can be established on short scales, but on scales where curvature is important, we should expect some corrections.

At present we have not yet developed an efficient method of initialising perturbations on the three-sphere and translating them to the coordinate system of our simulations. The procedure is clear at a conceptual level, but we leave the implementation to future work. In this work, we instead use an approximate method to initialise matter perturbations, ignoring the different nature of mode functions in different geometries. We therefore assume that the mode functions on the scales of interest are well approximated by the Fourier modes of the flat exterior geometry. The only thing we then need to do is to convert comoving wavenumbers between the two different coordinate systems.

In practice, our method of setting initial conditions consists of the following steps. First, we compute a linear prediction for the transfer functions of density, velocity, and gravitational potential for the curved model. This is done with the Einstein–Boltzmann code CLASS [29, 30]. Then we convert the wavenumber units from the CLASS output, which are for the curved FLRW model, into the units used in our flat exterior geometry. Unit conversion takes into account the difference in the expansion factor (relevant for comoving units) and the Hubble parameter (relevant because wavenumbers are given in units of $h\,\mathrm{Mpc}^{-1}$). In essence, this rescaling follows the philosophy of the separate universe approach. The resulting transfer functions in the flat coordinates are then used to generate a displacement field and a velocity potential for the particles, as well as the corresponding gravitational potential $\phi$ of the matter perturbations: for each Fourier mode, an initial random value of the primordial curvature perturbation is drawn from a Gaussian whose variance is given by the inflationary power spectrum, and the random value is multiplied by each transfer function to obtain the Fourier modes of the respective perturbation fields. The inflationary power spectrum follows a nearly scale invariant power law parameterised, as usual, with an amplitude $A_{\mathrm{s}}$ and spectral index $n_{\mathrm{s}}$, and we are mindful about unit conversion for the pivot scale at which $A_{\mathrm{s}}$ is defined. We refer the interested reader to Section 6 of [24] for the general picture and to Appendix A of [31] for details of our implementation. Since we do not want any matter perturbations outside of the curved domain, the final step of our procedure consists of apodising the displacement field, velocity potential, and gravitational field perturbations. This means that we smoothly trim these fields to zero as the distance from the center of the curved domain approaches the inner matching radius $r_{1}$ from the inside. Given our initial conditions for the simulation without matter perturbations, the perturbations are then linearly added.

3 Numerical results

We implement our methodology in the relativistic particle-mesh $N$-body code gevolution [31, 32]. Using a weak-field expansion in the Poisson gauge, this code includes relativistic corrections up to second order in the gravitational potentials $\phi$ and $\psi$, as well as frame dragging to its leading order. It is therefore ideally suited for our purpose. To verify our methodology, we performed three numerical experiments. The first two are designed to validate the approach at the background level. We run an Einstein–de Sitter model with large curvature parameter, $\tilde{\Omega}_{K}=-0.25$, to see if our implementation is robust for a wide range of parameters. We also run a $\Lambda$CDM model with moderate curvature, $\tilde{\Omega}_{K}=-0.1$, to demonstrate that a cosmological constant can be included and treated correctly. For both simulations, we measure the observed distance–redshift relations for different observers inside the curved domain. In the third experiment, we add matter perturbations in the curved $\Lambda$CDM model and measure observed angular clustering statistics in two redshift bins, again for different observers. This experiment shows that structure formation is consistent with the curved model and that the curved patch satisfies the cosmological principle, that is, that all observers see the same clustering statistics independent of location and direction.

3.1 Einstein–de Sitter simulation

For this test, we run a simulation for an Einstein–de Sitter model with $\tilde{\Omega}_{K}=-0.25$, $\tilde{\Omega}_{\mathrm{m}}=1.25$, $\tilde{h}=0.5$ and no matter perturbations in the curved patch since we want to test the performance of our approach at the level of a curved background cosmology. We use a simulation box of $6000\,h^{-1}\,\mathrm{Mpc}$ that embeds the curved patch inside a region with $r_{2}=2400\,h^{-1}\,\mathrm{Mpc}$ — note that distance units use the exterior value of $h\approx 0.512$ in our numerical setup. We start the simulation at $\tilde{z}_{\mathrm{in}}=25$, which corresponds to $z_{\mathrm{in}}\approx 25.464$ for our choice of parameters. The simulation has a resolution of $2048^{3}$ grid points for the mesh and uses $4096^{3}$ particles.

Using the existing code infrastructure of gevolution to build light cones on the fly, we generate data for the past light cones of two observers. The first one, observer A in Fig. 1, is located at the center of the curved patch and remains at rest with respect to exterior comoving coordinates. The second one, observer B in Fig. 1, is at a distance of about $2250\,h^{-1}\mathrm{Mpc}$ from the center and therefore close to the perimeter. The observed field of view points into the curved patch as illustrated in Fig. 1. As both observers are assumed to be comoving in the curved FLRW background, this observer has a nonvanishing peculiar velocity in the comoving frame of the simulation that, as we remind the reader, is following the expansion of the flat exterior. The peculiar velocity vector of the second observer points straight towards the center of the curved domain, following the flow of matter. We also note that, as discussed earlier, the two observers reach the hypersurface of $\tilde{z}=0$ at different simulation times (see Fig. 2). This is taken into account when we construct our light cones. For reference, the observation of the central observer is taken at $z=0$ according to our construction, but the second observer takes the observation at $z\approx 0.01524$. These redshift values correspond to the expansion factor of the simulation box, that is, to some coordinate time of the simulation, not to any observed values of redshifts in the data, $z_{\mathrm{obs}}$, to which we turn next.

The light cones are post-processed with a ray tracer that solves the geodesic equation and Sachs equations non-perturbatively, following the light rays backward in time from the observation event to the sources, which we take here to be a random sample of simulation particles. We find the unique geodesic to each source using a shooting method; see [33] for more details on our relativistic ray tracer. In this way, we obtain observed positions, redshifts, and angular distances to our sample of matter tracers. These are initially computed in the comoving frame of the simulation. For the second observer, we therefore need to apply a special relativistic boost to obtain the observables in the comoving frame of the matter. We do this by measuring the peculiar velocity at the observer position (here, we use the particle closest to the observer) and applying the appropriate boost transformations to observed redshift (relativistic Doppler boost), position, and angular distance (both are affected by relativistic aberration).

The final result is independent of any coordinate system used in the computation. In Fig. 3, we compare the observed angular distance as a function of the observed redshift with the exact relation predicted by the Einstein–de Sitter model. The left panel shows results for the observer at the center of the patch, where we find agreement at the level of one part in a thousand up to observed redshift $z_{\mathrm{obs}}\approx 1.7$, beyond which the past light cone crosses over into the flat embedding domain. This can be seen from the comoving distance in the curved patch, which we indicate in units of $r_{2}$ on the top horizontal axis. For reference, we also plot the exact distance–redshift relations predicted in the flat exterior model (red dashed curves) and the separate universe approximation (blue dashed curves). The latter consists of using the exact expansion history of the curved model but applying it to spatially flat comoving sections. This approximation leads to a $\sim 1\%$ error at redshift $z_{\mathrm{obs}}\simeq 1$ at our value of the curvature parameter, originating from the purely geometrical effect of curvature that is independent of the expansion.

The right panel of Fig. 3 shows corresponding results for the second observer, who is located close to the perimeter of the curved domain. Due to this vantage point, the lines of sight for a certain fraction of the sky can extend to much higher redshifts before the past light cone crosses over into the flat embedding domain (see Fig. 1 for an illustration). The reduced symmetry of this configuration also means that not all lines of sight take similar paths through the curved patch. To show that the distance–redshift relation remains isotropic from the point of view of the observer as long as the observed sources are contained within the curved patch, we separately plot samples that lie close to the central ray, pointing directly inward, and samples that lie close to $30$ degrees away from that axis. Both samples give measurements consistent with the curved model to one part in a thousand, this time up to redshifts beyond $z_{\mathrm{obs}}=8$. At such high redshifts, the separate universe approximation would be off by several percent.

3.2 Curved $\Lambda$CDM simulation without matter perturbations

For this test, we choose somewhat less extreme parameters that are closer to what is required for the forward modelling of real observations. Our curvature parameter, $\tilde{\Omega}_{K}=-0.1$, is still ruled out by observational constraints, but our aim here is to illustrate the effect of curvature. Importantly, however, this simulation includes a cosmological constant. The cosmological parameters are $\tilde{\Omega}_{K}=-0.1$, $\tilde{\Omega}_{\mathrm{m}}=0.4$, $\tilde{\Omega}_{\Lambda}\simeq 0.7$, $\tilde{h}=0.7$, and the remaining energy density is in a small amount of residual radiation (photons and massless neutrinos) which we keep for later convenience but which have no noticeable effect on the background evolution. We use a simulation box of $4500\,h^{-1}\,\mathrm{Mpc}$ that embeds a curved patch inside a region with $r_{2}=1800\,h^{-1}\,\mathrm{Mpc}$. Starting at redshift $\tilde{z}_{\mathrm{in}}=15$, we find that the exterior model has $z_{\mathrm{in}}\approx 15.719$ and $h\approx 0.716$. The simulation has a resolution of $1536^{3}$ grid points and uses $3072^{3}$ particles. Note that although these numbers are smaller than for the Einstein–de Sitter run, the physical resolution is significantly higher due to the choice of other parameters.

We place observer A at the center of the curved patch as before and observer B at a distance of about $1650\,h^{-1}\,\mathrm{Mpc}$ from the center. The coordinate time at which the latter makes the observation corresponds to $z\approx 0.00376$ in the exterior model. Apart from these details, the construction of the observables proceeds along the same steps as outlined for the Einstein–de Sitter case. Fig. 4 shows the results compared to the exact distance–redshift relation predicted in the curved $\Lambda$CDM input cosmology. We find excellent agreement, better than one part in a thousand, for both observers and for the entire range of redshifts for which the respective past light cones are contained within the curved domain. For observer A (left panel of Fig. 4), the usable light cone extends to about $z_{\mathrm{obs}}=0.6$, whereas observer B (right panel of Fig. 4) has access to much higher redshifts. As before, we find that the Hubble diagram remains isotropic within the curved domain. Across all samples, the small residual is a nearly constant offset of about $0.05\%$, which would correspond to an error in the fiducial Hubble parameter of the same order: changing the Hubble parameter of the reference model shifts a curve up or down in these plots without affecting the shape. It is also interesting to note that, despite the smaller curvature parameter of $\tilde{\Omega}_{K}=-0.1$, the separate universe approximation (blue dashed curves) still exceeds an error of $1\%$ beyond $z_{\mathrm{obs}}=1$.

3.3 Curved $\Lambda$CDM simulation with matter perturbations

In our final numerical test, we run the curved $\Lambda$CDM simulation another time, adding matter perturbations in the initial conditions. These are constructed as described in Sec. 2.5, with a primordial scalar amplitude $A_{\mathrm{s}}=1.8113\times 10^{-9}$ at the pivot scale $\tilde{k}_{\mathrm{pivot}}=0.05\,\mathrm{Mpc}^{-1}$, where the comoving wavenumber is defined with respect to the curved background, and scalar spectral index $n_{\mathrm{s}}=0.9619$. Since the Einstein–Boltzmann solver requires sensible values for the baryon and photon densities, we split $\tilde{\Omega}_{\mathrm{m}}$ into a baryonic contribution and a CDM contribution, with $\tilde{\Omega}_{\mathrm{b}}=0.05$ and $\tilde{\Omega}_{\mathrm{cdm}}=\tilde{\Omega}_{\mathrm{m}}-\tilde{\Omega}_{\mathrm{b}}=0.35$, respectively. The photon density is fixed by setting the temperature of the cosmic microwave background to $2.7255\,\mathrm{K}$. For completeness, we also include three standard massless neutrino species. The density of these radiation components has a negligible effect on the evolution in the redshift range covered by the simulation, but their presence is of course important for predicting the linear initial conditions of the matter perturbations.

The aim of this test is to measure the observed clustering statistics of matter in the curved patch. In the spirit of our previous test, we want to construct a pure observable that is independent of the coordinate system of the simulation. A suitable observable is given by the angular two-point function of matter within a redshift bin. Its harmonic transform, the angular power spectrum $C_{\ell}$, can be predicted in linear theory using CLASS. To estimate the angular power spectrum from simulations, we first construct light cones as before, but using a dense sampling of particles. One in every sixteen particles is used as a tracer, arranged in a face-centered cubic packing embedded in the initial particle template. This regular pattern avoids the introduction of shot noise. Note that we keep the observers comoving with the curved background, that is, we use their respective peculiar velocities from the simulation without matter perturbations.

Our observer B uses a cone-shaped field of view covering about 4500 square degrees, so in order to gather more statistics, we place four equivalent copies of this observer inside the curved domain, located at the vertices of an inscribed regular tetrahedron. The idea is that although their fields of view overlap, the data within some thin redshift bins are largely independent and can be treated as separate realisations.

Once the observed positions and redshifts of the matter tracers have been computed, we select two thin redshift bins, centered on $z_{\mathrm{obs}}=0.5$ and $z_{\mathrm{obs}}=1.0$, using a top-hat selection window of width $\Delta z_{\mathrm{obs}}=0.05$. Note that the curvature parameter reaches its maximum value of about $12.33\%$ around $\tilde{z}=0.5$, and starts decaying in the later Universe due to the presence of the cosmological constant. We create pixel maps of the particle distribution in these redshift bins using the HEALPix package [34], and compute their pseudo-$C_{\ell}$ estimator with NaMaster111[

 https://github.com/LSSTDESC/NaMaster](https://github.com/LSSTDESC/NaMaster) [35]. The redshift bin at $z_{\mathrm{obs}}=0.5$ is also available on the past light cone of observer A at the center of the curved domain. We want to check if the clustering statistics are independent of observer position, meaning that the cosmological principle is satisfied. We therefore extract four separate fields of view from the central observer’s sky, matching the cone-shaped selection of observer B, and apply the same pseudo-$C_{\ell}$ analysis. To minimise their correlations, the fields of view are again arranged at tetrahedral angles.

The left panel of Fig. 5 shows the results of our analysis. The four pseudo-$C_{\ell}$ estimates for each observer type and each redshift bin are averaged to reduce the sample variance. We confirm that the observed clustering amplitude is consistent with the linear prediction from CLASS (dashed curves) everywhere, except for high multipoles (small scales) at low redshift where we expect an enhancement due to nonlinear structure. We also find that the angular spectra measured by observer A and observer B are statistically indistinguishable, in line with the cosmological principle. For illustration, the right panel of Fig. 5 shows the matter density on a cross-section of the simulation volume at $z=0$.

Although possible in principle, we did not attempt a comparison to the separate universe approximation for this observable, as this may depend on the way the light cone is constructed. Treating the separate universe model as a flat geometry with a modified expansion history that matches the one of a curved model, one can think of curvature as a ‘dark fluid’ with equation of state $w=-^{1}\!/\!_{3}$. The linear prediction in this scenario has been discussed in [36], where it is directly compared to the linear prediction with geometric curvature. For $\Omega_{K}=0.01$, their study indicates a $\sim 1\%$ error on the angular power spectra of matter over a large range of multipoles and redshifts if the geometric effect of curvature is neglected; see their Fig. 2. If this result extrapolates to larger values of curvature, and specific details of their analysis such as redshift binning, linear bias and magnification bias do not change the picture, we might expect $\sim 10\%$ errors on the angular power spectra of matter for the separate universe approximation with our value of $\tilde{\Omega}_{K}=-0.1$. It is also worth pointing out that the geometric effects of curvature, at least for the examples studied in [36], turn out to be somewhat larger than the purely dynamical effects that are related to the modification of the expansion history and captured by the separate universe approximation.

4 Conclusion

In this paper, we have developed a methodology for running cosmological $N$-body simulations of structure formation in the presence of global spatial curvature. Our relativistic approach takes non-Euclidean distance and volume effects fully into account for the first time. These are mostly relevant for the construction of observables on the light cone that span cosmological distances. By embedding a curved patch of finite size into a periodic embedding space, we can work within the technical constraints of existing $N$-body simulation codes, which commonly rely on periodic boundary conditions. Starting from an exact LTB solution in the synchronous gauge, we present the gauge transformation to Poisson gauge (the gauge used in our relativistic simulation code) to second order, ensuring high accuracy even for large values of curvature. We demonstrate that our approach reproduces some key predictions for the curved models with high accuracy: in the matter-comoving frame, the distance–redshift relation is isotropic around any point in the curved patch and agrees with the one of the curved input model; the angular two-point correlation function of matter in observed reshift bins is also consistent with expectations and satisfies the cosmological principle.

Our work advances the state of the art in forward-modelling observables of models with curvature in several important ways. While the separate universe approximation gives accurate predictions for the clustering statistics on small scales, it lacks a rigorous prescription for constructing large contiguous survey footprints. We have shown that the distance–redshift relation becomes inaccurate in this approximation when taken over large enough distances. This is in line with previous findings that modifying the expansion history alone does not capture the full impact of curvature for observables on our past light cone. In some cases, for example cosmic shear, the effect has been shown to be below sample variance of current surveys for realistic values of curvature [37]. Angular clustering statistics, on the other hand, can be more strongly affected [36] and can therefore provide an interesting probe of curvature.

The LTB model has been used previously in simulations of structure formation; see [38, 39] for two relevant examples. Although these studies focus on modelling large spherical inhomogeneities in a standard flat cosmology, it is clear that their methodology could be adapted to our use case. This suggests that Newtonian $N$-body codes should be able to reproduce our results to a very good approximation; see [39] for a discussion. However, by not taking the Newtonian limit, our approach is more rigorous, and we clarify several subtleties. For example, we take into account the relativistic mass defect, which is almost $1\%$ in our curved $\Lambda$CDM simulation with $\tilde{\Omega}_{K}=-0.1$. The size of this effect, like most of the corrections we consider, scales with the difference between Euclidean and curved volume and therefore depends not only on $\tilde{\Omega}_{K}$ but also on the survey footprint we want to model. The volumes needed to forward-model future surveys may exceed the scale of previous LTB simulations by orders of magnitude, degrading the accuracy of the Newtonian approximation.

While our current implementation only allows for the simulation of a closed universe, the methodology can be generalised to the open case. The minor complication is the need for a compensating mass shell to obtain a plausible sequence of matching radii for the exact LTB solution. We plan to implement this in a future iteration of our code. Furthermore, setting the initial conditions for matter perturbations on scales approaching the curvature scale would require more careful consideration. Our present approach essentially follows the separate universe approximation in assuming that the mode functions on the curved geometry are given by the standard Fourier series instead of hyperspherical harmonics. However, if we were given some initial data on the three-sphere, we would already have computed the coordinate transformation for translating them into our simulation framework. We plan to revisit this aspect in future work.

Acknowledgments

We thank Selim Oueslati for discussions. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID uzh34. We acknowledge financial support from the Swiss National Science Foundation.

Code availability

Our implementation is available in the new version 1.3 of gevolution, which can be found at https://github.com/gevolution-code/gevolution-1.3.

Appendix A Second-order gauge transformation

Here we present the gauge transformation from the LTB metric (2.3) in comoving synchronous gauge to Poisson gauge, where it takes the form given in Eq. (2.1). Let us expand the metric functions in both gauges around flat FLRW,

[TABLE]

Using the expressions for $E(r)$ and $R(t,r)$ from Eqs. (2.9) and (2.10), we find

[TABLE]

where we defined the shorthand

[TABLE]

Since the line element $ds^{2}$ is invariant, we have

[TABLE]

where both metrics have to be evaluated at the same space-time event (which has different coordinate values in the two coordinate systems). The first order of Eq. (A.10) gives the following relations,

[TABLE]

Combining the latter two, we can solve for $L_{1}$,

[TABLE]

Inserting this solution into Eq. (A.12) fixes $T_{1}$,

[TABLE]

where we use the fact that $H^{2}a^{3}$ is independent of time. We can finally read off the solutions for $\phi_{1}$ and $\psi_{1}$,

[TABLE]

The second order of Eq. (A.10) gives the following relations,

[TABLE]

Subtracting the two last equations gives

[TABLE]

Using the first-order solutions and after some algebra, this simplifies to

[TABLE]

and hence

[TABLE]

We may now compute $a^{2}\partial_{t}L_{2}$ that appears in Eq. (A.19). Using the definition of $h(r)$ together with the fact that $H^{2}a^{3}$ is independent of time, we can derive

[TABLE]

We therefore find

[TABLE]

and hence

[TABLE]

where we use integration by parts on some of the terms to simplify the expression.

Next, we compute $\partial_{t}T_{2}$, as it enters Eq. (A.18),

[TABLE]

and therefore

[TABLE]

Lastly, we solve for $\phi_{2}$, e.g. using Eq. (A.21),

[TABLE]

Appendix B An exact solution in the Schwarzschild region

Outside of the curved patch, i.e. in the vacuum region and the flat FLRW embedding spacetime, the model corresponds exactly to the one of Einstein & Straus [20]. It describes a Schwarzschild spacetime embedded in an expanding FLRW solution, with appropriate junction conditions satisfied at the interface. Here we present the metric of this entire manifold, valid on both sides of the junction hypersurface, $\left\{r=r_{2}\right\}$, written in terms of FLRW comoving coordinates and using the Poisson gauge. Specifically, our metric takes the form given in Eq. (2.1), with $a(t)\propto t^{2/3}$ the FLRW scale factor and $\phi=\psi=0$ in the FLRW region. Our aim here is to obtain exact expressions for $\phi$ and $\psi$ in the vacuum (Schwarzschild) region.

Instead of attempting to derive the gauge transformation from the LTB metric in a non-perturbative way, we may instead consider solving Einstein’s equations directly, given that the right-hand side vanishes in vacuum. Due to symmetries, these reduce to only three coupled partial differential equations,

[TABLE]

The momentum constraint (B.2) implies

[TABLE]

which we may integrate on a constant-time hypersurface to get

[TABLE]

The constant of integration is fixed by the requirement that both $\psi$ and $\partial_{t}\phi$ must vanish identically at the junction. We can insert this solution into the Hamiltonian constraint (B.1) to get a nonlinear ordinary differential equation (ODE) for $\phi$ on a constant-time hypersurface,

[TABLE]

To get further analytical insight, we make the following change of variables,

[TABLE]

which casts the ODE into the convenient form

[TABLE]

where a prime denotes the derivative with respect to the argument, $x$ in this case. Multiplying with $u^{\prime}$ yields

[TABLE]

which integrates to

[TABLE]

We can fix the constant of integration by noting that at $x=0$ ($r=r_{2}$) we have $u(0)=\sqrt{3}aHr_{2}\equiv u_{0}$ and $u^{\prime}(0)=-u_{0}$ because both $\phi$ and $\partial_{r}\phi$ must vanish. We find $C=-\frac{1}{3}u_{0}^{3}$, and therefore

[TABLE]

According to a result by Biermann & Weierstrass (see Appendix A of [40] for details, including original references), $u$ can be expressed in terms of the Weierstrass elliptic function $\wp(-x;g_{2},g_{3})$,

[TABLE]

with the Weierstrass invariants $g_{2}=\frac{1}{12}$ and $g_{3}=-\frac{1}{432}(2+u_{0}^{6})$.

After changing variables back from $x$ to $r$, an exact expression for $\phi$ is given by

[TABLE]

where we note that the time dependence is encoded entirely in $u_{0}$, which also enters the definition of $g_{3}$. We remind the reader that the external FLRW solution is given by $a\propto t^{2/3}$ and hence $H=\frac{2}{3}t^{-1}$. An expression for $\psi$ can now be derived from Eq. (B.5). The solutions for $\phi$ and $\psi$ are illustrated in Fig. 6.

To our knowledge, this is the first time that this exact form of the Schwarzschild metric has been reported.