Robustness of slow contraction to cosmic initial conditions

Tetrad formulations of the Einstein-scalar field equations assign each spacetime point a family of unit basis 4-vectors (or vierbeins) $\{e_{0},e_{1},e_{2},e_{3}\}$ (as opposed to coordinates $\{x_{0},x_{1},x_{2},x_{3}\}$) that describe local Lorentz frames with the spacetime metric being given by the dot product of the basis vectors. The starting point is a timelike vector field $e_{0}$ that defines a future-directed timelike reference congruence, to which it is tangent. It is supplemented by a triad of spacelike unit 4-vectors $\{e_{1},e_{2},e_{3}\}$ that lie in the rest 3-spaces of $e_{0}$.

In an orthonormal tetrad formulation that we shall employ, the spacetime metric is everywhere given by

$$e_{\alpha}\cdot e_{\beta}=\eta_{\alpha\beta},$$

(2.1)

where $\eta_{\alpha\beta}={\rm diag}(-1,1,1,1)$ the Minkowski metric and "$\cdot$" is the spacetime inner product. Throughout, spacetime indices $(0-3)$ are Greek and spatial indices (1-3) are Latin. The beginning of the alphabet ($\alpha,\beta,\gamma$ or $a,b,c$) is used for tetrad indices and the middle of the alphabet ($\mu,\nu,\rho$ or $i,j,k$) is used for coordinate indices. Tetrad frame indices are raised and lowered with $\eta_{\alpha\beta}$.

The geometric variables of the formulation are the sixteen tetrad vector components and the twenty-four Ricci rotation coefficients

$$\gamma_{\alpha\beta\gamma}\equiv e_{\alpha}\cdot\nabla_{\gamma}e_{\beta}=-\gamma_{\beta\alpha\gamma},$$

(2.2)

which define the deformation of the tetrad when moving from point to point. Here, $\nabla_{\gamma}$ is the spacetime covariant derivative projected onto a tetrad $e_{\gamma}{}^{\lambda}\,\nabla_{\lambda}$. Note that the $\gamma_{\alpha\beta\gamma}$ are the ‘tetrad components’ of the Christoffel symbols $\Gamma_{\mu\nu\rho}$.

Similar to coordinate-based formulations of the field equations, the tetrad formulation greatly simplifies when making a space-time split. Unlike in the 3+1 (coordinate based) ADM formulation, where the split is defined by the constant-time spacelike hypersurface, here the split is relative to the timelike congruence defined by $e_{0}$. Note, though, that the timelike congruence does not uniquely define the auxiliary spatial congruence. Rather, the spatial triad vectors are fixed by imposing gauge conditions.

To perform the spacetime split, we first write the fifteen connection coefficients that have at least one timelike index in terms of 3-dimensional quantities $b_{a},\Omega_{a}$, and $K_{ab}$, reflecting the antisymmetry of $\gamma_{\alpha\beta\gamma}$ in its first two indices:

	$\displaystyle\gamma_{a00}$	$\displaystyle=$	$\displaystyle-\gamma_{0a0}=b_{a},$		(2.3)
	$\displaystyle\gamma_{ab0}$	$\displaystyle=$	$\displaystyle-\gamma_{ba0}=\epsilon_{abc}\Omega^{c},$		(2.4)
	$\displaystyle\gamma_{0ab}$	$\displaystyle=$	$\displaystyle-\gamma_{a0b}=-K_{ba};$		(2.5)

where $\epsilon_{abc}$ is the Levi-Civita symbol. As shown by Ehlers et al. [9, 19], the fifteen quantities describe kinematic fields associated with the timelike congruence tangent to $e_{0}$: the 3-vector $b_{a}$ is the local proper acceleration; the 3-vector $\Omega_{a}$ is the local angular velocity of the space-like triads $\{e_{1},e_{2},e_{3}\}$ relative to Fermi-propagated axes; and $K_{ba}$ is the local rate-of-strain (or shear) tensor. (A spatial triad is ‘Fermi-propagated’ if $\Omega_{a}\equiv 0$, i.e., it is a local, inertially non-rotating reference frame.)

Similarly, we express the remaining nine purely spatial connection coefficients $\gamma_{abc}$ that describe the induced curvature of the auxiliary 3-congruence using a 3-tensor,

$$N_{ab}=\textstyle{\frac{1}{2}}\epsilon_{b}{}^{cd}\gamma_{cda}.$$

(2.6)

The spatial connection coefficients $N_{ab}$ and the components of the shear tensor $K_{ab}$ are the eighteen dynamical variables. The acceleration and angular velocity vectors, $b_{a},\Omega_{a}$, are the six (tetrad) gauge source functions.

For completeness, we note that some authors use as basic variables the commutation (or structure) coefficients ${\cal C}_{\alpha\beta\gamma}$ rather than the Ricci rotation coefficients $\gamma_{\alpha\beta\gamma}$ [10, 28, 29]. Here, the ${\cal C}_{\alpha\beta\gamma}$ are given through

$$[e_{\alpha},e_{\beta}]\equiv\nabla_{\alpha}e_{\beta}-\nabla_{\beta}e_{\alpha}\equiv{\cal C}_{\alpha\beta\gamma}e^{\gamma}.$$

(2.7)

It is straightforward to relate the two conventions: With the definitions in Eqs. (2.2) and (2.7),

$$\gamma_{\alpha\beta\gamma}=\frac{1}{2}\Big{(}{\cal C}_{\gamma\beta\alpha}+{\cal C}_{\alpha\gamma\beta}-{\cal C}_{\beta\alpha\gamma}\Big{)},$$

(2.8)

or, alternatively,

$${\cal C}_{\alpha\beta\gamma}=\gamma_{\gamma\beta\alpha}-\gamma_{\gamma\alpha\beta},$$

(2.9)

i.e., expressed in terms of the twenty-four 3D quantities, the full set of commutation coefficients takes the following form,

	$\displaystyle{\cal C}_{a00}$	$\displaystyle=$	$\displaystyle-{\cal C}_{0a0}=b_{a},$		(2.10)
	$\displaystyle{\cal C}_{ab0}$	$\displaystyle=$	$\displaystyle-{\cal C}_{ba0}=2\epsilon_{abc}\omega^{c},$		(2.11)
	$\displaystyle{\cal C}_{0ab}$	$\displaystyle=$	$\displaystyle-{\cal C}_{a0b}=-K_{(ab)}+\epsilon_{abc}\left(\omega^{c}-\Omega^{c}\right),$		(2.12)
	$\displaystyle{\cal C}_{abc}$	$\displaystyle=$	$\displaystyle-{\cal C}_{bac}=\epsilon_{cb}{}^{d}N_{ad}.$		(2.13)

Here, the 3-vector $\omega_{a}$ is the antisymmetric part of $K_{ab}$,

$$\omega_{a}\equiv{\textstyle\frac{1}{2}}\epsilon_{a}{}^{bc}K_{bc};$$

(2.14)

it measures the vorticity (or twist) vector of the $e_{0}$-congruence.

Finally, the geometric variables have to be supplemented by the dynamical quantities that describe the matter source. In our case, this is the canonical scalar field $\phi$ which is specified in the Einstein-scalar field equations through its potential energy density $V(\phi)$.

Next we present the tetrad evolution and constraint equations for the (eighteen) dynamical variables $K_{ab},N_{ab}$. The remaining (six) gauge variables are fixed by our tetrad frame gauge choice.

A natural gauge choice is a frame with

1. i.

   Fermi-propagated axes ($\Omega_{a}\equiv 0$); and

2. ii.

   hypersurface orthogonal timelike congruence ($\omega_{a}\equiv 0$ or, equivalently, $K_{ab}\equiv K_{(ab)}$).

Here and throughout, parentheses denote symmetrization, i.e., $K_{(ab)}\equiv{\textstyle\frac{1}{2}}(K_{ab}+K_{ba})$. In this (frame) gauge, the time-like vierbein $e_{0}$ is the future-directed unit normal to the spacelike hypersurfaces $\Sigma_{t}$ of constant time, and the spatial tetrad vectors are tangent to $\{\Sigma_{t}\}$. Furthermore, $K_{ab}$ is the extrinsic curvature of $\Sigma_{t}$ and the $N_{ab}$ are the nine (intrinsic) spatial curvature variables. All Ricci rotation coefficients, $K_{ab},N_{ab}$ act as scalars on $\Sigma_{t}$.

Employing these gauge conditions, the tetrad evolution and constraint equations take the following form:

	$\displaystyle D_{0}K_{ab}$	$\displaystyle=$	$\displaystyle\epsilon_{a}{}^{cd}D_{c}N_{db}+D_{a}b_{b}+b_{a}b_{b}-\epsilon_{b}{}^{cd}N_{ac}b_{d}+N_{c}{}^{c}N_{ab}-K_{a}{}^{c}K_{cb}-N_{ca}N^{c}{}_{b}$
		$\displaystyle+$	$\displaystyle{\textstyle\frac{1}{2}}\epsilon_{a}{}^{df}\epsilon_{b}{}^{ce}\left(K_{dc}K_{fe}-N_{dc}N_{fe}\right)+s_{ab}-{\textstyle\frac{1}{2}}\delta_{ab}\big{(}\varrho+3p\big{)},$
	$\displaystyle D_{0}N_{ab}$	$\displaystyle=$	$\displaystyle-\epsilon_{a}{}^{cd}D_{c}K_{db}+\epsilon_{b}{}^{cd}K_{ac}b_{d}-N_{c}{}^{c}K_{ab}+2N_{c[a}K_{b]}{}^{c}+\epsilon_{a}{}^{df}\epsilon_{b}{}^{ce}N_{dc}K_{fe}$
		$\displaystyle-$	$\displaystyle\epsilon_{ab}{}^{c}\,j_{c},$
	$\displaystyle 2D_{b}A^{b}$	$\displaystyle=$	$\displaystyle N_{ab}N^{ab}+{\textstyle\frac{1}{2}}\left(K_{ab}K^{ab}-N_{ab}N^{ab}-(K_{a}{}^{a})^{2}-(N_{a}{}^{a})^{2}\right)+\varrho,$		(2.17)
	$\displaystyle D_{b}K_{a}{}^{b}$	$\displaystyle-$	$\displaystyle D_{a}K_{c}{}^{c}=\epsilon_{a}{}^{bc}K_{b}{}^{d}N_{dc}+2K_{a}{}^{c}A_{c}-j_{a},$		(2.18)
	$\displaystyle D_{b}N_{a}{}^{b}$	$\displaystyle-$	$\displaystyle D_{a}N_{c}{}^{c}=-\epsilon_{a}{}^{bc}N_{b}{}^{d}N_{dc},$		(2.19)

where

$$A_{b}\equiv{\textstyle\frac{1}{2}}\epsilon_{b}{}^{cd}N_{cd}.$$

(2.20)

is the antisymmetric part of $N_{ab}$; $D_{0}$ is the Lie derivative along the timelike vierbein $e_{0}$; and $D_{a}$ denotes the directional derivative along the spatial vierbein $e_{a}$. The matter variables associated with the stress-energy $T_{\alpha\beta}$, such as the energy density $\varrho$, pressure $p$, 3-momentum flux $j_{a}$, and (spatial) stress tensor $s_{ab}$, are defined as follows:

	$\displaystyle\varrho$	$\displaystyle\equiv$	$\displaystyle e_{0}{}^{\alpha}e_{0}{}^{\beta}T_{\alpha\beta},$		(2.21)
	$\displaystyle j_{a}$	$\displaystyle\equiv$	$\displaystyle-e_{0}{}^{\alpha}e_{a}{}^{\beta}T_{\alpha\beta},$		(2.22)
	$\displaystyle s_{ab}$	$\displaystyle\equiv$	$\displaystyle e_{a}{}^{\alpha}e_{b}{}^{\beta}T_{\alpha\beta},$		(2.23)
	$\displaystyle p$	$\displaystyle\equiv$	$\displaystyle{\textstyle\frac{1}{3}}s_{c}{}^{c}.$		(2.24)

Notably, there is no evolution equation for the acceleration 3-vector $b_{a}$, which reflects the fact that it is a (frame) gauge source function.

The tetrad evolution and constraint equations (2.2-2.19) were first obtained in Ref. [7]. Their derivation is straightforward when using the tetrad form of the Riemann tensor,

$$R_{\alpha\beta\gamma\delta}=D_{\gamma}\gamma_{\alpha\beta\delta}-D_{\delta}\gamma_{\alpha\beta\gamma}+\gamma_{\alpha\epsilon\gamma}\gamma^{\epsilon}{}_{\beta\delta}-\gamma_{\alpha\epsilon\delta}\gamma^{\epsilon}{}_{\beta\gamma}+\gamma_{\alpha\beta\epsilon}\left(\gamma^{\epsilon}{}_{\gamma\delta}-\gamma^{\epsilon}{}_{\delta\gamma}\right).$$

(2.25)

More exactly, substituting Eq. (2.25) into the (trace-reversed) Einstein-scalar field equations,

$$R_{\alpha\beta}=T_{\alpha\beta}-{\textstyle\frac{1}{2}}\eta_{\alpha\beta}\eta^{\gamma\delta}T_{\gamma\delta},$$

(2.26)

where $R_{\alpha\beta}\equiv R^{\gamma}{}_{\alpha\gamma\beta}$ is the Ricci tensor, yields the evolution equation (2.2) for $K_{ab}$ as well as the Hamiltonian and momentum constraints, Eqs. (2.17) and (2.18), respectively. The evolution and constraint equations for the spatial curvature variables $N_{ab}$, Eqs. (2.2, 2.19), follow from the Riemann identities,

	$\displaystyle R_{\alpha\beta\gamma\delta}=R_{\gamma\delta\alpha\beta},$		(2.27)
	$\displaystyle R_{\alpha\beta\gamma\delta}+R_{\alpha\delta\beta\gamma}+R_{\alpha\gamma\delta\beta}=0.$		(2.28)

For a single scalar field with canonical kinetic energy density and (non-zero) potential $V(\phi)$, the hydrodynamical (macroscopic) matter variables $\varrho,p,j_{a}$ and $s_{ab}$ are given by

	$\displaystyle\varrho$	$\displaystyle=$	$\displaystyle{\textstyle\frac{1}{2}}D_{0}\phi D_{0}\phi+{\textstyle\frac{1}{2}}D^{a}\phi D_{a}\phi+V(\phi),$		(2.29)
	$\displaystyle s_{ab}$	$\displaystyle=$	$\displaystyle D_{a}\phi D_{b}\phi+\left({\textstyle\frac{1}{2}}D_{0}\phi D_{0}\phi-D_{c}\phi D^{c}\phi-V(\phi)\right)\delta_{ab},$		(2.30)
	$\displaystyle p$	$\displaystyle=$	$\displaystyle{\textstyle\frac{1}{2}}D_{0}\phi D_{0}\phi-{\textstyle\frac{1}{6}}D^{a}\phi D_{a}\phi-V(\phi),$		(2.31)
	$\displaystyle j_{a}$	$\displaystyle=$	$\displaystyle-D_{0}\phi D_{a}\phi.$		(2.32)

Note that, in general, $j_{a}$ is non-zero and $s_{ab}$ is non-diagonal, which reflects the fact that choosing a hypersurface-orthogonal tetrad frame gauge generally does not coincide with the co-moving frame of the (scalar field) matter source. In fact, a hypersurface-orthogonal frame-gauge is co-moving only in the homogeneous (ultra-local) limit.

The system (2.2-2.19) is completed by adding the scalar-field evolution and constraint equations:

	$\displaystyle D_{0}\phi$	$\displaystyle=$	$\displaystyle W,$		(2.33)
	$\displaystyle D_{a}\phi$	$\displaystyle=$	$\displaystyle S_{a},$		(2.34)
	$\displaystyle D_{0}W$	$\displaystyle=$	$\displaystyle-\delta^{ab}K_{ab}W+D_{a}S^{a}+\left(b_{a}-2A_{a}\right)S^{a}-V,_{\phi},$		(2.35)
	$\displaystyle D_{0}S_{a}$	$\displaystyle=$	$\displaystyle D_{a}W+b_{a}W-K_{(ab)}S^{b}.$		(2.36)

Here, Eqs. (2.33) and (2.34) are the defining relations for the auxiliary variables $W$ and $S_{a}$, which denote the velocity and gradient of the scalar field $\phi$, respectively. Eq. (2.35) is obtained by expanding the Laplacian of the Klein-Gordon equation ($\Box\phi=V,_{\phi}$) using the Ricci rotation coefficients; and Eq. (2.36) is obtained by evaluating the commutation relation $[e_{0},e_{a}]\phi=-{\cal C}_{0a0}W+{\cal C}_{0ab}S^{b}$ using Eqs. (2.10-2.13).

In order to evolve the tetrad equations (2.2-2.19, 2.33-2.36) numerically, we must write them as a system of partial differential equations. That means, we must give a representation of the tetrad vector components $\{e_{\alpha}\}$ using a particular set of local coordinates $\{x^{\mu}\}$ and then convert the directional derivatives $D_{\alpha}$ in the tetrad equations to partial derivatives along these coordinates.

To this end, we introduce the transformation matrix $\{\lambda_{\alpha}^{\mu}\}$ between coordinate and tetrad basis vectors defined thru

$$e_{\alpha}=\lambda_{\alpha}^{\mu}e_{\mu}.$$

(2.37)

The coordinate metric components are then given by

$$g^{\mu\nu}=\eta^{\alpha\beta}\lambda_{\alpha}^{\mu}\lambda_{\beta}^{\nu};$$

(2.38)

and directional derivatives along tetrads can now be written as partial derivatives along coordinate directions,

$$D_{0}=N^{-1}\left(\partial_{t}-N^{i}\partial_{i}\right),\quad D_{a}=E_{a}{}^{i}\partial_{i},$$

(2.39)

where $N$ is the tetrad lapse function and the $N^{i}$ are the three coordinate components of the tetrad shift vector. Both the tetrad lapse function and the tetrad shift vector describe the evolution of the coordinates relative to the tetrad congruence (as opposed to the ADM lapse and shift that describe the evolution of the proper time and co-moving spatial coordinates relative to the coordinates of a particular foliation). The nine coordinate components $E_{a}{}^{i}$ describe projections of the spatial tetrads tangent to the constant-time hypersurface $\Sigma_{t}$. Note that the spatial triad vectors have zero time component, since we have chosen our tetrad frame gauge to be hypersurface-orthogonal. (For arbitrary tetrad frame gauge choices, this is not the case.)

The $E_{a}{}^{i}$ are dynamical variables determined by the evolution and constraint equations

	$\displaystyle N^{-1}\partial_{t}E_{a}{}^{i}$	$\displaystyle=$	$\displaystyle-K_{a}{}^{c}E_{c}{}^{i},$		(2.40)
	$\displaystyle\epsilon_{c}{}^{ab}E_{a}{}^{i}\partial_{i}E_{b}{}^{j}$	$\displaystyle=$	$\displaystyle N^{d}{}_{c}E_{d}{}^{j}-N^{d}{}_{d}E_{c}{}^{j},$		(2.41)

which are derived from applying the commutators of the basis vectors to the spatial coordinates $\{x^{i}\}$.

The lapse function $N$ and the shift vector $N^{i}$ are gauge variables that we can freely specify. A natural coordinate gauge choice when studying the evolution of cosmological spacetimes is to have co-moving coordinates, such that the $x^{i}$ are constant along the congruence and

$$N^{i}\equiv 0.$$

(2.42)

For the lapse function, we will impose the condition that hypersurfaces $\Sigma_{t}$ of constant time be constant mean curvature (CMC) hypersurfaces, i.e., the trace of the extrinsic curvature $K_{ab}$ is spatially uniform for each $\Sigma_{t}$,

$$\Theta^{-1}\equiv-{\textstyle\frac{1}{3}}K_{a}{}^{a}\big{|}_{\Sigma_{t}}={\rm const}>0.$$

(2.43)

Imposing CMC slicing, the trace of Eq. (2.2) reduces to a linear, elliptic equation for the lapse $N$,

$$\Big{(}-\big{(}D_{a}-2A_{a}\big{)}D^{a}+\Sigma_{ab}\Sigma^{ab}+3\Theta^{-2}+W^{2}-V(\phi)\Big{)}N=0,$$

(2.44)

where

$$\Sigma_{ab}\equiv K_{ab}-{\textstyle\frac{1}{3}}K_{c}{}^{c}\delta_{ab}$$

(2.45)

is the trace-free part of the extrinsic curvature. Few works that rigorously prove well-posedness treat elliptic gauge conditions, and we are not aware of any for the particular tetrad formulation we use; however, see Ref. [1] for a proof in a closely related coordinate based formulation. We note, though, that if a system of partial differential equations is not well posed it would be challenging, if not impossible, to obtain stable, convergent numerical solutions with any discretization scheme. That our code is stable and convergent can thus be viewed as ‘empirical evidence’ that the underlying equations are well posed.

Note that the (elliptic) equation for the (tetrad) lapse also determines the tetrad gauge source function $b_{a}$, which denotes the acceleration of the tetrad congruence worldlines. This can be seen, e.g., by computing the commutator $[e_{0},e_{a}]$ as applied to the time coordinate $x^{0}$,

$$b_{a}=N^{-1}E_{a}{}^{i}\partial_{i}N.$$

(2.46)

For the time coordinate $t$, we choose a particular (re-)scaling,

$${e^{t}}=3\Theta,$$

(2.47)

that is consistent with the CMC slicing condition. If Eq. (2.43) is satisfied, the inverse trace of the extrinsic curvature (here and throughout denoted by $\Theta$) is the well-known Hubble radius (as measured by the proper time coordinate $\tau$), i.e.,

$${\textstyle\frac{1}{3}}e^{-t}=\frac{d\ln a(\tau)}{d\tau}.$$

(2.48)

During contraction, the Hubble radius decreases. Accordingly, the time coordinate $t\leq 0$, running from small negative to large negative values. Yet, due to CMC slicing, all curvature variables remain finite and non-zero for any evolution of finite duration. This is a necessary condition to stably evolve contracting phases that last several hundreds of $e$-folds, which is required to study the robustness of slow contraction.

In addition, for a sufficiently long evolution, we must ensure that no two dynamical variables grow (or shrink) at significantly different rates to avoid a stiffness problem. Since in the case of slow contraction, some spatial metric variables, such as the spatially averaged scale factor, decrease at a significantly lower rate than some curvature variables, such as the Hubble radius, we eliminate the latter by introducing dimensionless Hubble-normalized variables,

	$\displaystyle N$	$\displaystyle\to$	$\displaystyle{\cal N}\equiv N/\Theta,$		(2.49)
	$\displaystyle\Big{\{}E_{a}{}^{i},\Sigma_{ab},A_{a},n_{ab},W,S_{a},V\Big{\}}$	$\displaystyle\to$	$\displaystyle\Big{\{}{\bar{E}}_{a}{}^{i},{\bar{\Sigma}}_{ab},{\bar{A}}_{a},{\bar{n}}_{ab},{\bar{W}},{\bar{S}}_{a},{\bar{V}}\Big{\}},$		(2.50)

where ${\cal N}$ is the Hubble-normalized lapse function; bar denotes multiplication by the Hubble radius $\Theta$; and

$$n_{ab}\equiv N_{(ab)}$$

(2.51)

is the symmetric part of $N_{ab}$. Substituting into Eq. (2.44), yields an elliptic equation for the Hubble-normalized lapse ${\cal N}$

$\displaystyle-$

$\displaystyle\bar{E}^{a}{}_{i}\partial^{i}\left(\bar{E}_{a}{}^{j}\partial_{j}{\cal N}\right)+2\bar{A}^{a}\bar{E}_{a}{}^{i}\partial_{i}{\cal N}+{\cal N}\left(3+\bar{\Sigma}_{ab}\bar{\Sigma}^{ab}+\bar{W}^{2}-\bar{V}\right)=3\,.$

(2.52)

Imposing CMC slicing as defined in Eqs. (2.43, 2.47) and using Hubble-normalized variables in Eqs. (2.2-2.2, 2.33-2.36, and 2.40), the gravitational quantities ${\bar{E}}_{a}{}^{i},{\bar{\Sigma}}_{ab},{\bar{n}}_{ab}$, and ${\bar{A}}_{a}$ as well as the scalar field matter variables $\phi,\bar{W},\bar{S}_{a}$ satisfy the hyperbolic evolution equations

	$\displaystyle\partial_{t}\bar{E}_{a}{}^{i}$	$\displaystyle=$	$\displaystyle-\Big{(}{\cal N}-1\Big{)}\bar{E}_{a}{}^{i}-{\cal N}\,\bar{\Sigma}_{a}{}^{b}\bar{E}_{b}{}^{i},$		(2.53)
	$\displaystyle\partial_{t}\bar{\Sigma}_{ab}$	$\displaystyle=$	$\displaystyle-\Big{(}3{\cal N}-1\Big{)}\bar{\Sigma}_{ab}-{\cal N}\Big{(}2\bar{n}_{\langle a}{}^{c}\,\bar{n}_{b\rangle c}-\bar{n}^{c}{}_{c}\bar{n}_{\langle ab\rangle}-\bar{S}_{\langle a}\bar{S}_{b\rangle}\Big{)}+\bar{E}_{\langle a}{}^{i}\partial_{i}\Big{(}\bar{E}_{b\rangle}{}^{i}\partial_{i}{\cal N}\Big{)}$
		$\displaystyle-$	$\displaystyle{\cal N}\left(\bar{E}_{\langle a}{}^{i}\partial_{i}\bar{A}_{b\rangle}-\epsilon^{cd}{}_{(a}\Big{(}\bar{E}_{c}{}^{i}\partial_{i}\bar{n}_{b)d}-2\bar{A}_{c}\bar{n}_{b)d}\Big{)}\right)+\epsilon^{cd}{}_{(a}\bar{n}_{b)d}\bar{E}_{c}{}^{i}\partial_{i}{\cal N}+\bar{A}_{\langle a}\bar{E}_{b\rangle}{}^{i}\partial_{i}{\cal N},$
	$\displaystyle\partial_{t}\bar{n}_{ab}$	$\displaystyle=$	$\displaystyle-\Big{(}{\cal N}-1\Big{)}\bar{n}_{ab}+{\cal N}\Big{(}2\bar{n}_{(a}{}^{c}\bar{\Sigma}_{b)c}-\epsilon^{cd}{}_{(a}\bar{E}_{c}{}^{i}\partial_{i}\bar{\Sigma}_{b)d}\Big{)}-\epsilon^{cd}{}_{(a}\bar{\Sigma}_{b)d}\bar{E}_{c}{}^{i}\partial_{i}{\cal N},$		(2.55)
	$\displaystyle\partial_{t}\bar{A}_{a}$	$\displaystyle=$	$\displaystyle-\Big{(}{\cal N}-1\Big{)}\bar{A}_{a}-{\cal N}\Big{(}\bar{\Sigma}_{a}{}^{b}\bar{A}_{b}-{\textstyle\frac{1}{2}}\bar{E}_{b}{}^{i}\partial_{i}\bar{\Sigma}_{a}{}^{b}\Big{)}-\bar{E}_{a}{}^{i}\partial_{i}{\cal N}+{\textstyle\frac{1}{2}}\bar{\Sigma}_{a}{}^{b}\bar{E}_{b}{}^{i}\partial_{i}{\cal N},$		(2.56)
	$\displaystyle\partial_{t}\phi$	$\displaystyle=$	$\displaystyle{\cal N}\,\bar{W},$		(2.57)
	$\displaystyle\partial_{t}\bar{W}$	$\displaystyle=$	$\displaystyle-\Big{(}3{\cal N}-1\Big{)}\bar{W}-{\cal N}\Big{(}\bar{V}_{,\phi}+2\bar{A}^{a}\bar{S}_{a}-\bar{E}_{a}{}^{i}\partial_{i}\bar{S}^{a}\Big{)}+\bar{S}^{a}\bar{E}_{a}{}^{i}\partial_{i}{\cal N},$		(2.58)
	$\displaystyle\partial_{t}\bar{S}_{a}$	$\displaystyle=$	$\displaystyle-\Big{(}{\cal N}-1\Big{)}\bar{S}_{a}-{\cal N}\Big{(}\bar{\Sigma}_{a}{}^{b}\bar{S}_{b}-\bar{E}_{a}{}^{i}\partial_{i}\bar{W}\Big{)}+\bar{W}\bar{E}_{a}{}^{i}\partial_{i}{\cal N}.$		(2.59)

(Here angle brackets denote traceless symmetrization defined as $X_{\langle ab\rangle}\equiv X_{(ab)}-{\textstyle\frac{1}{3}}X_{c}{}^{c}\delta_{ab}$.) The system of Eqs. (2.52-2.59) will serve as our numerical scheme. Notably, this system is well-posed, as was shown, e.g., in [4], and, by construction, it admits stable numerical evolution that involves several hundreds of $e$-folds of contraction.

In addition, the same variables satisfy the constraint equations

	$\displaystyle{\cal C}_{\rm G}$	$\displaystyle\equiv$	$\displaystyle 3+2\bar{E}_{a}{}^{i}\partial_{a}\bar{A}^{a}-3\bar{A}^{a}\bar{A}_{a}-{\textstyle\frac{1}{2}}\bar{n}^{ab}\bar{n}_{ab}+{\textstyle\frac{1}{4}}(\bar{n}^{c}{}_{c})^{2}-{\textstyle\frac{1}{2}}\bar{\Sigma}^{ab}\bar{\Sigma}_{ab}$
		$\displaystyle-$	$\displaystyle{\textstyle\frac{1}{2}}\bar{W}^{2}-{\textstyle\frac{1}{2}}\bar{S}^{a}\bar{S}_{a}-{\bar{V}}=0\,,$
	$\displaystyle({\cal C}_{\rm C})_{a}$	$\displaystyle\equiv$	$\displaystyle\bar{E}_{b}{}^{i}\partial_{i}{\bar{\Sigma}}_{a}{}^{b}-3{\bar{\Sigma}}_{a}{}^{b}\bar{A}_{b}-\epsilon_{a}{}^{bc}\bar{n}_{b}{}^{d}\bar{\Sigma}_{cd}-{\bar{W}}{\bar{S}}_{a}=0\,,$		(2.61)
	$\displaystyle({\cal C}_{\rm J})_{a}$	$\displaystyle\equiv$	$\displaystyle\bar{E}_{b}{}^{i}\partial_{i}\bar{n}^{b}{}_{a}+\epsilon^{bc}{}_{a}\bar{E}_{b}{}^{i}\partial_{i}\bar{A}_{c}-2\bar{A}_{b}\bar{n}^{b}{}_{a}=0\,,$		(2.62)
	$\displaystyle({\cal C}_{\rm S})_{a}$	$\displaystyle\equiv$	$\displaystyle{\bar{S}}_{a}-{\bar{E}}_{a}{}^{i}\partial_{i}\phi=0\,,$		(2.63)
	$\displaystyle({\cal C}_{\rm E})_{a}^{i}$	$\displaystyle\equiv$	$\displaystyle\epsilon^{bc}{}_{a}\Big{(}\bar{E}_{b}{}^{j}\partial_{j}\bar{E}_{c}{}^{i}-\bar{A}_{b}\bar{E}_{c}{}^{i}\Big{)}-\bar{n}_{a}{}^{d}\bar{E}_{d}{}^{i}=0.$		(2.64)

The constraints are the Hubble-normalized version of Eqs. (2.17-2.19, 2.34, 2.41), again imposing CMC slicing conditions as in Eqs. (2.43, 2.47). (The subscripts $G,C$ and $J$ stand for Gauss, Codazzi, and Jacobi, respectively, referring to the commonly used terminology.) As detailed in the following two sections, we shall use the constraint equations to specify the initial conditions and to ensure constraint damping and numerical convergence.

To specify the initial conditions, we first choose a particular time $t_{0}$. With the tetrad and coordinate gauge conditions remaining the same as specified above for the evolution scheme, the $t_{0}$-hypersurface has constant mean curvature $\Theta^{-1}_{0}$, the value of which we can freely choose, and zero shift. We note that, using Hubble-normalized variables in our evolution scheme, the natural units are set by this initial Hubble radius $\Theta_{0}$ rather than the Planck scale.

Next, we must fix the variables

$$\big{\{}\bar{E}_{a}{}^{i},\bar{n}_{ab},\bar{A}_{a},\bar{\Sigma}_{ab}\big{\}}$$

(3.1)

that describe the geometry of the $t_{0}$-hypersurface. Not all of these variables are freely specifiable, though, as they must satisfy the constraint equations (2.3-2.64). (Notably, the evolution equations (2.53-2.59) propagate the constraints, i.e., ensure that the constraints are satisfied at later times provided they are satisfied on the initial time slice.)

Adapting the York method [31] commonly used in numerical relativity computations, we choose the initial metric to be conformally flat,

$$g_{ij}\equiv\psi^{4}(x,t_{0})\delta_{ij},$$

(3.2)

where $\psi$ is the conformal factor. The conformal factor is not a free function but fixed by the Hamiltonian (or Gauss) constraint (2.3), as we will detail below.

With Eq. (2.38), this choice for $g_{ij}$ simultaneously fixes the coordinate components of the spatial triad:

$${\bar{E}}_{a}{}^{i}=\psi^{-2}\Theta^{-1}_{0}\delta_{a}{}^{i}.$$

(3.3)

Substituting into the spatial commutator,

$$[\bar{e}_{a},\bar{e}_{b}]\equiv\bar{{\cal C}}_{abc}\bar{e}^{c}=\Big{(}2\bar{A}_{[a}\delta_{b]c}+\epsilon_{abd}\bar{n}^{d}{}_{c}\Big{)}\bar{e}^{c},$$

(3.4)

as defined in Eq. (2.7), we find the components of the intrinsic curvature tensor,

	$\displaystyle{\bar{n}}_{ab}$	$\displaystyle=$	$\displaystyle 0,$		(3.5)
	$\displaystyle{\bar{A}}_{a}$	$\displaystyle=$	$\displaystyle-2\psi^{-1}\bar{E}_{a}{}^{i}\partial_{i}\psi;$		(3.6)

and the constraint equations (2.62) and (2.64) are trivially satisfied by this combination of $\bar{E}_{a}{}^{i},\bar{n}_{ab},\bar{A}_{a}$.

We stress that, in general, ${\bar{A}}_{a}$ is non-zero, reflecting the fact that the anti-symmetric part of the intrinsic curvature tensor does not transform trivially under conformal rescaling, as pointed out in Ref. [4]. That means, most especially, assuming the initial metric to be conformally flat does not imply uniform spatial curvature on the initial slice and, hence, does not impose a real restriction on the initial data set.

Having set the geometric variables $\{\bar{E}_{a}{}^{i},\bar{n}_{ab},\bar{A}_{a}\}$ through specifying the spatial metric for the initial slice, it remains to determine the components of the Hubble-normalized shear tensor $\bar{\Sigma}_{ab}$, as defined in Eq. (2.45), to close the set of variables describing the geometry of the initial $t_{0}$-hypersurface.

It is a particular advantage of the conformal rescaling as suggested by the York method that the constraint equations significantly simplify. For this reason, we will determine the initial shear contribution $\bar{\Sigma}_{ab}(x,t_{0})$ by first specifying its conformally rescaled counterpart

$$Z_{ab}(x,t_{0})\equiv\psi^{6}\bar{\Sigma}_{ab}(x,t_{0}),$$

(3.7)

using the momentum constraint (2.61) as evaluated for $Z_{ab}(x,t_{0})$,

$$\bar{E}^{a}{}_{i}(x,t_{0})\,\partial^{i}Z_{ab}(x,t_{0})=Q(x,t_{0})\bar{E}_{b}{}^{i}(x,t_{0})\partial_{i}\phi(x,t_{0}).$$

(3.8)

Here, $Q(x,t_{0})$ is the conformally rescaled scalar field kinetic energy defined by

$\displaystyle Q(x,t_{0})$

$\displaystyle\equiv$

$\displaystyle\psi^{6}(x,t_{0})\bar{W}(x,t_{0}).$

(3.9)

Then, solving the Hamiltonian constraint (2.3) for the conformal factor $\psi(x,t_{0})$ yields $\bar{\Sigma}_{ab}(x,t_{0})$.

It is immediately obvious that we can follow two strategies in specifying $Z_{ab}(x,t_{0})$: either we freely specify the initial shear contribution $Z_{ab}(x,t_{0})$ that then fixes parts of the initial field and velocity distribution $\{\phi(x,t_{0}),Q(x,t_{0})\}$, or we freely specify only parts of the initial shear contribution $Z_{ab}(x,t_{0})$ so we can freely choose $\{\phi(x,t_{0}),Q(x,t_{0})\}$ that together fix the rest of $Z_{ab}(x,t_{0})$ using the momentum constraint. In this analysis, we have chosen the latter option.

By definition, the vacuum contribution $Z_{ab}^{0}(x,t_{0})$ of the conformally rescaled shear tensor, i.e., the divergence-free part of $Z_{ab}(x,t_{0})$, is independent of any matter source. Accordingly, we will freely specify only $Z_{ab}^{0}(x,t_{0})$ and constrain the rest of $Z_{ab}(x,t_{0})$ by the choice of the initial scalar field and velocity distribution $\{\phi(x,t_{0}),Q(x,t_{0})\}$.

For the numerical simulations, we use periodic boundary conditions $0\leq x\leq 2\pi$ with $0$ and $2\pi$ identified. In addition, we restrict ourselves to deviations from homogeneity along a single spatial direction $x$ so that the spacetimes have two Killing fields. Since the variables depend only on $x$ and since $x$ is periodically identified, we specify their spatial variation as a sum of Fourier modes with period $2\pi$.

A simple yet generic divergence-free and trace-free choice for $Z_{ab}^{0}(x,t_{0})$ is given by

$$Z_{ab}^{0}(x,t_{0})=\left(\begin{array}[]{ccc}{b_{2}}&\xi&0\\ \xi&a_{1}\cos x+{b_{1}}&a_{2}\cos x\\ 0&a_{2}\cos x&-{b_{1}}-{b_{2}}-a_{1}\cos x\end{array}\right),$$

(3.10)

where $\xi,\,a_{1},\,a_{2},\ b_{1}$ and $b_{2}$ are constants. Note that since $\bar{\Sigma}_{ab}$ must be trace-free, so must $Z^{0}_{ab}$. The rest of the initial shear distribution, $Z_{ab}-Z^{0}_{ab}$ is obtained by solving the momentum constraint (3.8), which turns into an algebraic relation for the Fourier coefficients of this non-zero divergence piece of $Z_{ab}$ upon specifying the initial scalar field matter variables.

In setting the initial velocity and field distribution, we freely specify the Fourier coefficients of $Q(x,t_{0})$ and $\phi(x,t_{0})$ via

	$\displaystyle\phi(x,t_{0})$	$\displaystyle=$	$\displaystyle f_{0}\cos\big{(}m_{0}x+d_{0}\big{)},$		(3.11)
	$\displaystyle Q(x,t_{0})$	$\displaystyle=$	$\displaystyle\Theta\Big{(}f_{1}\cos\big{(}m_{1}x+d_{1}\big{)}+Q_{0}\Big{)}.$		(3.12)

where $f_{0},m_{0},d_{0},f_{1},m_{1},d_{1}$, and $Q_{0}$ are constants.

To compute the initial value of the (Hubble-normalized) scalar field gradient $\bar{S}_{a}(x,t_{0})$, we substitute $\phi(x,t_{0})$ into the constraint equation (2.63).