Can gravitational waves halt the expansion of the universe?

We numerically investigate the propagation of plane gravitational waves in the form of an initial boundary value problem with de Sitter initial data. The full non-linear Einstein equations with positive cosmological constant $\lambda$ are written in the Friedrich-Nagy gauge which yields a wellposed system. The propagation of a single wave and the collision of two with colinear polarization are studied and contrasted with their Minkowskian analogues. Unlike with $\lambda=0$, critical behaviours are found with $\lambda>0$ and are based on the relationship between the wave profile and $\lambda$. We find that choosing boundary data close to one of these bifurcations results in a"false"steady state which violates the constraints. Simulations containing (approximate) impulsive wave profiles are run and general features are discussed. Analytic results of Woodard and Tsamis, which describe how gravitational waves could affect an expansion rate at an initial instance of time, are explored and generalized to the entire space-time. Finally we put forward boundary conditions that, at least locally, slow down the expansion considerably for a time.


Introduction
Space-times containing plane gravitational waves have seen extensive analytical study over the years and many closed form solutions, which necessarily assume certain symmetries or wave profiles, now exist and their properties are known (see [3] for an excellent overview.) While there are a number of analytic solutions for the propagation and collision of waves assuming a vanishing cosmological constant [4,5,6,7,8,9], the non-vanishing cosmological constant analogues pale in number, and there are no closed form solutions for colliding waves in this case.
Penrose's cut-and-paste method [8,10], which cuts Minkowski space-time along a null hyperplane, shunts one half along the same surface and then pastes the two halves back together gives rise to a space-time with one impulsive gravitational wave (i.e. with a Dirac delta function wave profile.) This has been generalized to non-zero, constant curvature backgrounds [11,12,13,14,15,16] where the wave fronts are topologically spherical for λ > 0 and hyperboloidal for λ < 0. There do not exist however, closed form solutions to the full non-linear Einstein equations with λ = 0 that contain gravitational waves with plane symmetric wave fronts.
De Sitter space-time, the unique solution to the Einstein vacuum equations with constant positive scalar curvature, can be thought of as a model of a universe which is expanding at an accelerated rate from the positive λ contribution. Quantum gravitational back-reaction on inflation [17] allows for the creation of cosmic scale gravitational radiation, where, if one does not account for their creation, can be modelled completely classically through gravitational perturbations of de Sitter space-time [1,2]. It is theorized that such a background of radiation may weaken the expansion and even halt it completely. Analytical calculations have been done to explore this hypothesis by studying how an expansion parameter and its time derivative could be manipulated through such a field at an initial instance of time. The question of what happens away from this surface remains unanswered, and attempting to answer this in the full nonlinear regime analytically would be very difficult, if not impossible.
In this paper, we numerically evolve the Einstein vacuum equations with positive cosmological constant in plane symmetry with the goal of shedding light on the above topics. To do so, we implement an initial boundary value problem following Friedrich and Nagy [18], which is wellposed, and allows us to generate gravitational perturbations through boundary conditions rather than solving the constraints. This framework has already been implemented and numerically validated in previous work [19] for λ = 0. We generalize this to an arbitrary cosmological constant as well as the inclusion of matter terms through components of Φ ab = −(1/2)R ab + (1/8)Rg ab and scalar curvature Λ = (1/24)R for completeness. We follow the conventions of Penrose and Rindler [20,21] throughout.

Review of plane gravitational waves with λ = 0
Here we briefly present the space-times of a single impulsive gravitational plane wave and the collision of two, which are colinearly polarized, with λ = 0. This can be accomplished by summarizing the Khan-Penrose solution [7], which describes the latter.   1 showcases the Khan-Penrose solution in null coordinates u, v, where the two spatial dimensions that span the planes are suppressed, so that each point represents a plane. Null curves are represented by lines with slope ±1 and the impulsive waves are given by Ψ 0 = δ(v), Ψ 4 = δ(u), where δ is the Dirac delta function, so their path is given by the dashed lines. These lines split the space-time into four regions. The lower region is Minkowski space-time, the two side regions are space-times containing one propagating wave only, and the top region is the interaction region after scattering. All four regions can be represented by the single line element where p := u Θ(u) and q := v Θ(v). The interaction region contains a spacelike curvature singularity on the surface u 2 + v 2 = 1 and can be seen as such due to the divergence of, for example, the Weyl invariant I. The region containing only the Ψ 0 wave is where u < 0 and v ≥ 0 and the line element Eq. (1) reduces to This region and its Ψ 4 counterpart contain a fold singularity along v = 1 resp. u = 1. As Eq. (2) can be transformed to Minkowski space-time by a coordinate transformation, one would think this is merely a coordinate singularity. However looking closer one sees that this is not the case as there does not exist a C 1 extension from this region to v = 1 resp. u = 1 [22]. Further, it is found that a certain projection of the u = constant, v = constant surfaces into Minkowski space in standard null coordinates converge at v = 1. This has the consequence, which is discussed in more detail in Sec. 5.1, that the spin-coefficients ρ and ρ , which when positive, represent the converging of a null geodesic congruence along l a and n a respectively, diverge to positive infinity, showing an ever strengthening contraction of null rays in both null directions.

General setup
We write the Einstein equations in the form of an IBVP following Friedrich and Nagy [18] with the additional imposition of a pair of commuting space-like Killing vectors that represent our plane symmetry. Further, we include matter coming from an energy momentum tensor T ab . A detailed explanation of this process in vacuum with vanishing cosmological constant has been laid out in [19]. We only give a brief summary here, emphasising the differences when including a non-vanishing cosmological constant and matter. The Einstein equations take the form where and Λ and Φ ab correspond to the trace and tracefree part of the Ricci tensor Φ ab and λ is the cosmological constant.
To start setting up our gauge, we first assume our space-time can be foliated by planes. We then define the coordinates t, z for time and the direction of wave propagation respectively, both being constant within the planes. Using the holonomic basis we define the null tetrad where A, B, ξ, η are functions of (t, z) only. This leads to the metric To obtain equations for the metric functions and find algebraic relations for the spincoefficients (due to the plane symmetry assumption) we apply the commutator equations (see [20] Eq. (4.11.11)) to the coordinates. To obtain equations for the spin-coefficients we use the curvature equations (see [20] Eq. (4.11.12)). To obtain equations and algebraic relations for the components of the Weyl tensor C abcd , Φ ab and Λ, we use the equations coming from the Bianchi identity (see [20] Eqs (4.12.36-4.12.41)).
The algebraic conditions are found to be Following Friedrich and Nagy, we set where the free function F = χ + if is a freely specifiable gauge source function and µ is taken as a system variable. χ is the mean extrinsic curvature of the z = constant hypersurfaces and f determines the rotation of the m a frame vector along (∂ t ) a . The geometrical interpretation of the new variable µ can be explained in the gauge F = ρ − ρ, which is the gauge used for most of our results and turns out to be the Gauß gauge. Although predisposed to develop caustics, an expanding universe, which we consider here, acts to counter this. The fact we are in the Gauß gauge can be seen immediately by noticing that the only non-vanishing component of the acceleration of the unit time-like vector (∂ t ) a along itself is proportional to γ +γ + +¯ = F +F + 2(ρ − ρ ) = 0 (13) for this choice of F . The "acceleration" z a ∇ a z b of the space-like unit vector z a := A(∂ z ) a along itself is proportional to µ +μ, which gives an interpretation for the real part of µ. The imaginary part just corresponds to a phase change of m a . It is found that the equations for η, ξ decouple from the others, and as they are superfluous to the results subsequently presented we do not include them in the system. The evolution equations are while the constraints take the form To supplement the above, the divergence free condition on the energy-momentum tensor (equivalently the Bianchi identity, which are given in [20], see Eq. 4.12.40) gives Considering only the vacuum equations with cosmological constant, i.e. Φ ab = 0, Eqs (16a)-(16b) are identically satisfied and Eqs (14a)-(14i), Eqs (15a)-(15d) comprise a closed system of equations, where the evolution equations are symmetric hyperbolic and the constraints propagate. When matter terms are present and one takes into account Eqs (16a)-(16b), it is still found that the above system is symmetric hyperbolic and the constraints propagate. The resulting subsidiary system is In order to close the system, one must in general couple it to equations describing the evolution of matter. There is a lot of freedom in this choice and it depends very much on the physical situation one wants to model. In general this choice will alter the principal part and, as a consequence, symmetric hyperbolicity and constraint propagation could be lost. Two useful quantities are now introduced for monitoring the behaviour of the evolved space-time. First we note that the extrinsic curvature of our t = constant is the induced 3-metric on the surfaces and t a = (1 − B 2 ) −1/2 (dt) a is the unit conormal. We then define a local expansion parameter proportional to the mean extrinsic curvature K a a as which is used to monitor the expansion rate of the space-time along the time coordinate vector field. The Weyl scalar curvature invariants are useful tools for identifying whether a singularity is a curvature singularity. In the absence of matter and with our plane symmetry assumptions, the real part of C abcd C abcd is the Weyl scalar curvature invariant We define the wave profile where b = 35π/4 so that the area of the profile is a = π/b 0 p(x)dx and the amplitude is 32a. We take a as a measure of the strength of the wave. The boundary conditions for Ψ 0 and Ψ 4 will make use of p(x) and are chosen in the subsequent sections.

De-Sitter space-time
We investigate a variety of cases of plane gravitational waves propagating in de Sitter space-time (dS). The unperturbed metric in inflationary coordinates can be written [23] which covers half of the space-time and matches our setup for plane symmetry. This represents an expanding universe of the FLRW type. The appearance of A 0 := A(0, z) is used to scale the spatial directions and will be useful later. It is useful to write dS in terms of null coordinates as with transformations The Minkowskian analogue of the above can be found in the limit H → 0. In our formalism Eq. (21) gives the initial data with the remaining system variables, gauge quantities and matter terms vanishing. We will use the negative non-vanishing initial data, corresponding to a future expanding universe, and set Φ ab = 0.
We incorporate into the system null coordinates u(t, z), v(t, z) which satisfy l a ∇ a u = 0 and n a ∇ a v = 0 respectively. Their initial and boundary data are fixed by Eqs (24), (25) so that when no wave is present we reproduce the same null coordinates as in Eq. (21) when F is chosen appropriately. The above expressions for u, v were chosen so that initially u(0, −1) = 0 = v(0, 1). Having u, v available allows us to define the semi-invariant coordinates (T, Z) by T := √ 2(v + u) and Z := √ 2(v − u) with which we can produce Penrose-Carter diagrams, i.e. diagrams where null curves are lines with slope ±1. When in exact dS, as t → ∞ we obtain

Numerical setup
We utilize the Python package COFFEE [24], which contains all the necessary functionality to perform a numerical evolution using the method of lines. We discretize the z-direction into equi-distant points in the interval [−1, 1] and approximate the zderivative using Strand's finite difference stencil [25] which is fourth order in the interior, third order on the boundary and has the summation-by-parts property [26]. We march in time using the explicit fourth order Runge-Kutta scheme with a timestep determined by ∆t = c ∆z, where ∆z is the step size in the z-direction and c is the CFL constant. Unless otherwise stated we take c = 0.5. Boundary conditions are imposed using the Simultaneous Approximation Term (SAT) method [27] with τ = 1. This particular selection of numerical methods within COFFEE has proven to be numerically sound for a variety of different systems (see for example [19,28]). In the subsequent situations, all constraints are verified to converge at the expected order everywhere.

An analytical view
Before analyzing the numerical results, it is worthwhile to perform a small analytic study of the propagation of one wave when either Minkowski or de Sitter initial data are taken. Firstly, the evolution equations for ρ and σ (Eqs (14c) and (14e)), which have a close relationship to Sachs' optical equations, give with vanishing Φ ab For the case of Minkowski initial data, which is obtained by setting λ = 0 in the de Sitter initial data, and where we choose F (t, z) = 0 to extend the exact gauge of dS to the whole space-time, we find the following: The introduction of a non-zero Ψ 0 on the right boundary causes σ to become non-zero there. This in turn causes ρ to become non-zero. As ∂ t ρ > 0, we find that ρ will inevitably diverge. Further, one can see by looking at the evolution equations for the primed spin-coefficients, all primed spin-coefficients stay zero throughout the space-time, due to the forever zero Ψ 4 . Further, this implies that Ψ 2 = 0 everywhere and thus the Weyl invariant I given by Eq. (19) also remains zero everywhere. These are well known result for propagation of a single plane gravitational wave in Minkowski space-time, see [3] for an overview (in a different gauge). The case of expanding de Sitter initial data, with non-vanishing λ and again choosing F (t, z) = 0, is quite different. A non-zero Ψ 0 leads to a non-zero σ as before, but now a non-zero σ leads to a non-zero σ as well as a non-zero ρ. This non-zero σ then makes ρ and even Ψ 4 non-zero, implying the non-linear back-reaction effect is realized. This in turn leads to a non-zero Weyl invariant. The added complexity of the non-zero λ, which couples all system variables together in a complicated, non-linear way, stops us from concluding statements analogous to the Minkowski case as above, emphasising the need for numerics.

Numerical analysis
We now fix λ = 3 and choose the boundary conditions to be where p(v) has the area of the wave packet as a parameter, and the change of area is realized by a change in amplitude. We perform evolutions with wave areas a taking the values 1.67, 1.6765, 1.6769105, 1.6769106, 1.676912, 1.67695, 1.68 for reasons that will become apparent shortly. In all these cases, once the wave has entered and subsequently left the computational domain, the space-time is fully excited in that all system variables have evolved away from their original values. For the case of four smallest values of a we find that the space-time asymptotes back to the de Sitter space-time everywhere. This indicates that the wave has been wiped out by the accelerated expansion, already in stark contrast to the Minkowskian analogue where a future singularity is guaranteed. Fig. 2 shows a contour plot of Ψ 0 and H over the entire space-time. It is clear that H decreases due to the addition of the gravitational wave, but then settles back down to its original value of one. The only remaining effect after the wave has passed is the time delay between different regions of space-time, such as the left and right boundaries. To see how the representation of the null directions l a and n a in the coordinate basis change during the simulation, we look at the metric functions A and B. It is seen that A → 0 as in the exact de Sitter case, representing the exponential expansion, and although initially B increases to some value less than one, it asymptotes back to zero. Notably, the rate at which A → 0 and B → 0 causes the dtdz metric coefficient to asymptote to a constant non-zero value and the dz 2 metric coefficient to diverge to positive infinity. The fact that A and B never actually reach zero implies our gauge remains regular. Figure 3: Two diagrams showcasing how the vectors l a , n a , ∂ a t and ∂ a z behave on the right boundary when a future singularity occurs.
For the four largest values of a we find that the simulation crashes after some time due to A → 0 and B → 1 in finite time on the right boundary, the same as in the Minkowskian case. Fig. 3 shows how this affects the relevant frame vectors there, where we note the relationships where t a and z a are normalised. The left diagram is with respect to the {l a , n a } null basis defined in the tangent space and exemplifies the fact that t a = ∂ a t and is always normalised to one. It also showcases that the evolution of z a can cause trouble. This can be seen by noting that as B → 1 the z = constant surfaces become characteristic. The right diagram looks at another potential issue, this time in our (T, Z)-coordinates. In this case t a is no longer given by a vertical line, but l a and n a remain as lines with slope ±1 from the definition of u and v. The "shrinking" of the n a and the "growing" of l a is due to both coefficients of n a in the coordinate basis approaching zero, and enforces that t a is proportional to the sum of the two and that their normalisation conditions are maintained.
This behaviour affects the expansion rate H and we find it decreases and actually diverges to −∞ on the right boundary, as shown in Fig. 4. The features discussed above indicate that for these larger wave areas, the expansion rate of the space-time is not strong enough to overcome the contractivity of the wave, and a future singularity is formed. In the Minkowski case the analogue is a fold singularity as discussed in Sec. 5.1. In our de Sitter case, Fig. 4 and Fig. 5 show that the Weyl invariant I is diverging on the right boundary (and similarly close to the right boundary), unlike the Minkowski case, and adds emphasis to the classification of a curvature singularity. A final note is that changing the polarization of the wave, implemented by replacing p(x) with e iφ p(x) for some real constant φ, does not affect the expansion rate or Weyl invariant as seen in Fig. 2, Fig. 4 and Fig. 5.

Critical behaviour
An obvious question has been raised: What is the critical behaviour when the ingoing wave has the critical wave area a c that separates these two distinct futures? One can obtain a c by using a simple binary search. For λ = 3 this is found to be 1.67691055 < a c < 1.67691056. Fig. 6 shows ρ and σ along the right boundary for various wave areas close to a c . It is clear that as a → a c an interval appears where ρ and σ are constant in time and the interval becomes longer the closer a is to a c . This indicates that a special critical behaviour may exist. All system variables except A become constant in a finite t-interval on the boundary which becomes larger the closer to a c we take our wave area, and µ, ρ, ρ , σ, σ take on values different than their initial ones. This implies a steady state solution, different from the de Sitter space-time. It turns out we can solve for unconstrained steady state solutions (but with A a function of time) algebraically by setting all time derivatives except A to zero in our evolution system, as well as taking Ψ 0 = Ψ 4 = F = 0. One of these solutions is found to match the values we see numerically. However, this exact solution does not satisfy the constraint equations, and is thus a "false" steady state. This can be seen explicitly during our evolution, by noticing that the constraints do not converge and are wildy violated during this steady state period, see Fig. 7. This is a consequence of our free evolution scheme, which by definition, is "free" from enforcing the constraints to be satisfied. It is found that the only free steady state solution (with A varying in time) that also satisfies the constraints in the case of a positive cosmological constant with Ψ 0 = Ψ 4 = F = 0 is the de Sitter space-time. The fact that no critical behaviour exists for this wave profile ansatz will be important when attempting to find a solution where the expansion is halted with gravitational radiation, and will be discussed in detail in Sec. 7.

An impulsive wave
Many analytical solutions describing gravitational waves in the literature have an impulsive wave profile, i.e. Ψ 0 = δ(v) where v is a null coordinate, which is a consequence of the cut-and-paste method of Penrose [8]. An example is the propagation of a single impulsive gravitational plane wave with λ = 0 given by Eq. (2), where Ψ 2 = 0 = Ψ 4 . To date, an exact solution for a single propagating plane gravitational wave with λ > 0 has not been found. One cannot use Penrose's cut-and-paste method to find such a solution because this leads to wavefronts that are spherical or hyperboloidal when λ > 0 or λ < 0 respectively [16]. Thus, to try shed some light toward an analytic solution, we numerically evolve our system with λ > 0 and with one ingoing wave, whose wave profile approximates the Dirac delta function. We set Ψ 0 (v, 1) = q(v) where where b = 35πa/128 and q(x) has the property that lim a→∞ q(x) = δ(x). We also change our gauge and fix F by the condition that ∂ t B = 0. This matches the gauge of the exact solution given by Eq. (2) and yields F = ρ − ρ. We choose a = 128, 256, 512, 1024, populate our spatial interval z ∈ [−1, 1] with 6401 equi-distant points to accurately resolve these steep wave profiles and choose λ = 0.6 and λ = 1.2 to exemplify futures that do and do not have a singularity respectively.
For λ = 1.2, to see the effect of the limit a → ∞, Fig. 8 shows the Weyl components along z = 0. These seem to indicate that in this limit, they all vanish for v > 0 along z = 0. By inspection it is clear that this happens along any z = constant curve once the wave has past and thus in the whole region v > 0. Further, all system variables asymptote back to dS after the wave has past, and no singularity is formed. Note the numerical error in Fig. 8 (a) just before t = 1. This is due to the steep wave profile and its interaction with the left boundary propagating back into the computational domain. This phenomenon is discussed in detail in Sec. 6.3.   Like in the λ = 0 case this is a curvature singularity. It is much easier to see this in the λ > 0 case as the Weyl invariant I diverges to positive infinity.

Two waves
We now present results pertaining to the scattering of two colinearly polarized gravitational waves. The setup is analogous to the single wave case of Sec. 5 with the exception of the boundary condition for Ψ 4 , which is now taken to be Ψ 4 (t, −1) = p(u(t)). We continue to use the gauge F = ρ − ρ which corresponds to the gauge used in the Khan-Penrose solution for colliding colinearly polarized impulsive gravitational plane waves with λ = 0 [7]. It is found that many features are similar to the case of one wave.

Comparison against λ = 0
The general behaviour can be explained by looking at contour plots of I in Fig. 10 for varying λ (so that we can see how λ > 0 differs from λ = 0) and fixing a = 1 in the wave profiles. If λ is small enough (λ = 0 or λ = 0.06), we obtain a future curvature singularity. As λ gets larger (λ = 0.6), the expansion increases the time before this singularity occurs. If we increase λ more (λ = 6), we get to the situation where the expansion has wiped out the waves and the effect of their scattering on the curvature, and we asymptote back to dS again.

Critical behaviour
Now we fix λ = 0.6 and see how varying the area of the wave profile a affects things. We find the following three scenarios, where a 1 and a 2 are given later in the section, and are found with binary search: 1 a a 1 : asymptote back to dS.
Only in case 3 do ρ, ρ , σ, σ diverge, in the other two they asymptote back to their initial values. Due to the evolution equation √ 2∂ t A = (µ +μ)A, in cases 2 and 3 we have that A → ∞ also, causing the t, z portion of the line element to approach dt 2 , causing an infinite contraction in the z-direction. This is represented in l a and n a as shown in Fig. 12, where the t = constant surfaces approach being null. Further, as we discovered in Sec. 3.1, the real part of µ is essentially the acceleration of the unit conormal to the z = constant surfaces and the fact that this acceleration diverges to negative infinity agrees with the contraction in this direction.
We are in the Gauß gauge, and along spatially constant curves, which are in this case geodesics, the proper time and the time t are equivalent. Our gauge can then be thought of as adapted to free falling observers. This then lends the physical interpretation of the caustic singularity in case 2. The three possible futures occurring after the interaction of the gravitational waves with these observers can then be described as follows: • Case 1: The gravitational contraction is not strong enough to cause the timelike geodesics to converge or the curvature to diverge.
• Case 2: The gravitational contraction is strong enough to cause the timelike geodesics to converge and create a coordinate singularity. However, it is not strong enough to cause the curvature to diverge and this goes back to zero.
• Case 3: The gravitational contraction is strong enough to cause both the timelike geodesics to converge and the curvature to diverge, resulting in a physical curvature singularity. Figure 12: The effect of A → ∞ as t → ∞ on the null vectors along a z = constant curve.
It is noted that in the gauge B = 0 the characteristic speeds of the waves are ±A. In the cases where A → ∞ we decrease the CFL number c dynamically to avoid instabilities and settle with smaller timesteps instead.
As we now have two bifurcations, which we call a 1 and a 2 , it remains to be seen whether these will have critical behaviours.
We find, again using binary search, that approximately a 1 ≈ 0.852548. Unlike in Sec. 5.3, the constraints do not diverge as our simulations use a wave profile with area a closer to a 1 . Fig. 13 and Fig. 14 show that the expansion rate drops to around 25% of its original value at its minimum, for a long time, before asymptoting back to dS again. This implies that we can cause, with just two colliding waves, the expansion rate to locally decrease substantially for a certain period, without causing a future singularity. Further, it is noticed that although µ differs substantially in the above cases, ρ, ρ , σ and σ change very little, and if drawn differ by an amount smaller than the drawn curve.  We find, again using binary search, that approximately a 2 ≈ 0.9595, and find that taking a close to this value results in the constraints remaining well behaved. Fig. 15 shows that the Weyl invariant I diverges for a > a 2 , goes to zero for a < a 2 and goes to some other value when a ≈ a 2 . In all these cases µ diverges to infinity and thus so does A. This implies that to maintain a stable evolution our timestep must decrease to compensate, and the simulations shown in Fig. 15 stop when the timestep becomes smaller than 1e-8. It is likely that the simulation with a = 0.9595 does not converge to some constant value other than zero, but rather we cannot march in time far enough to see it either diverge to infinity or approach zero.

Impulsive waves
As in Sec. 5.4 we mimic the Dirac delta function wave profiles of the λ = 0 solutions. For colliding waves, this is when Ψ 0 = δ(v) and Ψ 4 = δ(u). We thus choose our wave profiles as Ψ 0 (v, 1) = q(v), Ψ 4 (u, −1) = q(u), where q(x) is given in Eq. (31) and approximates the Dirac delta function. Our results in Sec. 6.2 indicate that we should explore three possible regions, namely regions where we asymptote back to dS, µ diverges but not I, and where I diverges. These still exist for the approximately impulsive wave profiles and are exemplified by choosing λ = 6, 0.72 and 0.6 respectively.  In particular, we see that they do not converge to zero for u, v > 0 as a → ∞ as in the case of one wave. This is to be expected from comparison with the Khan-Penrose solution, which already has non-vanishing Ψ 0 , Ψ 2 and Ψ 4 in the region after scattering, as well as a theorem by Szekeres [29].
Of particular note is the abrupt change in sign of the first time derivative of Ψ 4 for λ = 0.72. This sharp turn, which is smooth with a small enough timestep, does not appear this distinctly in any other system variables, except for Ψ 0 due to symmetry. Fig. 17 shows that this turning point occurs not only at some point along u = v but along an entire null surface which follows the characteristic of Ψ 4 from the point where the left boundary hits v = 0. This is the result of the vanishing boundary condition for Ψ 4 on the left boundary being in disagreement with the non-vanishing Ψ 4 tail generated by Ψ 0 as it passes through the boundary. While at first sight it makes sense to impose a no ingoing radiation condition, this is blatantly unphysical when the evolution itself creates ingoing modes. Between the boundary condition and the evolution equation it is the latter which is fundamental. The boundary condition is nearly completely free to choose and is put in "by hand". A common question in a non-linear regime with boundaries that contain both ingoing and outgoing modes is then: How does one make consistent the "corner condition", i.e. the physical compatibility between data on a timeslice induced via evolution and the boundary data to yield a physically meaningful result? The answer is simply that there is no clear way to prescribe boundary conditions that match the values in the interior unless one already has an exact solution.

Suppressing the expansion with a train of waves
In [2], the rate of change of an expansion rate parameter H T W with respect to time on a space-like initial value surface (IVS) was calculated to be N (H 2 − (1/3)K ab K ab ), where N is the lapse in their coordinate system, K ab is the extrinsic curvature to the IVS and H is as per our definition of dS in inflationary coordinates. They hypothesize that there should be no reason why initial data cannot be chosen to satisfy K ab K ab > 3H 2 so that the expansion is slowed down and even completely halted ‡. We can investigate this numerically without solving the constraints by simply choosing dS initial data together with a variety of boundary conditions and seeing how the space-time evolves. We thus explore how a train of waves, generated by choosing the boundary conditions for Ψ 0 and Ψ 4 appropriately, might accomplish this.
To do so, we fix λ = 0.6 and define a new function where a = 0.894, b = 3129π/128000, c = 1/3 and we choose Ψ 0 (v, 1) = p(v), Ψ 4 (u, −1) = p(u). These constants were chosen through trial and error to give the ‡ N is a lapse and should always be positive.
largest decrease in the expansion while maximizing the interval of time this occurred, before either a singularity is formed or the space-time starts to approach dS again. The cosine factor has the effect of decreasing the amplitude of the wave until it completely vanishes at ct 2 = π/2. This is to hold off a future singularity forming, while still decreasing the expansion rate H.     19 shows just how long we can decrease the expansion for, while holding off forming a singularity. Note that along a spatially constant curve the simulation time t is the proper time of a free falling observer along this curve, and so when we talk about trying to maximize the length of time before a singularity occurs, it is inherently physical. In previous sections we have found that if a singularity was to form after a collision of two waves (with our wave profile), this happens after a few t. However here, we can decrease the expansion considerably, without forming a singularity, for up to t = 13. We cannot however find boundary conditions that lower the expansion rate to zero without very quickly forming a singularity. It is certainly possible that such finely tuned boundary conditions exist, but our studies suggest that they would be very special.

Summary and discussion
In this paper we have put forward the full non-linear Einstein equations with cosmological constant and non-vanishing energy momentum with the assumption of plane symmetry. These equations were realized through the Newman-Penrose formalism and the imposition of the Friedrich-Nagy gauge, leading to a wellposed initial boundary value problem with timelike boundaries. We specialized to vacuum where λ > 0 and chose initial data to be that induced by the de Sitter space-time in inflationary coordinates. This allowed the exploration of how this space-time is affected by gravitational perturbations, which we generated through appropriate boundary conditions for Ψ 0 and Ψ 4 .
It was found that when only one of the waves was non-vanishing the space-time either wiped out the wave via expansion, or the wave was too strong and a future singularity was produced. The bifurcation was studied and did not produce any critical behaviour. The wave profile was taken to approximate the Dirac delta function to analogize with a known exact solution for λ = 0.
With both waves non-vanishing, and in the physically motivated Gauß gauge, we found three distinct situations: The waves were not strong enough to cause a contraction of our timelike curves to create a singularity, a coordinate singularity is formed but the curvature remains finite, or a curvature singularity is formed. The second case shows that we can create a singularity where the Weyl invariant I does not diverge, but our expansion parameter diverges to negative infinity along, and close to, the surface u = v. The critical behaviour of the two bifurcations separating these futures was explored. Impulsive wave profiles were approximated and it was shown that two bifurcations occur in this case as well.
We encountered two numerical pitfalls during our exploration. Firstly, our free evolution resulted in a false steady state solution close to the bifurcation of the single wave case. As we chose our wave area closer to the bifurcation value, our free evolution approached a steady state (while A was still evolving in time) that did not satisfy the constraints. This happened even though the constraints were satisfied and converged above and below this critical value, showing how careful one must be in monitoring constraints during a free evolution. Secondly, we found that the combination of the evolution system and our non-radiating boundary conditions became unphysical in the colliding wave case after the waves left the computational domain through the boundaries. This was due to the backreaction of the waves creating tails of ingoing radiation, at odds with the boundary conditions. This is however, independent of the fact that our system is wellposed and numerically stable. The question as to how one could "guess" the right boundary conditions is delicate and creates a problem that all non-linear simulations, in particular in numerical relativity, face.
We presented how the above situations affected the local expansion rate H, which was taken to be the mean extrinsic curvature of our timeslices up to a constant factor. It was shown that for the case of two waves colliding, we could decrease H substantially for a long period of time, where the cut-off was determined by numerical limitations, before the space-time asymptoted back to dS. We could do a similar thing with a continuous stream of waves, making the expansion rate drop more uniformly over the computational domain. This showcased the potential to lower the expansion rate over a wider spatial interval.
Although we were not able to find boundary conditions that completely halted expansion for a period before either asymptoting back to de Sitter space-time or forming a singularity, we could still lower it substantially for a long time. Even so, our results do not violate the hypothesis of Woodard and Tsamis', namely that our universe may be in an unstable gravitationally bound state. It would be interesting to see whether, with further testing, we may be able to find boundary conditions that do completely halt expansion for a period. Now that exploration toward the behaviour of plane gravitational waves with λ > 0 has started and details have been uncovered, it would also be interesting to see whether one can use the results as hints toward an exact solution for impulsive waves. For the case of one propagating impulsive wave, knowing that the Weyl components vanish in the region after the wave has passed should already be a good start.