# Benchmarks of Dark Energy

Eric V. Linder  Berkeley Center for Cosmological Physics & Berkeley Lab,
University of California, Berkeley, CA 94720, USA

###### Abstract

Dark energy is a premier mystery of physics, both theoretical and
experimental. As we look to develop plans for high energy physics over the
next decade, within a two decade view, we consider benchmarks for revealing
the nature of dark energy. We conclude, based on fundamental physical
principles detailed below, that understanding will come from experiments
reaching key benchmarks:

  * •

$\sigma(w_{a})<2.5\sigma(w_{0})$

  * •

$\sigma(w_{0})<0.02$

  * •

$\sigma(\rho_{\rm de}/\rho_{\rm crit})<(1/3)\rho_{\Lambda}/\rho_{\rm crit}$
for all redshifts $z<5$

where the dark energy equation of state $w(a)=w_{0}+w_{a}(1-a)$. Beyond the
cosmic expansion history we also discuss benchmarks for the cosmic growth
history appropriate for testing classes of gravity theories. All benchmarks
can be achieved by a robust Stage 5 program, using extensions of existing
probes plus the highly complementary, novel probe of cosmic redshift drift.

##  I Introduction

Cosmic acceleration, in its late time instantiation known as dark energy, was
highlighted at Snowmass 2001 as a premier problem, with quotes such as “the
problem of dark energy will be a roadblock on our path to a comprehensive
fundamental physical theory”, “the biggest and most profound gap in our
current understanding of the physical world”, “it is vitally important to know
if the cosmological ‘constant’ is truly constant” [1]. These comments are as
true now.

In this article, as succinctly as possible, we present the outstanding
questions about the nature of dark energy, and the benchmarks that will still
need to be met after the current (“Stage 4”) dark energy experiments complete.
It is challenging to set benchmarks – revelatory measurement levels – for a
physical phenomenon that is not understood. To do so we rely on basic physical
principles and observations, e.g. our universe is many e-folds of expansion
old, it had a long period of matter domination before the current
acceleration, etc. plus an appeal to avoidance of fine tuning to whatever
extent is viable.

Another guiding principle is that for phenomenology we choose as few
parameters as can robustly describe the impact of theory on observations, and
for theory we choose as few parameters as can showcase major characteristics
of classes of theory (in particular for modified gravity). That is, we are
looking for an alert, a signal, a direction, and want the highest significance
indication as to signatures of dark energy. This minimalist approach likely
shortchanges the rich tapestry of the dark sector, but has precedent in that
most of late time cosmic expansion and growth for example can neglect the
intricacies of the Standard Model of particle physics and simply investigate
the matter density. Once observations give key hints based on the benchmarks,
then deeper investigation will burst forth.

Section II examines the cosmic expansion history over the last 10 billion
years, while Section III motivates the investigation of dark energy before it
came to dominance. Section IV explores the cosmic growth history, the
importance of testing consistency with the expansion history, and how modified
gravity manifests. Section V provides some thoughts on how to focus
experiments (“Stage 5”) to achieve the benchmarks.

##  II Cosmic Acceleration

Cosmic expansion is determined by the Friedmann equations, giving the
expansion rate or Hubble parameter $H=\dot{a}/a$ in terms of the effective
energy-momenta of the components. For a homogeneous and isotropic universe
these resolve to the energy density and pressure of each component, and the
ratio $w_{i}(a)=P_{i}(a)/\rho_{i}(a)$ is known as the equation of state
parameter (EOS). This holds even outside general relativity by regarding the
energy density and pressure of an effective component. The late time matter
component is well treated as pressureless, $w_{m}=0$. A cosmological constant
$\Lambda$, with time independent energy density, has
$P_{\Lambda}=-\rho_{\Lambda}$ and so $w_{\Lambda}=-1$.

A central question is what is the EOS of dark energy – is it a cosmological
constant so $w_{\rm de}=-1$, does it have a different, constant value, or is
it time dependent? For a constant $w_{\rm de}$, benchmarks can be given from
basic physical principles. To cause cosmic acceleration a component must have
$w_{\rm de}<-1/3$. Topological defects such as dominating domain walls have
$w_{\rm de}=-2/3$. If we wish to distinguish between these fundamental
physical mechanisms of a cosmological constant and domain walls, we need say
$5\sigma(w_{\rm const})<({}^{-}2/3\,-\,{}^{-}1)$ or $\sigma(w_{\rm
const})<0.067$. This has been achieved by cosmic observations.

There is essentially no non-finetuned mechanism that has a constant EOS that
can cause acceleration other than topological defects, so measuring $w$ as a
constant (e.g. is $w_{\rm const}=-0.95$?) is no longer a goal. In addition,
dynamical $w(a)$ can mimic constant $w$ through the time averaging inherent
over the light propagation of a cosmic signal, so measuring $w_{\rm const}$,
even $w_{\rm const}=-1$ says very little. This is called the mirage of
$\Lambda$ [2].

Thus we need a benchmark for $w(a)$. Describing a full function with just a
small number of parameters might appear difficult, but in fact it is a solved
problem. In 2002 [3] demonstrated through study of the exact solutions of the
scalar field equation of motion (Klein-Gordon equation) that a wide variety of
models for dark energy had their exact predictions for observables (distances,
Hubble parameter) reconstructed to $\sim 0.1\%$ accuracy by a remarkably
simple form,

$w(a)=w_{0}+w_{a}(1-a)\ .$ (1)

The simplicity of the form should not obscure that this incorporates the full
physics (and does not arise from e.g. any artificial Taylor expansion).

Even models with rapid variations in $w(a)$ have their observables matched by
the $w_{0}$–$w_{a}$ parametrization: for example despite a phase transition
the vacuum metamorphosis model [4, 5] is fit by $(w_{0},w_{a})=(-1.24,-1.5)$
to 0.55% in the observables [6]; despite oscillations the viable range of
Albrecht-Skordis models [7] are fit by the $w_{0}$–$w_{a}$ parametrization to
$\lesssim 0.1\%$.

In 2008 [8] demonstrated that $w_{0}$–$w_{a}$ actually functions as a
calibration relation on the dark energy phase EOS phase space
$w$–$w^{\prime}$, where prime denotes an e-fold derivative $d/d\ln a$. It
cleanly separates different classes of physics behind dark energy. This
important property leads us to well motivated benchmarks for testing dark
energy. The two distinct classes of dark energy dynamics are freezers and
thawers [9], where the (effective) field either starts rolling or frozen,
respectively, in the high redshift, high Hubble friction universe. That is,
the field evolves toward or away from cosmological constant behavior
respectively.

The class is determined by the steepness of the (effective) potential. Each
class forms a wedge region in the $w$–$w^{\prime}$ phase space, with the outer
edges determined by basic physics: that the universe was matter dominated at
high redshift for the upper boundary of the thawing region, and that the
friction from the Hubble expansion limits the acceleration of the field for
the lower boundary of the freezing region. The constraints are well
understood, and quantified, in terms of the flow formalism [10] analogous to
that for early universe inflation.

What is of particular importance for benchmarking is that there is a gap
between the allowed regions, corresponding to the field acceleration. Only
fields that are fine tuned to have zero acceleration can live more than
briefly in this region, so there is a natural scale to how well
$w$–$w^{\prime}$ needs to be determined to separate which of the two physics
classes describes dark energy. We also emphasize that this phase space
structure has been shown to hold beyond physical scalar fields, e.g. for
modified gravity interpreted as an effective EOS.

The benchmark separation is between the lower bound of the thawing region,
along $w^{\prime}\approx 1+w$, and the upper bound of the freezing region,
along $w^{\prime}\approx w(1+w)$. See [11] for detailed discussion on the
origin of all the boundaries. Thus understanding of the physical mechanism
behind cosmic acceleration, i.e. distinction between the two classes, comes
when $\sigma(w^{\prime})<2(1+w)$, for $w\approx-1$. We can translate this into
a $w_{0}$–$w_{a}$ benchmark using the calibration relation
$w_{a}=-w^{\prime}(a_{\star})/a_{\star}$ for $a_{\star}=0.8$ from [8]. Thus
our first benchmark criterion is

$\sigma(w_{a})<2.5\sigma(w_{0})\qquad{\rm Benchmark\ Criterion\ 1.}$ (2)

Note that for Stage 4 dark energy experiments we expect roughly
$\sigma(w_{0})\approx 0.08$–0.1, $\sigma(w_{a})\approx 0.2$–0.3 including
systematics so we are close to BC1. However, this is a relative not an
absolute criterion; it does not say how accurately either individual parameter
needs to be measured.

We can go further by providing benchmarks for $w_{0}$ and $w_{a}$ separately
from either basic physics or motivated models. Recall that $w^{\prime}=dw/d\ln
a$ measures the scalar field evolution, which is driven by the steepness of
the potential and retarded by Hubble friction. Within the minimal parameter
approach (no extra scales in the problem) and in the absence of fine tuning,
the natural timescale for the scalar field to evolve on is the Hubble e-fold
scale. That is, $\dot{w}\sim H^{-1}$. This translates to $w^{\prime}\approx 1$
or $w_{a}\approx 1$. The other natural value is no evolution – that of the
cosmological constant. Thus we want to distinguish between these two natural
scales. If we ask for this important distinction be made at the $5\sigma$
level then

$\sigma(w_{a})<\frac{1}{5}\times 1=0.2\qquad{\rm Benchmark\ Criterion\ 2a.}$
(3)

This would then give $\sigma(w_{0})=0.08$, $\sigma(w_{a})=0.2$ as the
criterion, and indeed this is at the Stage 4 level. However, distinction
between evolution at a characteristic e-folding of 0 and of 1 is a relatively
modest goal. While it would give a broad impression of dark energy, it would
not yield critical clues to its specific nature.

Therefore we want to obtain more accuracy, through Stage 5 experiments. It is
more difficult to obtain a general threshold for $w_{0}$, as it can approach
arbitrarily close to the cosmological constant value of $-1$. However, an
interesting and well motivated model views cosmic acceleration in a unified
manner, across both early universe inflation and late time acceleration. This
“quintessential inflation” arises from $\alpha$-attractor high energy physics,
with connections to supergravity and string theory [12]. In particular, it
predicts a relation between $w_{0}$ and the inflation tensor-to-scalar power
ratio $r$. The parameter $r$ is a key goal of next generation cosmic microwave
background experiments detecting primordial gravitational waves through B-mode
polarization.

The $\alpha$-attractor model predicts $r=12\alpha/N^{2}$, where $N$ is the
number of inflationary e-folds, about 50–60 depending on reheating and model
details. Simultaneously it predicts that after inflation ends and the field
freezes again due to Hubble friction, late time acceleration returns as the
universe expands, the Hubble friction lessens, and the field thaws. The dark
energy then reaches an asymptotic future attractor with $w_{\rm
future}=-1+2/(9\alpha)$. Note the key promise of this theory: if $\alpha$ is
high, we will detect the gravitational wave signature in the CMB, if $\alpha$
is low we will detect the deviation of the dark energy EOS from the
cosmological constant value! An especially well motivated value is $\alpha=1$,
which corresponds to Starobinsky ($R^{2}$) inflation and Higgs inflation. This
would give $r\approx 4\times 10^{-3}$ and $w_{\rm future}=-0.78$. Today, with
$\Omega_{\rm de}\approx 0.7$ the field evolution has not fully reached the
future attractor and $1+w_{0}\approx 0.5(1+w_{\rm future})$, e.g.
$w_{0}\approx-0.9$. If we want to distinguish this from a cosmological
constant at $5\sigma$ (just as CMB experiments seek to detect $r$ for this
model at $5\sigma$) then we require

$\sigma(w_{0})<\frac{1}{5}\times 0.1=0.02\quad\ {\rm Benchmark\ Criterion\
2b.}$ (4)

Combining BC2b and BC1 together gives our absolute benchmark criteria for
investigating the nature of dark energy:

Expansion Benchmark: $\displaystyle\sigma(w_{0})=0.02,\ \sigma(w_{a})=0.05$

In summary, these two parameters can successfully describe dark energy for the
cosmic expansion history, and physics motivates specific benchmark criteria
for the measurement accuracy required for meaningful understanding. Stage 4
accuracy will answer some broadbrush questions about dark energy properties
but Stage 5 is required to meet the detailed benchmark criteria.

##  III The Case for $z\approx 2-5$

The $w_{0}$–$w_{a}$ parametrization of dark energy fits observations with high
accuracy, even when the scalar field behavior is more extreme, such as phase
transitions or dark energy density persisting to high redshift (freezing
fields and early dark energy). We can, however, draw further information about
the evolution by looking to higher redshifts, $z\gtrsim 2$.

For a late time phase transition, we expect dark energy density to vanish at
higher redshifts before the phase transition (e.g. the vacuum metamorphosis
model has its phase transition around $z\approx 1$–2). Conversely, freezing
fields and early dark energy should have a much higher dark energy density at
early times than the cosmological constant. At recombination, $z\approx
10^{3}$, the fraction of the critical density in the cosmological constant is
$\sim 10^{-9}$. However, current constraints on early dark energy density are
at the $10^{-2.5}$ level. Thus there are more than six orders of magnitude of
unexplored range.

While redshifts $z=5$–1000 are fairly inaccessible to precision measurement of
dark energy, even a cosmological constant contributes 8% (1%) of the critical
density at $z=2$ (5), making the dark energy density amenable to constraint.
If we detected no dark energy density at $z=2$, say, this would point us
toward phase transition models, while if we detected greater than 1%
fractional dark energy density at $z=5$ this would favor freezing fields and
early dark energy. A $3\sigma$ distinction between the cosmological constant
behavior and no dark energy density (or double the cosmological constant
density) provides a clear direction for understanding cosmic acceleration
physics. Thus the benchmark criterion is

Density Benchmark: $\displaystyle\sigma\left(\frac{\rho_{\rm de}}{\rho_{\rm
crit}}\right)\equiv\sigma(\Omega_{\rm
de}(z))<\frac{1}{3}\,\Omega_{\Lambda}(z)\,,{\rm for\ }z<5$

##  IV Growth and Gravity

Large scale structure growth can provide incisive constraints on dark energy,
apart from expansion probes. Within general relativity the growth history is
essentially determined by the expansion history, and so $w_{0}$–$w_{a}$ works
as well for cosmic growth. However, deviations of the growth history behavior
relative to the expansion history can reveal extra physics beyond general
relativity. Again, a very successful, single parameter measure of aspects of
growth separate from that predicted by the cosmic expansion was developed in
[13]. The gravitational growth index $\gamma$, given through the growth factor
behavior

$\frac{D(a)}{D(a_{i})}=e^{\int_{a_{i}}^{a}(da/a)\,\Omega_{m}(a)^{\gamma}}\ ,$
(7)

is accurate to 0.2% relative to the exact (subhorizon, linear perturbation,
scale independent) result. For general relativity $\gamma\approx 0.55$ with
only a very mild (and well calibrated) dependence on the dark energy EOS. Thus
$\gamma$ tests the theory of gravity independent of the expansion history.

Establishing a benchmark for $\gamma$ is difficult to do from basic physical
principles, since we have a scarcity of well motivated theories of gravity
apart from general relativity, and many theories can approach general
relativity arbitrarily closely. We do expect that gravity restores to general
relativity at high redshift, and one can define a phase space evolution for
the gravitational coupling $G_{\rm matter}$–$G_{\rm matter}^{\prime}$ [14].
However, it does not separate into distinct regions like the $w$–$w^{\prime}$
phase space. At best we can compare a couple of common models: $f(R)$ scalar-
tensor gravity and DGP braneworld gravity, having $\gamma=0.42$ and
$\gamma=0.68$ respectively. One could then argue that a $5\sigma$ distinction
from general relativity requires $\sigma(\gamma)\approx 0.026$. In fact, due
to the phase space evolution, the $\gamma$ value for $f(R)$ theories is closer
to general relativity’s 0.55 in the past (roughly half the deviation at $z=1$
on larger scales), so for a (admittedly less basic) benchmark criterion we
adopt

$\sigma(\gamma)<\frac{1}{2}\times\frac{1}{5}\times 0.13=0.013\quad\ {\rm
Growth\ Criterion\ 1.}$ (8)

While the gravitational growth index is powerful for a single parameter,
giving a critical alert of growth deviating from general relativity, there are
two other important avenues providing distinction from general relativity.
Growth in the linear or quasilinear perturbation regime can be scale dependent
in modified gravity, and light propagation as well as cosmic growth can
deviate. These can be described (in the subhorizon, quasistatic regime) as
gravitational couplings entering modified Poisson equations for matter and for
light [15] in place of Newton’s constant. Thus they are often called $G_{\rm
matter}$ and $G_{\rm light}$, normalized to unity when they equal Newton’s
constant (there are many other names also, e.g. $\mu$ and $\Sigma$, but
$G_{\rm matter}$ and $G_{\rm light}$ are the most easily understood).

We can explore deviations from general relativity in a model independent
manner (i.e. without adopting a specific theory) by treating the time and
scale dependent $G_{\rm matter}(a,k)$ and $G_{\rm light}(a,k)$ as values
defined in bins of scale factor $a$ and density perturbation wavenumber $k$.
Two bins in $k$ and three bins in $a$, so six parameters total, have been
demonstrated to be accurate to 0.1%–0.3% rms in the growth rate of large scale
structure (better on the growth amplitude itself) [16]. In many cases only two
bins in $a$ are needed for 0.2% accuracy, giving four parameters.

As with $\gamma$ it is difficult to develop from first principles a benchmark
for required accuracy, as gravity theories can approach general relativity
arbitrarily closely. Again looking to the $f(R)$ gravity model, we can seek a
clear signature of deviation from general relativity in scale and time
dependence of gravitational coupling, adopting

$\displaystyle{\rm Growth\ Criterion\ 2.}$ $\displaystyle\sigma(G_{{\rm
matter,low\ }k,{\rm high\ }a})$ $\displaystyle<$ $\displaystyle 0.02$ (9)
$\displaystyle\sigma(G_{{\rm matter,high\ }k,{\rm high\ }a})$ $\displaystyle<$
$\displaystyle 0.05\ ,$ (10)

which each correspond to roughly $3\sigma$ constraints for a canonical $f(R)$
model. Here low (high) $k$ corresponds to $k=0.055$ ($0.125$) $h$/Mpc and high
$a$ is $a\approx 0.75$ ($z\approx 0.35$). One can use similar criteria on
$G_{\rm light}$.

Another important aspect of modified gravity becomes apparent as one
approaches the nonlinear perturbation regime: a screening mechanism to restore
the theory to general relativity behavior in regions of high density or
density gradients (to satisfy solar system tests for example). Here the
physics becomes complex and model dependent enough that it seems worthwhile
abandoning parametrization and dealing with a full theory, starting from its
Lagrangian. There are in-between methods such as effective field theory, but
these do not capture the fully nonlinear aspects. At this point we have to
choose some benchmark models, rather than benchmark principles. This is
difficult as few to no models are regarded as completely compelling by the
community as a whole.

We therefore consider two illustrative models, with the aim of keeping them as
simple as possible while showcasing essential elements of physics. One is
$f(R)$ gravity, but it is important to recognize that stating this by itself
is insufficient: one must adopt a particular functional form for $f(R)$, with
particular parameter values. A common choice is the Hu-Sawicki form [17],
which has two parameters (apart from the density) – an amplitude $f_{R0}$ and
an evolution scaling (power law index $n$). We instead adopt for simplicity a
single parameter exponential $f(R)$ model,

$f(R)=-b\left(1-e^{-cR/b}\right)\ ,$ (11)

where $b$ is determined by the present matter density $\Omega_{m,0}$ and $c$
is the one parameter [18]. Thus our first full gravitational theory benchmark
is

$\displaystyle{\rm Benchmark\ Gravity\ 1.}$ $\displaystyle{\rm Exponential\
}f(R)\ {\rm with\ }c$ $\displaystyle=$ $\displaystyle 4$ (12)

Note that neither exponential nor Hu-Sawicki $f(R)$ have a strong theoretical
foundation for their functional form. The benchmark is simply a point in the
theory landscape that is sufficiently distinct from general relativity to be
interesting and currently observationally viable.

One of the drawbacks of using $f(R)$ gravity is that it has essentially no
effect on light propagation different from general relativity. Therefore we
adopt a second gravity benchmark theory that does affect light propagation.

It is useful to briefly consider the effective field theory (EFT) or property
function approach [19]. The two key property functions are $\alpha_{M}$,
giving the running of the Planck mass arising from coupling of the scalar
field to the Ricci scalar, and $\alpha_{B}$, describing the braiding between
the scalar field kinetic structure and the metric. All $f(R)$ gravity theories
have $\alpha_{B}=-\alpha_{M}$, which is also the condition for negligible
effect on light propagation (there is a residual effect of order $f_{R0}$, but
this is usually $\lesssim 10^{-5}$). This arises from the conformal nature of
the theory (so null geodesics are unaffected), and the coupling to the Ricci
scalar can give rise to the chameleon mechanism for screening.

Thus, as a second gravity benchmark theory we will choose a theory relying on
a different mechanism, called Vainshtein screening, which can arise due to the
kinetic structure, or braiding. Since we also want to have a minimal number of
free parameters we could simply change the constant of proportionality
relating $\alpha_{B}$ to $\alpha_{M}$ – for example No Slip Gravity has
$\alpha_{B}=-2\alpha_{M}$ (and then $G_{\rm matter}=G_{\rm light}$). But to
emphasize the distinction, and maintain as few parameters as possible, we will
choose a class of theories with standard gravitational coupling to the Ricci
scalar, so a constant Planck mass (i.e. $8\pi G_{N}=m_{p}^{2}$) hence
$\alpha_{M}=0$. Similarly, the benchmark will have a canonical kinetic term
(though possibly of opposite sign) with no potential, and a cubic Horndeski
term only dependent on the field kinetic energy $X$, i.e. $G_{3}=G_{3}(X)$, so
the action is shift symmetric.

Kinetic gravity braiding [20, 21] is such a class of theories. The specific
form $G_{3}(X)$ to choose, keeping in mind minimal number of parameters, is
unclear. We know a form $G_{3}=bX^{n}$, say, requires $n\geq 0$ for stability.
The most natural choices might be $n=0$ – but this gives $\alpha_{B}=0$ and no
light propagation effect; $n=1/2$, i.e. $G_{3}\sim\dot{\phi}$ – but this can
have stability issues; $n=1$, i.e. $G_{3}\sim X$ – but this can have
observational issues with a wrong sign Sachs-Wolfe effect. Also, $G_{3}\sim\ln
X$ does not work as it does not restore to general relativity at high
redshift. One could instead work with $\alpha_{B}$, and take a simple form for
it such as $\alpha_{B}(a)=A/[1+(4a)^{-3}]$, which has the physically expected
behaviors of vanishing at high redshift (restoring general relativity) and
freezing in the de Sitter future. However, property functions come from a
linear level EFT and do not fully describe the nonlinear theory (e.g. the
screening). Thus we leave for future work the determination of a one
parameter, fully nonlinear theory in the manner that $f(R)$ is, that has an
effect on light propagation and possibly minimal coupling to the Ricci scalar.
The second full gravitational theory benchmark has a placeholder

$\displaystyle{\rm Benchmark\ Gravity\ 2.}$ (13) $\displaystyle{\rm Kinetic\
Gravity\ Braiding\ with\ 1\ parameter?}$

##  V Achieving Benchmarks

As mentioned, Stage 4 experiments underway this decade will make significant
inroads toward quantifying dark energy. To achieve the physics derived or
motivated benchmarks described here will require Stage 5 experiments however.
While in part such experiments can scale up from Stage 4 techniques, some
different approaches are necessary as well.

For the Density Benchmark, surveys must extend to $z\approx 5$. A well tested
approach is that of spectroscopic galaxy clustering, and plans for a
Spectroscopic Stage 5 experiment are in design [22, 23]. This, as well as
lower redshift techniques such as peculiar velocity surveys, high density
surveys, and spectrophotometric supernova surveys [24, 25, 26] will also
improve the Expansion, Growth, and Gravity Benchmarks.

Significant progress on the challenging benchmark criteria for the Expansion
Benchmarks (the most motivated from physical first principles) likely will
rely on more than scaling up existing probes. A novel technique coming toward
maturity is direct measurement of acceleration through redshift drift. While
other probes involve one or two integrals over the dark energy equation of
state, redshift drift $dz/dt_{\rm obs}=H_{0}\,(1+z)-H(z)$ is related directly
to the acceleration $\ddot{a}$ and hence to $w(a)$.

Two major advances in the last decade make it viable. On the theoretical front
it was shown that low redshift observations ($z\lesssim 0.5$) of the well
characterized class of emission line galaxies, using differential motion of
bright forbidden transition doublet lines, was optimal. Furthermore, in
combination with CMB observations redshift drift could achieve the Expansion
Benchmark with five measurements of 1% accuracy from $z=0.1$–0.5 [27]. On the
experimental front, interferometers in tandem with spectrographs have now used
differential measurements to demonstrate control of systematics $1000\times$
better than conventional spectrographs [28, 29, 30]. Exoplanet searches are
driving similar requirements and the technique has been demonstrated on the
Keck Planet Finder. Improvements are being tested currently using Fabry-Perot
interferometry to further reduce systematics. Most impressively, the redshift
drift probe simply adds onto existing infrastructure: a suitcase size
interferometer is used in conjunction with an existing spectrograph at a 10
meter telescope (though larger telescopes give further S/N improvement).

Between a cosmic redshift drift experiment and extensions of existing probes,
a Stage 5 program can achieve the necessary benchmarks to answer fundamental
physics questions about dark energy.

##  VI Conclusions

While we do not have clear clues as to the physical origin of cosmic
acceleration and dark energy, we do know the questions to ask to obtain
insights:

  1. 1.

There are two natural e-fold timescales for dark energy dynamics: 0 and
$\mathcal{O}(1)$ – can we distinguish these observationally at $5\sigma$?

  2. 2.

There are two natural phase space regions for dark energy dynamics: freezing
fields moving toward cosmological constant behavior and thawing fields moving
away from it – can we distinguish these observationally at $2\sigma$
(comparing mean of one region to another)?

  3. 3.

If cosmic acceleration from inflation is connected to late time acceleration
through an $\alpha$-attractor origin, can we detect its deviation from a
cosmological constant at $5\sigma$, matching the $5\sigma$ detection of
primordial gravitational waves in the CMB for $\alpha=1$ Starobinsky/Higgs
inflation?

For further dark energy dynamics we must look beyond the region of its
domination but where it is still observationally relevant, at $z\approx 2$–5.
Since dark energy is allowed to either vanish at high redshift, or contribute
more than a million times more energy density than the cosmological constant
by recombination, we should measure the dark energy density to allow at least
a $3\sigma$ determination of whether dark energy dies off or significantly
exceeds $\rho_{\Lambda}$ out to $z\approx 5$.

Beyond the Expansion Benchmarks and Density Benchmark, we must test whether
the cosmic growth history is governed by the cosmic expansion history in the
manner dictated by general relativity. In the first instance this can be
accomplished by separating expansion from growth through the gravitational
growth index $\gamma$. The next key element of understanding comes from the
time and space dependence of the gravitational coupling, tested both for
matter and for light. Finally, to access the full predictions of a gravity
theory into the nonlinear regime we adopt two specific benchmark models of
gravity, using only a single parameter each but with a rich array of
observational implications.

All these goals are achievable, starting in the next decade and within the 20
year vision. The Stage 5 program to answer these fundamental physics questions
on the nature of dark energy merges extensions of current Stage 4 probes plus
the highly complementary and novel probe of cosmic redshift drift. The
redshift drift technique is rapidly maturing, with strong advances in
systematics control (due in part to synergy with exoplanet search goals) and
survey design.

High energy physics can advance significantly toward the “vitally important”,
“most profound”, and “comprehensive fundamental physical” nature of dark
energy.

###### Acknowledgements.

I thank Hitoshi Murayama for interesting discussions about dark energy
benchmarks. This work was supported in part by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics, under contract no. DE-
AC02-05CH11231.

## References

  * [1] These quotes are from Steven Weinberg, Frank Wilczek, and Edward Witten, respectively, in Resource Book on Dark Energy, ed. E. Linder, Snowmass 2001 https://supernova.lbl.gov/~evlinder/sci.html
  * [2] E.V. Linder, The Mirage of $\Lambda$, arXiv:0708.0024
  * [3] E.V. Linder, Exploring the Expansion History of the Universe, Phys. Rev. Lett. 90, 091301 (2003) [arXiv:astro-ph/0208512]
  * [4] L. Parker, A. Raval, New Quantum Aspects of a Vacuum-Dominated Universe, Phys. Rev. D 62, 083503 (2000) [arXiv:gr-qc/0003103]
  * [5] R.R. Caldwell, W. Komp, L. Parker, D.A.T. Vanzella, A Sudden Gravitational Transition, Phys. Rev. D 73, 023513 (2006) [arXiv:astro-ph/0507622]
  * [6] E. Di Valentino. E.V. Linder, A. Melchiorri, A Vacuum Phase Transition Solves $H_{0}$ Tension, Phys. Rev. D 97, 043528 (2018) [arXiv:1710.02153]
  * [7] A. Albrecht, C. Skordis, Phenomenology of a Realistic Accelerating Universe Using Only Planck-Scale Physics, Phys. Rev. Lett. 84, 2076 (2000) [arXiv:astro-ph/9908085]
  * [8] R. de Putter, E.V. Linder, Calibrating Dark Energy, JCAP 0810, 042 (2008) [arXiv:0808.0189]
  * [9] R.R. Caldwell, E.V. Linder, The Limits of Quintessence, Phys. Rev. Lett. 95, 141301 (2005) [arXiv:astro-ph/0505494]
  * [10] R.N. Cahn, R. de Putter, E.V. Linder, Field Flows of Dark Energy, JCAP 0811, 015 (2008) [arXiv:0807.1346]
  * [11] E.V. Linder, The Paths of Quintessence, Phys. Rev. D 73, 063010 (2006) [arXiv:astro-ph/0601052]
  * [12] Y. Akrami, R. Kallosh, A. Linde, V. Vardanyan, Dark Energy, $\alpha$-attractors, and Large-Scale Structure Surveys, JCAP 1806, 041, 2018 [arXiv:1712.09693]
  * [13] E.V. Linder, Cosmic Growth History and Expansion History, Phys. Rev. D 72, 043529 (2005) [arXiv:astro-ph/0507263]
  * [14] E. V. Linder, Model-Independent Tests of Cosmic Gravity, Phil. Trans. Roy. Soc. A 369, 4985 (2011) [arXiv:1103.0282]
  * [15] E. Bertschinger, P. Zukin, Distinguishing Modified Gravity from Dark Energy, Phys. Rev. D 78, 024015 (2008) [arXiv:0801.2431]
  * [16] M. Denissenya, E.V. Linder, Constraining Scale Dependent Growth with Redshift Surveys, JCAP 2211, 029 (2022) [arXiv:2208.10508]
  * [17] W. Hu, I. Sawicki, Models of $f(R)$ Cosmic Acceleration that Evade Solar-System Tests, Phys. Rev. D 76, 064004 (2007) [arXiv:0705.1158]
  * [18] E. V. Linder, Exponential Gravity, Phys. Rev. D 80, 123528 (2009) [arXiv:0905.2962]
  * [19] E. Bellini, I. Sawicki, Maximal Freedom at Minimum Cost: Linear Large-Scale Structure in General Modifications of Gravity, JCAP 1407, 050 (2014) [arXiv:1404.3713]
  * [20] C. Deffayet, O. Pujolas, I. Sawicki, A. Vikman, Imperfect Dark Energy from Kinetic Gravity Braiding, JCAP 1010, 026 (2010) [arXiv:1008.0048]
  * [21] O. Pujolas, I. Sawicki, A. Vikman, The Imperfect Fluid behind Kinetic Gravity Braiding, JHEP 1111, 156 (2011) [arXiv:1103.5360]
  * [22] D.J. Schlegel et al., A Spectroscopic Road Map for Cosmic Frontier: DESI, DESI-II, Stage-5, arXiv:2209.03585
  * [23] D.J. Schlegel et al., MegaMapper: A Stage-5 Spectroscopic Instrument Concept for the Study of Inflation and Dark Energy, arXiv:2209.04322
  * [24] A. Palmese, A.G. Kim, Probing Gravity and Growth of Structure with Gravitational Waves and Galaxies’ Peculiar Velocity, Phys. Rev. D 103, 103507 (2021) [arXiv:2005.04325]
  * [25] K. Dawson, A. Hearin, K. Heitmann, M. Ishak, J.U. Lange, M. White, R. Zhou, Snowmass2021 Cosmic Frontier White Paper: High Density Galaxy Clustering in the Regime of Cosmic Acceleration, arXiv:2203.07291
  * [26] K. Boone et al., Twins Embedding of Type Ia Supernovae II: Improving Cosmological Distance Estimates, ApJ 912, 71 (2021) [arXiv:2105.02204]
  * [27] A.G. Kim, E.V. Linder, J. Edelstein, D. Erskine, Astropart. Phys. 62, 195 (2015) [arXiv:1402.6614]
  * [28] D. Erskine, E. Wishnow, E. Linder, E.J. Davies, M. Sirk, R. Ozer, D. Fratanduono, J. Edelstein, Spectrograph stabilization of 500x using a single-delay interferometer, SPIE 12184 (2022) https://doi.org/10.1117/12.2628388
  * [29] D. Erskine, J. Edelstein, E. Wishnow, E.J. Davies, M. Sirk, R. Ozer, D. Fratanduono, Externally dispersed interferometer testbed diagnosing Keck Planet Finder spectrograph high resolution performance, SPIE 12184 (2022) https://doi.org/10.1117/12.2627160
  * [30] D.J. Erskine, Method for boosting dispersive spectrograph stability 1000x using interferometry with crossfaded pairs of delays, JATIS 7, 025006 (2021) https://doi.org/10.1117/1.JATIS.7.2.025006


Numerical confirmations of joint spike transitions in $G_{2}$ cosmologies

W C Lim

Department of Mathematics, University of Waikato, Private Bag 3105, Hamilton
3240, New Zealand

<EMAIL_ADDRESS>
###### Abstract

We produce numerical evidence that the joint spike transitions between Kasner
eras of $G_{2}$ cosmologies are described by the non-orthogonally transitive
$G_{2}$ spike solution. A new matching procedure is developed for this
purpose.

Keywords: Spike, joint spike transition, numerical matching, area time gauge,
zooming technique

##  1 Introduction

Belinskii, Khalatnikov and Lifschitz [art:LK63, art:BKL1970, art:BKL1982]
describe the dynamics of a spacetime (with fluids with soft equation of state)
as it evolves toward a generic spacelike singularity as an infinite sequence
of transitions between anisotropic vacuum Kasner saddle states. The
transitions are described by vacuum Bianchi type II solutions. The time spent
near a Kasner state is called a Kasner epoch. Each Kasner state is
characterised by its BKL parameter $u\geq 1$, which decreases by 1 for each
Bianchi type II transition. A Kasner era consists of a sequence of Kasner
epochs with decreasing BKL parameter $u$. If $u<2$ before the transition, then
$u$ is mapped to $\frac{1}{u-1}$ after the transition, starting a new Kasner
era.

Due to spatial inhomogeneities, neighbouring worldlines experience slightly
different, but increasingly diverging, sequences. The dynamics are
asymptotically local, as if each typical worldline follows the dynamics of a
spatially homogeneous model. Exception to this locality was found when spikes
were discovered numerically by Berger and Moncrief [art:BergerMoncrief1993] in
the context of vacuum models with two commuting Killing vector fields that act
orthogonally transitively (OT $G_{2}$ models). Spike forms along worldlines
where the crucial variable responsible for the Bianchi type II transition
($N_{-}$ in our formulation below) is zero, as it changes sign due to spatial
inhomogeneities [art:ColeyLim2014]. When a spike forms, its width can become
temporarily narrower than the particle horizon, and its dynamics is described
by a spatially inhomogeneous solution. The exact OT $G_{2}$ spike solution was
found in 2008 [art:Lim2008], and is shown to match non-moving spikes in
numerical simulations in 2009 [art:Limetal2009], which introduces a zooming
technique as a cost-effective way to maintain adequate numerical resolution.
The exact solution describes a transition between Kasner states that are
connected by two consecutive Bianchi type II transitions. OT $G_{2}$ models
allow only a single Kasner era plus the first Kasner epoch of the next era, so
a spike forming near the end of the era ends up as a permanent spike. Relaxing
the orthogonal transitivity condition allows multiple Kasner eras and a non-
terminating sequence of transitions. The permanent spike is replaced by the
so-called joint spike transition [art:HeinzleUgglaLim2012], which is a more
elaborate spike transition that straddles two Kasner eras. The results of
[art:Limetal2009] is therefore incomplete until a new exact solution
describing the joint spike transition is found and numerically matched.

The exact non-OT $G_{2}$ spike solution that seems to describe such a
transition was found in 2015 [art:Lim2015]. However, the exact non-OT $G_{2}$
spike solution uses a different time parameterisation, and it was not known at
the time how to recover this time variable in numerical simulations. In
[art:ColeyLim2016], where the exact solution is generalised to the stiff fluid
case, the fluid-comoving or volume gauge seemed to be a more natural gauge
choice than the area time gauge. The numerical simulation in
[art:ColeyLim2016] uses the fluid-comoving gauge in conjunction with a dynamic
version of the zooming technique (developed in [art:LimRegisClarkson2013]).
Doing so comes at a cost of fine-tuning the excision boundary. The precision
of fine-tuning is determined by the size of the spatial domain at the end of
simulation relative to that at the start. Fine-tuning requires one to run
numerical simulations multiple times. Therefore, while the goal to extend the
result of the paper [art:Limetal2009] can be achieved, it would come at a
great cost, which is very inefficient and unsatisfactory. Since 2017, much
effort has been spent on improving the fine-tuning through better
understanding of the transition times of exact solution. A method of analysing
the transition times was developed in the doctoral thesis of Moughal
[thesis:Moughal2021, art:MoughalLim2021], and it was then applied to the non-
OT $G_{2}$ spike solution in [art:LimMoughal2022]. Despite the better
understanding, improvement in the fine-tuning is very little. Attempts to
modify the zooming technique to avoid the fine-tuning problem have been
unsuccessful.

The breakthrough came in 2021 when it was realised that we can actually keep
using the area time gauge and _evolve the time variable of the exact solution_
to accommodate numerical matching, thus avoiding fine-tuning. This paper will
show how this is done, develop a new matching procedure, and show that the
exact solution indeed matches the joint spike transitions. The reader should
read the paper [art:Limetal2009] and make frequent comparisons as the two
papers are similar in the approach.

##  2 $G_{2}$ spacetimes

We use the orthonormal frame approach [book:WainwrightEllis1997,
art:Ugglaetal2003] to formulate the Einstein field equations, adopting the
Iwasawa spatial frame [art:HeinzleUgglaRohr2009] and $\beta$-normalised
variables [art:vEUW2002] for numerical evolution (although we will plot some
Hubble-normalised variables). We represent the metric components of non-OT
$G_{2}$ models the same way as in [art:Lim2015], where indices $0$, $1$, $2$,
$3$ correspond to coordinates $t$, $x$, $y$, $z$, and the metric components
are given in terms of $N$, $b^{1}$, $b^{2}$, $b^{3}$, $n_{1}$, $n_{2}$,
$n_{3}$ as follows:

$\displaystyle g_{00}$ $\displaystyle=-N^{2},$ (1) $\displaystyle g_{11}$
$\displaystyle=\text{e}^{-2b^{1}},\quad g_{12}=\text{e}^{-2b^{1}}n_{1},\quad
g_{13}=\text{e}^{-2b^{1}}n_{2},$ (2) $\displaystyle g_{22}$
$\displaystyle=\text{e}^{-2b^{2}}+\text{e}^{-2b^{1}}n_{1}{}^{2},\quad
g_{23}=\text{e}^{-2b^{1}}n_{1}n_{2}+\text{e}^{-2b^{2}}n_{3},$ (3)
$\displaystyle g_{33}$
$\displaystyle=\text{e}^{-2b^{3}}+\text{e}^{-2b^{1}}n_{2}{}^{2}+\text{e}^{-2b^{2}}n_{3}{}^{2}.$
(4)

The metric components depend on $t$ and $z$ only. Note the change in alignment
and notation from $(\tau,x)$ in [art:Limetal2009] to $(t,z)$ here.

The expansion shear components $\Sigma_{\alpha\beta}$ and remaining nonzero
spatial curvature components $N_{\alpha\beta}$ here are decomposed as follows:

$\displaystyle\Sigma_{\alpha\beta}$
$\displaystyle=\left(\begin{matrix}\Sigma_{+}+\sqrt{3}\Sigma_{-}&\sqrt{3}\Sigma_{\times}&\sqrt{3}\Sigma_{2}\\\
\sqrt{3}\Sigma_{\times}&\Sigma_{+}-\sqrt{3}\Sigma_{-}&\sqrt{3}\Sigma_{1}\\\
\sqrt{3}\Sigma_{2}&\sqrt{3}\Sigma_{1}&-2\Sigma_{+}\end{matrix}\right),$ (5)
$\displaystyle N_{\alpha\beta}$
$\displaystyle=\left(\begin{matrix}2\sqrt{3}N_{-}&\sqrt{3}N_{\times}&0\\\
\sqrt{3}N_{\times}&0&0\\\ 0&0&0\end{matrix}\right).$ (6)

We use the frame rotation freedom to set $\Sigma_{2}=0$.

The $\beta$-normalised variables in terms of the metric components are:

$\displaystyle\beta=-\frac{1}{2}N^{-1}\partial_{t}(b^{1}+b^{2}),\quad\mathcal{N}=N\beta,\quad
E_{3}{}^{3}=\frac{\text{e}^{b^{3}}}{\beta},$ (7)
$\displaystyle\Sigma_{+}=\frac{1}{3}(1+\mathcal{N}^{-1}\partial_{t}b^{3}),\quad\Sigma_{1}=\frac{\text{e}^{b^{2}-b^{3}}}{2\sqrt{3}}\mathcal{N}^{-1}\partial_{t}n_{3},$
(8)
$\displaystyle\Sigma_{-}=-\frac{1}{2\sqrt{3}}\mathcal{N}^{-1}\partial_{t}(b^{1}-b^{2}),\quad
N_{\times}=\frac{1}{2\sqrt{3}}E_{3}{}^{3}\partial_{z}(b^{1}-b^{2}),$ (9)
$\displaystyle\Sigma_{\times}=\frac{\text{e}^{b^{1}-b^{2}}}{2\sqrt{3}}\mathcal{N}^{-1}\partial_{t}n_{1},\quad
N_{-}=\frac{\text{e}^{b^{2}-b^{1}}}{2\sqrt{3}}E_{3}{}^{3}\partial_{z}n_{1},$
(10)

where $\partial_{t}$ and $\partial_{z}$ denote partial differentiation with
respect to $t$ and $z$ respectively. $t$ tends to infinity toward the
singularity. The Hubble expansion scalar $H$ is related to $\beta$ through
$H=\beta(1-\Sigma_{+})$, so Hubble-normalised variables (denoted with a
superscript $H$) are related to $\beta$-normalised variables through

$\Sigma_{+}^{H}=\frac{\Sigma_{+}}{1-\Sigma_{+}},\quad\Sigma_{-}^{H}=\frac{\Sigma_{-}}{1-\Sigma_{+}},$
(11)

and so on.

We use the same temporal gauge and time parameterisation for numerical
evolution as in [art:Limetal2009], namely the area time gauge [art:vEUW2002]
and the time parameterisation such that

$\mathcal{N}=-\frac{1}{2}.$ (12)

To maintain adequate numerical resolution as spikes become narrow when they
form, we use the same zooming technique as in [art:Limetal2009], introducing
zooming coordinates $(T,Z)$ to zoom in on a specified worldline
$z=z_{\text{zoom}}$:

$T=t,\quad Z=\frac{2}{(E_{3}^{3})_{0}}\text{e}^{t}(z-z_{\text{zoom}}).$ (13)

Compare with Equation (20) of [art:Limetal2009]. In this paper we allow
$(E_{3}{}^{3})_{0}$ (the initial value of $E_{3}{}^{3}$) to take value other
than $2$. We fix the zoom rate to the natural zoom rate ($A=1$ in Equation
(20) of [art:Limetal2009]) here because we can now better estimate the growth
of the lower bound on the right excision boundary, which we will elaborate on
below. The differential operators in the new coordinates are

$\partial_{t}=\partial_{T}+Z\partial_{Z},\quad\partial_{z}=\frac{2}{(E_{3}^{3})_{0}}\text{e}^{T}\partial_{Z}.$
(14)

The evolution equations in the zooming coordinates are

$\displaystyle\partial_{T}E_{3}{}^{3}$
$\displaystyle=-Z\partial_{Z}E_{3}{}^{3}-(1-\tfrac{3}{2}\Sigma_{1}^{2})E_{3}{}^{3}$
(15) $\displaystyle\partial_{T}\Sigma_{-}$
$\displaystyle=-Z\partial_{Z}\Sigma_{-}+\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}\partial_{Z}N_{\times}+\frac{3}{2}\Sigma_{1}^{2}\Sigma_{-}-\sqrt{3}(\Sigma_{\times}^{2}-N_{-}^{2})+\frac{\sqrt{3}}{2}\Sigma_{1}^{2}$
(16) $\displaystyle\partial_{T}N_{\times}$
$\displaystyle=-Z\partial_{Z}N_{\times}+\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}\partial_{Z}\Sigma_{-}-(1-\tfrac{3}{2}\Sigma_{1}^{2})N_{\times}$
(17) $\displaystyle\partial_{T}\Sigma_{\times}$
$\displaystyle=-Z\partial_{Z}\Sigma_{\times}-\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}\partial_{Z}N_{-}+(\tfrac{3}{2}\Sigma_{1}^{2}+\sqrt{3}\Sigma_{-})\Sigma_{\times}+\sqrt{3}N_{\times}N_{-}$
(18) $\displaystyle\partial_{T}N_{-}$
$\displaystyle=-Z\partial_{Z}N_{-}-\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}\partial_{Z}\Sigma_{\times}-(1-\tfrac{3}{2}\Sigma_{1}^{2}-\sqrt{3}\Sigma_{-})N_{-}-\sqrt{3}N_{\times}\Sigma_{\times}$
(19) $\displaystyle\partial_{T}\Sigma_{1}$
$\displaystyle=-Z\partial_{Z}\Sigma_{1}+\tfrac{1}{2}(3\Sigma_{1}^{2}+3\Sigma_{+}-\sqrt{3}\Sigma_{-})\Sigma_{1}$
(20) where $\displaystyle\Sigma_{+}$
$\displaystyle=\tfrac{1}{2}(1-\Sigma_{-}^{2}-\Sigma_{\times}^{2}-\Sigma_{1}^{2}-N_{-}^{2}-N_{\times}^{2}).$
(21)

For numerical accuracy, we choose to evolve the logarithm of $E_{3}{}^{3}$ and
$\Sigma_{1}$. There is one constraint equation:

$\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}\partial_{Z}\Sigma_{1}=(3N_{-}\Sigma_{\times}-3N_{\times}\Sigma_{-}-\sqrt{3}N_{\times})\Sigma_{1}.$
(22)

Equations (15)–(22) are essentially the same as Equations (22)–(28) of
[art:Limetal2009].

The characteristic velocity for the evolution equations (15) and (20) is $Z$.
To find the characteristic velocities for the subsystem (16)–(19), write them
in the form

$\displaystyle\partial_{T}(\Sigma_{-}+N_{\times})$
$\displaystyle=-(Z-\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}})\partial_{Z}(\Sigma_{-}+N_{\times})+\cdots,$
(23) $\displaystyle\partial_{T}(\Sigma_{-}-N_{\times})$
$\displaystyle=-(Z+\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}})\partial_{Z}(\Sigma_{-}-N_{\times})+\cdots,$
(24) $\displaystyle\partial_{T}(\Sigma_{\times}-N_{-})$
$\displaystyle=-(Z-\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}})\partial_{Z}(\Sigma_{\times}-N_{-})+\cdots,$
(25) $\displaystyle\partial_{T}(\Sigma_{\times}+N_{-})$
$\displaystyle=-(Z+\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}})\partial_{Z}(\Sigma_{\times}+N_{-})+\cdots.$
(26)

Then we see that their characteristic velocities are
$Z\pm\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}$. So the maximum and
minimum characteristic velocities of the system (15)–(20) are

$v_{\text{max}}=Z+\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}},\quad
v_{\text{min}}=Z-\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}}.$ (27)

To avoid specifying boundary conditions, we want to ensure that all
characteristic velocities are outgoing at the excision boundary. So we want to
place the excision boundary far away enough from $Z=0$, so that
$v_{\text{min}}\geq 0$ at the right excision boundary $Z=Z_{r}$ and
$v_{\text{max}}\leq 0$ at the left excision boundary $Z=Z_{l}$ throughout the
duration of numerical evolution. How does $E_{3}{}^{3}$ behave? The joint
spike transition undergoes one $\Sigma_{1}$ frame transition, during which
$\Sigma_{1}$ is of order 1. $\Sigma_{1}$ is otherwise negligible. From
numerical observations, we see that

$E_{3}{}^{3}\approx(E_{3}{}^{3})_{0}\text{e}^{-T}$ (28)

before the $\Sigma_{1}$ transition, and

$E_{3}{}^{3}\approx\frac{4}{(w-1)^{2}}(E_{3}{}^{3})_{0}\text{e}^{-T}$ (29)

after the $\Sigma_{1}$ transition (for the case $0<w<1$, as we will make this
choice in (55) below), where $w$ is a parameter of the exact spike solution
below. This means $v_{\text{min}}$ drops from $Z-1$ to
$Z-\frac{4}{(w-1)^{2}}$. That is, the lower bound for the right-hand excision
boundary, obtained by solving $v_{\text{min}}=0$ for $Z$ in (27):

$Z_{\text{right bound}}=\text{e}^{T}\frac{E_{3}{}^{3}}{(E_{3}{}^{3})_{0}},$
(30)

grows from $Z_{\text{right bound}}=1$ to $Z_{\text{right
bound}}=\frac{4}{(w-1)^{2}}$. Therefore the right excision boundary $Z=Z_{r}$
must satisfy $Z_{r}\geq\frac{4}{(w-1)^{2}}$. Similarly, the left excision
boundary $Z=Z_{l}$ must satisfy $Z_{l}\leq-\frac{4}{(w-1)^{2}}$.

In the matching procedure below, the value of $w$ can be obtained before the
$\Sigma_{1}$ transition occurs, thus the excision boundary can be easily
adjusted after a first run. If $w$ is too close to $1$, then $Z_{\text{right
bound}}$ is too large and the numerical simulation becomes too expensive to
run, and one really needs a more sophisticated zooming technique. If the
numerical run covers a second $\Sigma_{1}$ transition into the third Kasner
era, then $Z_{\text{right bound}}$ will increase a second time, and again a
more sophisticated zooming technique should be used to reduce wastage of
numerical resources. We shall leave that to future research.

The Courant-Friedrichs-Lax condition for numerical stability requires that the
numerical timestep size $\Delta T$ satisfies

$\Delta T<\frac{\Delta Z}{v_{\text{max}}},$ (31)

evaluated at the excision boundary, where $\Delta Z$ is the numerical grid
size. We use

$\Delta T=0.9\frac{\Delta Z}{v_{\text{max}}}.$ (32)

We use the same numerical method as in [art:Limetal2009], namely the classical
fourth-order Runge-Kutta method, with fourth-order accurate spatial
derivatives. We use double precision for our numerical runs, and use quad
precision to check that our double-precision runs are accurate.

##  3 The exact non-OT $G_{2}$ spike solution

The exact non-OT $G_{2}$ spike solution in [art:Lim2015] with $K=0$ and
$\omega_{0}=0$ is given by:

$\displaystyle N$
$\displaystyle=-\text{e}^{-\frac{1}{4}(w^{2}+3)\tau}\sqrt{\omega^{2}+\lambda^{2}}$
(33) $\displaystyle\text{e}^{-2b^{1}}$
$\displaystyle=\frac{\lambda}{\omega^{2}+\lambda^{2}}$ (34)
$\displaystyle\text{e}^{-2b^{2}}$
$\displaystyle=\frac{\mathcal{A}^{2}}{\lambda}(\omega^{2}+\lambda^{2})$ (35)
$\displaystyle\text{e}^{-2b^{3}}$
$\displaystyle=\text{e}^{-\frac{1}{2}(w^{2}+3)\tau}\mathcal{A}^{-2}(\omega^{2}+\lambda^{2})$
(36) $\displaystyle n_{1}$
$\displaystyle=-2w(w-1)n_{30}z^{2}+\frac{\omega^{2}}{\lambda}(n_{30}\text{e}^{-(w+1)\tau}+n_{10}n_{20}\text{e}^{-\frac{1}{2}(w^{2}-1)\tau})$
$\displaystyle\quad-\Bigg{[}n_{30}w\text{e}^{-2\tau}+n_{10}n_{20}\frac{(w+3)}{(w-1)}\text{e}^{-\frac{1}{2}(w-1)^{2}\tau}$
$\displaystyle\quad\qquad+n_{20}n_{30}(n_{10}n_{30}-n_{20})\frac{(w-3)}{(w+1)}\text{e}^{-\frac{1}{2}(w+1)^{2}\tau}\Bigg{]}$
(37) $\displaystyle n_{2}$
$\displaystyle=n_{20}\Bigg{[}\frac{\omega^{2}}{\lambda}\text{e}^{-\frac{1}{2}(w^{2}-1)\tau}-\frac{(w+3)}{(w-1)}\text{e}^{-\frac{1}{2}(w-1)^{2}\tau}-n_{30}^{2}\frac{(w-3)}{(w+1)}\text{e}^{-\frac{1}{2}(w+1)^{2}\tau}\Bigg{]}$
(38) $\displaystyle n_{3}$
$\displaystyle=\mathcal{A}^{-2}\left[n_{10}\text{e}^{-\frac{1}{2}(w-1)^{2}\tau}+n_{30}(n_{10}n_{30}-n_{20})\text{e}^{-\frac{1}{2}(w+1)^{2}\tau}\right],$
(39) where $\displaystyle\mathcal{A}^{2}$
$\displaystyle=\text{e}^{-2\tau}+n_{10}^{2}\text{e}^{-\frac{1}{2}(w-1)^{2}\tau}+(n_{10}n_{30}-n_{20})^{2}\text{e}^{-\frac{1}{2}(w+1)^{2}\tau}$
(40) $\displaystyle\lambda$
$\displaystyle=\text{e}^{(w-1)\tau}+n_{20}^{2}\text{e}^{-\frac{1}{2}(w^{2}-1)\tau}+n_{30}^{2}\text{e}^{-(w+1)\tau}$
(41) $\displaystyle\omega$ $\displaystyle=2wn_{30}z.$ (42)

To keep $\Sigma_{2}=0$, the parameters must satisfy

$n_{20}=\frac{4w}{(w-1)(w+3)}n_{10}n_{30}.$ (43)

Note that the exact solution uses a different time parameterisation $\tau$,
with

$N=-\text{e}^{-b^{1}-b^{2}-b^{3}},$ (44)

and not $N=-\frac{1}{2\beta}$ as used in the numerical evolution. As mentioned
in the introduction, using the gauge condition (44) in numerical simulation
would come at a cost of fine-tuning. To avoid this, we shall keep using the
area time gauge in numerical simulation and evolve $\tau$ along the spike
worldline through (48) below to accommodate numerical matching.

##  4 Matching with exact solution

A new matching procedure is needed. We often rescale and shift the coordinates
variables to simplify the exact solution, but the numerical solution does not
necessarily appear in the simplified form. To accommodate numerical matching,
we relax the parameterisation of the coordinates for the exact solution, to
linear order:

$\tau\rightarrow k_{0}(\tau-\tau_{0}),\quad x\rightarrow k_{1}x,\quad
y\rightarrow k_{2}y,\quad z\rightarrow k_{3}z,$ (45)

where $k_{0}$, $k_{1}$, $k_{2}$, $k_{3}$ are positive constants, and
$\tau_{0}$ is a constant. This has the following effect on the metric
components:

$N\rightarrow k_{0}N,\quad b^{\alpha}\rightarrow b^{\alpha}-\ln
k_{\alpha},\quad n_{1}\rightarrow\frac{k_{2}}{k_{1}}n_{1},\quad
n_{2}\rightarrow\frac{k_{3}}{k_{1}}n_{2},\quad
n_{3}\rightarrow\frac{k_{3}}{k_{2}}n_{3}.$ (46)

But we shall maintain the time parameterisation (44), leading to the relation

$k_{0}=k_{1}k_{2}k_{3}.$ (47)

The spike solution is generated by the Geroch transformation, which leaves
arbitrary additive functions of $y$ and $z$ in $\beta_{2}$ and $\beta_{3}$,
and hence in $n_{1}$ and $n_{2}$ (see [art:Lim2015, Equations (27)–(29)]). But
because zooming reduces spatial dependence to essentially linear order, and
the fact that we will require the metric components to be even functions of
$Z$ in the section below, adding a constant in $n_{1}$ and $n_{2}$ will
suffice in most cases. To accommodate numerical matching, we must also evolve
the time variable $\tau$ of the exact solution along the spike worldline
$Z=0$, with evolution equation

$\frac{\text{d}\tau}{\text{d}T}=\frac{N_{\text{numerical}}}{N_{\text{exact}}}=\frac{\text{e}^{b^{1}+b^{2}+b^{3}}}{2\beta}=\frac{\text{e}^{b^{1}+b^{2}}}{2}E_{3}{}^{3}.$
(48)

We do not need to recover $\tau$ at places other than the spike worldline,
because the matching procedure below performs the matching _only along the
spike worldline_.

For completeness we shall evolve all the spatial metric components, which will
allow us to match every metric component. Their evolution equations in zooming
coordinates are111We take this opportunity to correct the errors in
[art:HeinzleUgglaLim2012, Equations (C1a)–(C1b)]). See
[art:HeinzleUgglaRohr2009, Equation (A.11b)] for the correct equations.

$\displaystyle\partial_{T}b^{1}$
$\displaystyle=-Z\partial_{Z}b^{1}+\tfrac{1}{2}(1+\sqrt{3}\Sigma_{-})$ (49)
$\displaystyle\partial_{T}b^{2}$
$\displaystyle=-Z\partial_{Z}b^{2}+\tfrac{1}{2}(1-\sqrt{3}\Sigma_{-})$ (50)
$\displaystyle\partial_{T}b^{3}$
$\displaystyle=-Z\partial_{Z}b^{3}+\tfrac{1}{2}(1-3\Sigma_{+})$ (51)
$\displaystyle\partial_{T}n_{1}$
$\displaystyle=-Z\partial_{Z}n_{1}-\sqrt{3}\Sigma_{\times}\text{e}^{b^{1}-b^{2}}$
(52) $\displaystyle\partial_{T}n_{2}$
$\displaystyle=-Z\partial_{Z}n_{2}-\sqrt{3}\Sigma_{\times}\text{e}^{b^{1}-b^{2}}n_{3}$
(53) $\displaystyle\partial_{T}n_{3}$
$\displaystyle=-Z\partial_{Z}n_{3}-\sqrt{3}\Sigma_{1}\text{e}^{b^{2}-b^{3}}.$
(54)

$\beta$ and $\Sigma_{+}$ have their own evolution equations, but we do not
need them since we can compute $\beta$ and $\Sigma_{+}$ algebraically using
$\beta=\text{e}^{b^{3}}/E_{3}{}^{3}$ and (21).

From [art:Lim2015], we see that the exact spike solution is multiply-
represented, where the same state-space orbit yields multiple values of $w$.
Here we shall give the new matching procedure using a small positive $w$, with

$0<w<1.$ (55)

Similar matching procedures can be developed for the other cases. We shall
focus on the joint spike transition only (described by the first alternative
in [art:Lim2015]), ignoring the uninteresting second alternative in
[art:Lim2015], which describes what is essentially an OT $G_{2}$ spike
transition preceded and succeeded by $\Sigma_{1}$ frame transitions. The new
matching procedure to determine the values of parameters $w$, $n_{10}$,
$n_{20}$, $n_{30}$, $\tau_{0}$, $k_{0}$, $k_{1}$, $k_{2}$, $k_{3}$ is as
follows.

Firstly, we determine the value of the parameter $w$. We plot the combination

$s_{1}(T)=(\Sigma_{-}+\tfrac{2}{\sqrt{3}})^{2}+\Sigma_{\times}^{2}$ (56)

along $Z=0$. For a joint spike transition, this combination behaves like a
sigmoid curve, transitioning from the value $\tfrac{1}{3}w^{2}$ to the value
$\frac{(w+3)^{2}}{3(w-1)^{2}}$. We thus obtain an approximate value for $w$
from $s_{1}$.

Next, the combination

$s_{2}(\tau)=(-6b^{1}-2b^{2}-2b^{3})|_{Z=0}=-\frac{1}{2}(w^{2}+3)k_{0}(\tau-\tau_{0})+\ln(k_{0}^{2}k_{1}^{4})$
(57)

allows us to obtain $k_{0}$ through its slope against $\tau$.

Thirdly, we want to match $e^{2b^{1}}|_{Z=0}$ with

$\displaystyle\frac{1}{k_{1}^{2}}\left[\text{e}^{(w-1)k_{0}(\tau-\tau_{0})}+n_{20}^{2}\text{e}^{-\frac{1}{2}(w^{2}-1)k_{0}(\tau-\tau_{0})}+n_{30}^{2}\text{e}^{-(w+1)k_{0}(\tau-\tau_{0})}\right]$
(58)
$\displaystyle=C_{1}\text{e}^{(w-1)k_{0}\tau}+C_{2}\text{e}^{-\frac{1}{2}(w^{2}-1)k_{0}\tau}+C_{3}\text{e}^{-(w+1)k_{0}\tau},$
(59)

where the coefficients $C_{1}$ , $C_{2}$ and $C_{3}$ are given by

$C_{1}=\frac{1}{k_{1}^{2}}\text{e}^{(w-1)k_{0}(-\tau_{0})},\
C_{2}=\frac{1}{k_{1}^{2}}n_{20}^{2}\text{e}^{-\frac{1}{2}(w^{2}-1)k_{0}(-\tau_{0})},\
C_{3}=\frac{1}{k_{1}^{2}}n_{30}^{2}\text{e}^{-(w+1)k_{0}(-\tau_{0})}.$ (60)

That is, we want to minimise the relative difference

$M(C_{1},C_{2},C_{3})=\sum\frac{(e^{2b^{1}}|_{Z=0}-C_{1}\text{e}^{(w-1)k_{0}\tau}-C_{2}\text{e}^{-\frac{1}{2}(w^{2}-1)k_{0}\tau}-C_{3}\text{e}^{-(w+1)k_{0}\tau})^{2}}{(e^{2b^{1}}|_{Z=0})^{2}},$
(61)

where the sum is done over selected numerical data points. $M$ is a sum of
squares and is quadratic in $C_{1}$, $C_{2}$ and $C_{3}$ with a single
critical point, which is a local minimum point. Its global minimum point is
located at the critical point, following the method of least squares. To find
the critical point, we solve the system

$\frac{\partial M}{\partial C_{1}}=0,\quad\frac{\partial M}{\partial
C_{2}}=0,\quad\frac{\partial M}{\partial C_{3}}=0.$ (62)

This yields the (unique) values for $C_{1}$, $C_{2}$ and $C_{3}$. $\ln
C_{1}^{2}$ plus a particular value for $s_{2}(\tau)$ eliminates $k_{1}$, and
can be solved to give $\tau_{0}$:

$\tau_{0}=\frac{\ln C_{1}^{2}+(-6b^{1}-2b^{2}-2b^{3})|_{Z=0}-\ln
k_{0}^{2}+\frac{1}{2}(w^{2}+3)k_{0}\tau}{[\tfrac{1}{2}(w^{2}+3)-2(w-1)]k_{0}}.$
(63)

Then $k_{1}$, $|n_{20}|$ and $|n_{30}|$ are obtained from $C_{1}$, $C_{2}$ and
$C_{3}$ through (60). $|n_{10}|$ is obtained from (43). Metric components
$n_{1}$ and $n_{2}$ are adjusted by an additive constant to match their
numerical counterpart:222This is because we have simplified the spike solution
by making $F_{2}(y,z)$ and $F_{3}(y,z)$ in [art:Lim2015] as simple as
possible. The numerical solution again does not necessarily take this simple
form. For our purpose here, a zeroth order adjustment (an additive constant)
is sufficient.

$n_{1}\rightarrow n_{1}+n_{1c},\quad n_{2}\rightarrow n_{2}+n_{2c},$ (64)

where the value of $n_{1c}$ is approximated by the final value of $n_{1}$
along $Z=0$, and similarly for $n_{2c}$. Next, we determine the sign of
$n_{10}$, $n_{20}$ and $n_{30}$. From (39), $n_{10}$ and $n_{3}$ have the same
sign, so we use $n_{3}$ to determine the sign of $n_{10}$. From (38), $n_{20}$
and $n_{2}$ have the same sign, so we use the adjusted $n_{2}$ above to
determine the sign of $n_{20}$. $n_{30}$ is obtained again from (43). $k_{2}$
is obtained from the relation

$\displaystyle k_{2}=k_{1}^{-1}\text{e}^{-b^{1}-b^{2}}\Big{[}$
$\displaystyle\text{e}^{-2k_{0}(\tau-\tau_{0})}+n_{10}^{2}\text{e}^{-\frac{1}{2}(w-1)^{2}k_{0}(\tau-\tau_{0})}$
$\displaystyle+(n_{10}n_{30}-n_{20})^{2}\text{e}^{-\frac{1}{2}(w+1)^{2}k_{0}(\tau-\tau_{0})}\Big{]}^{-1/2}.$
(65)

$k_{3}$ is obtained from the relation (47). This completes the matching
procedure for the case $0<w<1$. Similar matching procedures can be developed
for other ranges of $w$. Note that the matching procedure performs matching
along the spike worldline $Z=0$ only, which is the reason why we need to
evolve $\tau$ along $Z=0$ only in (48).

##  5 Results

As in [art:Limetal2009], we shall impose symmetry on the metric components to
hold the spike worldline fixed at $Z=0$. To form a non-moving true spike at
$Z=0$, we require that the metric components be even functions of $Z$.

We present a numerical run, in which a generic initial condition similar to
that in [art:Limetal2009, Section IV.B] is used. The goal is to show that the
joint spike transition observed in numerical simulation is described by the
non-OT $G_{2}$ spike solution.

The format for the initial condition is

$\displaystyle\tau=0,\quad
z=\frac{1}{2}(E_{3}{}^{3})_{0}Z+z_{\text{zoom}},\quad
E_{3}{}^{3}=(E_{3}{}^{3})_{0},$ (66)
$\displaystyle\Sigma_{1}=(\Sigma_{1})_{0}\left[\frac{(3a_{3}a_{6}-3a_{1}a_{5}-\sqrt{3}a_{5})z^{2}}{2(E_{3}{}^{3})_{0}}+\frac{3(a_{4}a_{6}-a_{2}a_{5})z^{4}}{4(E_{3}{}^{3})_{0}}\right]$
(67) $\displaystyle\Sigma_{-}=a_{1}+a_{2}z^{2},\quad N_{\times}=a_{5}z,$ (68)
$\displaystyle\Sigma_{\times}=a_{3}+a_{4}z^{2},\quad N_{-}=a_{6}z,$ (69)
$\displaystyle
b^{1}=(b^{1})_{0}-\frac{\sqrt{3}a_{5}z^{2}}{2(E_{3}{}^{3})_{0}},\quad
b^{2}=(b^{2})_{0}+\frac{\sqrt{3}a_{5}z^{2}}{2(E_{3}{}^{3})_{0}},$ (70)
$\displaystyle
b^{3}=(b^{3})_{0}-\frac{3(a_{3}a_{6}-a_{1}a_{5})z^{2}}{2(E_{3}{}^{3})_{0}}-\frac{3(a_{4}a_{6}-a_{2}a_{5})z^{4}}{4(E_{3}{}^{3})_{0}},$
(71) $\displaystyle
n_{1}=(n_{1})_{0}+\frac{a_{6}}{a_{5}}\text{exp}\left((b^{2})_{0}-(b^{1})_{0}+\frac{\sqrt{3}a_{5}z^{2}}{(E_{3}{}^{3})_{0}}\right),$
(72) $\displaystyle n_{2}=0,\quad n_{3}=0.$ (73)

We use the following values:

$\displaystyle a_{1}=1.3,\quad a_{2}=0.002,\quad a_{3}=0.3,$ (74)
$\displaystyle a_{4}=-0.001,\quad a_{5}=0.004,\quad a_{6}=-0.005,$ (75)
$\displaystyle(E_{3}{}^{3})_{0}=1,\quad(\Sigma_{1})_{0}=10^{-10},\quad
z_{\text{zoom}}=0,$ (76)
$\displaystyle(b^{1})_{0}=0,\quad(b^{2})_{0}=0,\quad(b^{3})_{0}=0,\quad(n_{1})_{0}=0,$
(77)

with $2001$ grid points over the $Z$ interval $[0,10]$ (exploit symmetry and
only simulate the right half of the spatial domain), and $T$ interval
$[0,40]$. In particular, $a_{1}$ is chosen such that we end up with a $w$
value that is not close to $1$, and $(\Sigma_{1})_{0}$ is chosen to be small
enough that we have a distinctive joint spike transition. It takes only a few
minutes to run on a single computer.

Figure 1: Plots of several variables relevant to the numerical run and
matching procedure.

A joint spike transition is observed after the generic data undergoes a
$\Sigma_{\times}$ frame transition. We use the data from $T\approx 9.8182$ to
$T=40$ for matching. Figure 1 plots six variables involved with the numerical
run and matching procedure, namely $\tau$ along the spike worldline $Z=0$ as
evolved by Equation (48), the lower bound for the right-hand excision boundary
$Z_{\text{right bound}}$ given by Equation (30) (which shows that the right
excision boundary $Z_{r}=10$ is greater than $Z_{\text{right bound}}$ for the
duration of simulation), the combinations $s_{1}(T)$ and $s_{2}(\tau)$ in
Equations (56)–(57), $\tau_{0}$ and $k_{2}$ expressions in Equations (63) and
(65).

Figure 2: The metric components for the numerical solution, the matching exact
spike solution, and their difference. Figure 3: The $\beta$-normalised
variables for the numerical solution, the matching exact spike solution, and
their difference. Figure 4: The Hubble-normalised Weyl scalars for the
numerical solution, the matching exact spike solution, and their difference.

The matching procedure yields the following values (rounded to 4 decimal
points):

$\displaystyle w=0.3134,\quad k_{0}=2.000,\quad\tau_{0}=-5.5372,$ (78)
$\displaystyle n_{10}=-2.7539\times 10^{-19},\quad n_{20}=-1.1273\times
10^{-14},\quad n_{30}=-7.4287\times 10^{4},$ (79) $\displaystyle
k_{1}=5.8354,\quad k_{2}=1.1053\times 10^{4},\quad k_{3}=3.1007\times 10^{-5}$
(80) $\displaystyle n_{1c}=-10.0260,\quad n_{2c}=4.3635\times 10^{-10}.$ (81)

Figures 2, 3, 4 plot the metric components, $\beta$-normalised variables and
Hubble-normalised Weyl scalars (see [art:Limetal2009, Appendix C] for their
formulas) for the numerical solution, the matching exact spike solution, and
their difference. They show that the relative difference is at the order of
$10^{-9}$ (except for $\ln\Sigma_{1}$, with difference growing to order
$10^{-7}$ at late times), which is quite good. This is a strong evidence that
the joint spike transition observed numerically, even when starting with a
generic initial condition, is well-matched by the exact spike solution. This
extends and completes the results of [art:Limetal2009]. While in
[art:Limetal2009] only the Weyl scalars are matched, here we also match the
metric components and the $\beta$-normalised variables.

Figure 5: State space orbits projected onto the Hubble-normalised
$(\Sigma_{+}^{H},\Sigma_{-}^{H})$ plane, showing a distinctive orbit for each
of the four groups of worldlines. The representative worldlines used are $\ln
z=-50,-41,-30,-20$. A red circle marks the start of the orbit, a red star
marks the end.

Recently in [art:LimMoughal2022], it was found that the non-OT $G_{2}$ spike
solution has four groups of worldlines with distinctive state space orbits. To
see this in the numerical solution, four additional numerical runs with $\ln
z_{\text{zoom}}=-50,-41,-30,-20$ and $T$ interval $[0,45]$ are made to produce
Figure 5, which shows the distinctive orbits for each of the four groups,
projected onto the Hubble-normalised $(\Sigma_{+}^{H},\Sigma_{-}^{H})$ plane.
Compare with Figure 6 of [art:LimMoughal2022].

##  6 Conclusion

When used in conjuction with the zooming technique for numerical simulations,
the area time gauge avoids the fine-tuning problem encountered in the fluid-
comoving or volume gauge used in [art:ColeyLim2016]. To recover the time
variable of the exact non-OT $G_{2}$ spike solution along the spike worldline,
it is evolved using Equation (48). A new matching procedure is developed for
the case $0<w<1$. Similar procedures can be developed for the other cases. We
have used a generic initial condition, where the metric components are even
functions of $z$ to hold the spike worldline fixed at $z=0$. Numerical
evolution from this initial condition shows a joint spike transition, which is
well-matched by the exact spike solution. This gives a strong evidence that
the non-OT $G_{2}$ spike solution indeed describes the joint spike transition.
This extends and completes the results of [art:Limetal2009], strengthening the
evidence that spike transitions are part of the generalised BKL dynamics. The
results of this paper can in principle be extended to stiff fluid $G_{2}$
models. The numerical code runs efficiently in a typical case ($w$ not too
close to 1). An alternative to zooming technique is the adaptive mesh
refinement technique, which requires a high performance computing cluster to
run, but which can deal with moving spikes and spikes beyond $G_{2}$ models,
as demonstrated in [art:GarfinklePretorius2020].