# Kink-Antikink Interaction Forces and Bound States in a $\phi^{4}$ Model with
Quadratic and Quartic dispersion

G. A. Tsolias  Department of Mathematics and Statistics, University of
Massachusetts,Amherst, MA 01003-4515, USA  Robert J. Decker  Mathematics
Department, University of Hartford, 200 Bloomfield Ave., West Hartford, CT
06117, USA  A. Demirkaya  Mathematics Department, University of Hartford, 200
Bloomfield Ave., West Hartford, CT 06117, USA  T.J. Alexander  Institute of
Photonics and Optical Science (IPOS), School of Physics, The University of
Sydney, NSW 2006, Australia  P. G. Kevrekidis  Department of Mathematics and
Statistics, University of Massachusetts,Amherst, MA 01003-4515, USA

###### Abstract

We consider the interaction of solitary waves in a model involving the well-
known $\phi^{4}$ Klein-Gordon theory, but now bearing both Laplacian and
biharmonic terms with different prefactors. As a result of the competition of
the respective linear operators, we obtain three distinct cases as we vary the
model parameters. In the first the biharmonic effect dominates, yielding an
oscillatory inter-wave interaction; in the third the harmonic effect prevails
yielding exponential interactions, while we find an intriguing linearly
modulated exponential effect in the critical second case, separating the above
two regimes. For each case, we calculate the force between the kink and
antikink when initially separated with sufficient distance. Being able to
write the acceleration as a function of the separation distance, and its
corresponding ordinary differential equation, we test the corresponding
predictions, finding very good agreement, where appropriate, with the
corresponding partial differential equation results. Where the two findings
differ, we explain the source of disparities. Finally, we offer a first
glimpse of the interplay of harmonic and biharmonic effects on the results of
kink-antikink collisions and the corresponding single- and multi-bounce
windows.

##  I Introduction

The study of nonlinear Klein-Gordon models is a topic that has a rich history.
Many of the early developments on the subject have focused on the
mathematically appealing theory of the inverse scattering transform and
integrable systems ablowitz ; ablowitz2 ; drazin , such as the famous sine-
Gordon equation gibbon ; oursg . However, more recently, the intriguing
features stemming from non-integrable dynamics have been at the center of
numerous studies centered around, e.g., the $\phi^{4}$ model ourp4 . The
latter has often been considered to be a prototypical system for phase
transitions, ferroelectrics, and high-energy physics among other themes gibbon
; ourp4 . Moreover, it has been a central point of both analytical and
numerical explorations, involving kink interactions, collective coordinates,
resonant dynamics (including with impurities) starting from the 1970’s and
extending over nearly 5 decades kud ; aub ; get ; AKL ; Sugiyama ; Campbell ;
belova ; Ann ; goodman ; goodman2 ; weigel ; weigel2  and even reaching to
this day weig1 ; clisthenis ; bazeia ; see also the recent recap of lizunova .

On the other hand, more recently, a diverse set of variants of the so-called
nonlinear beam (or biharmonic) wave equation have been considered also; a
collection of relevant examples includes, e.g., levandosky ; champneys ; CM ;
karageorgis . The corresponding models also span a diverse array of contexts,
including, e.g., suspension bridges and the propagation of traveling waves
therein. Another setting that has emerged very recently and has substantially
promoted the relevance of biharmonic models has been the form of generalized
nonlinear Schrödinger (NLS) type models in nonlinear optics in the context of
the so-called “pure-quartic solitons” pqs . Not only has this type of
dispersion engineering been realized in the lab, but it has also been used in
the context of the so-called pure-quartic soliton laser pqs3 . Mathematical
studies of the existence and stability of solitary waves of such equations are
also ongoing beam_demirkaya ; beam1 ; atanas .

It has, however, been recognized that the typical scenario involves (e.g., in
the NLS experimental settings discussed above) both regular quadratic and
quartic dispersion and that experimental settings have the ability to tune the
interplay between the two pqs2 . This results in the form of a generalized NLS
equation incorporating both quadratic and quartic dispersion and thus
presenting the potential for engineering a situation involving competition
between the two pqs2 . While the generalized NLS setting is the most canonical
one to consider in the realm of nonlinear optics, in the present work we opt
to consider the slightly simpler, yet highly informative, setting of a
corresponding Klein-Gordon model. The rationale behind the latter choice
involves the fact that the two models share the same existence properties, at
least in one spatial dimension, yet the nature of the real field-theory
renders the analytical calculations somewhat simpler, especially as regards
the stability and dynamical implications of the inter-wave interactions.

Our aim in the present work is to formulate the existence problem of a single
kink in a model incorporating quadratic and quartic dispersion in the presence
of a $\phi^{4}$ potential (this part will be entirely analogous to the
corresponding generalized NLS case). Subsequently, we intend to explore the
interaction of two such waves and identify their pairwise interaction force
and how it depends on the model parameters. Subsequently, conclusions of the
analytical theory will be tested against full numerical computations of the
interaction dynamics. Lastly, collisions between two coherent structures will
be simulated, and the possible scenarios thereof will be considered. Our aim
is to reveal the possibility that either the biharmonic effect may dominate
(yielding oscillatory tails and forces, equilibrium steady states of
alternating stability etc.) or the harmonic effect will prevail (featuring
exponential interactions and forces). The critical case between the two and
its own intriguing behavior will be revealed as well. In our study of
collisions and, in particular, in the case kinks and antikinks interact and
eventually separate, we create velocity-out versus velocity-in curves. These
curves show windows of velocity-in values for which we see different numbers
of bounces before the coherent structures separate. We compare this behavior
to both the “pure $\phi^{4}$” case of Eq. 1 and the “pure biharmonic” special-
case limits of the present model interpolating between them.

##  II Model Setup & Kink-Antikink Tail Behavior

The standard $\phi^{4}$ Klein-Gordon theory yields the field equation

$u_{tt}=u_{xx}-V^{\prime}(u)$ (1)

where $V(u)=\frac{1}{2}(u^{2}-1)^{2}$. In beam1 ; manton-decker , a variant of
this equation was explored where the harmonic spatial derivative term was
replaced by a biharmonic term of the form:

$u_{tt}=-u_{xxxx}-V^{\prime}(u).$ (2)

Here, as indicated in the above section, motivated by the corresponding
generalized NLS of pqs2 , we explore a model incorporating the competition of
the features of the two models:

$u_{tt}=\alpha u_{xx}-\beta u_{xxxx}-V^{\prime}(u)$ (3)

where $\alpha$ and $\beta$ are assumed positive (to ensure the competition
referred to above) and the potential function $V(u)$ is taken as before. When
we pick $\alpha=1$ and $\beta=0$, we get Eq. (1) and when we pick $\alpha=0$
and $\beta=1$, we get Eq. (2). Notice that while one of the coefficients could
be scaled out via a rescaling of space, we maintain both coefficients, in
order to maintain the tractability of the special case limits of $(0,1)$ and
$(1,0)$, i.e., biharmonic and harmonic respectively.

A central consideration of the present work is to explore both the features of
a single solitary wave, but also to examine the interaction between two such
waves, a kink and an antikink. We will use a method developed by Manton (as in
Manton_nuclear ; manton-decker ) to find the force between a separated kink
and antikink as a function of the separation distance. To do this we must
first determine the tail behavior for a single kink or antikink. Once the
force is determined, we can use the corresponding acceleration to generate an
ordinary differential equation (ODE), whose behavior can then be compared to
the soliton trajectories of Eq. (3), i.e., the corresponding partial
differential equation (PDE). As long as the separation distance between kink
and antikink remains sufficiently large, the agreement between ODE and PDE
should be quite good. However, in cases where the kink and antikink approach
each other at distances comparable to their respective widths, then it is no
longer obvious that the ODE model should be an adequate description of the
full PDE dynamics and the exchanges of energy between the different modes
present in the latter ourp4 . We will explore both the former agreement (at
large distances) and the latter deviations (at short ones) in the numerical
results below. In order to determine the tail behavior of a single kink we
proceed as follows. Substituting $\phi(x)=u(t,x)$ into Eq. (3) we get the
steady-state equation

$\alpha\phi^{\prime\prime}-\beta\phi^{(iv)}-V^{\prime}(\phi)=0$ (4)

To examine the relevant asymptotics, we substitute $\phi=1-\varepsilon
e^{\lambda x}$ into the Eq. (4). Neglecting terms of $\varepsilon^{2}$ and
higher (i.e., linearizing), for the above mentioned $\phi^{4}$ potential, we
get

$-\alpha\lambda^{2}+\beta\lambda^{4}+4=0$ (5)

It is easy to show that the roots of this equation are real for $\alpha\geq
4\sqrt{\beta}$ and complex for $\alpha<4\sqrt{\beta}$. In particular, for
$\alpha<4\sqrt{\beta}$:

$\lambda_{1}=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}+i\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}},\;\;\;\lambda_{2}=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}-i\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}$
$\lambda_{3}=-\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}+i\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}},\;\;\;\lambda_{4}=-\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}-i\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}$

for $\alpha=4\sqrt{\beta}$ (critical case), the degenerate roots are:

$\lambda_{1,2}=\sqrt{\frac{\alpha}{2\beta}},\;\;\;\lambda_{3,4}=-\sqrt{\frac{\alpha}{2\beta}}$

and for $\alpha>4\sqrt{\beta}$:

$\lambda_{1}=\sqrt{\frac{\alpha-\sqrt{\alpha^{2}-16\beta}}{2\beta}},\;\;\;\lambda_{2}=\sqrt{\frac{\alpha+\sqrt{\alpha^{2}-16\beta}}{2\beta}}$
$\lambda_{3}=-\sqrt{\frac{\alpha-\sqrt{\alpha^{2}-16\beta}}{2\beta}},\;\;\;\lambda_{3}=-\sqrt{\frac{\alpha+\sqrt{\alpha^{2}-16\beta}}{2\beta}}$

For real $\lambda$, similarly to the pure $\phi^{4}$ case, our model for the
tail behavior is

$be^{-ax}+de^{-cx},$ (6)

for the critical case with the double roots the model is

$be^{-ax}(x-d)$ (7)

(accounting for the relevant generalized eigenvector) and for the complex
$\lambda$ case (similarly to the pure biharmonic one), the model is

$be^{-ax}\cos(c(x-d)).$ (8)

We also know that in the real case $a=\lambda_{1}$
$=\sqrt{\frac{\alpha-\sqrt{\alpha^{2}-16\beta}}{2\beta}}$ and $c=\lambda_{2}$
$=\sqrt{\frac{\alpha+\sqrt{\alpha^{2}-16\beta}}{2\beta}}$ and in the complex
case $a={\rm Re}(\lambda)$
$=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}$ and $c={\rm
Im}(\lambda)$ $=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}$, while
in the critical case $a=\sqrt{\frac{\alpha}{2\beta}}$. We use curve fitting to
get the other parameters $b$ and $d$ . The results are in Table 1 for a
sequence of prototypical case examples that we have considered.

$\alpha$ | $\beta$ | $\lambda$ | Tail behavior
---|---|---|---
0 | 1 | 1-1i | $0.9700e^{-x}\cos(x-0.4083)$
1 | 1 | 1.1180-0.8660i | $1.205e^{-1.118x}\cos(0.8660(x-0.9909))$
2 | 1 | 1.2247-0.7071i | $1.793e^{-1.225x}\cos(0.7071(x-1.824))$
3 | 1 | 1.3229-0.5000i | $3.614e^{-1.323x}\cos(0.5(x-3.299))$
3.5 | 1 | 1.3693-0.3536i | $6.662e^{-1.369x}\cos(0.3536(x-4.942))$
4 | 1 | 1.4142 | $3.363e^{-1.414x}(x-0.9786)$
4.5 | 1 | 1.1042, 1.8113 | $4.451e^{-1.104x}-12.92e^{-1.811x}$
5 | 1 | 1, 2 | $3.354e^{-x}-25.64e^{-2x}$
6 | 1 | 0.8740, 2.2882 | $2.679e^{-0.8740x}-157.4e^{-2.288x}$
1 | 0 | 2 | $2e^{-2x}$

Table 1: Single Kink Tail Behavior for different model parameter
$(\alpha,\beta)$ in columns 1 and 2. Column 3 yields the corresponding
(spatial) eigenvalues and column 4 the functional form providing the optimal
fit to the tail behavior. One can read off the values of $a$, $b$, $c$, $d$ in
column 4 by referring to Equations (6), (7), (8).

.

In Figure 1 we show a single kink, $\phi_{K}(x)$, in the top three panels and
the tail-behavior of each kink in the bottom three panels for the cases
$\alpha=1$, $\alpha=4$, $\alpha=5$ with $\beta=1$ (respectively, left to
right). For the tail behavior we graph the right tail of $1-\phi_{K}(x)$
multiplied by $\exp(kx)$ as well as a model fitted to the tail of
$1-\phi_{K}(x)$, also multiplied by $\exp(kx)$ (appropriate for $x$
sufficiently large). The value of $k$ is equal to the real part Re$(\lambda)$
in the complex case, $\lambda$ in the critical case, and the smaller (in
absolute value) $\lambda$ in the real case (corresponding to the slow decay).
We can observe an excellent agreement in the oscillatory case (especially
factoring in that we have multiplied the expression by an exponential, hence
any deviation in the exponent would lead to an exponential growth). Similarly,
also a remarkable fit can be discerned even in the critical case, revealing
the underlying linear dependence modulating the exponential decay of the
generalized eigenvector in this setting. The exponential case (originally
doubly exponential turned into a single exponential upon multiplication by
$\exp(kx)$) is found to be less accurate. In the latter case of two real
$\lambda$ an improved fit can be obtained to the model $be^{-ax}+de^{-cx}$ if
$c$ is left as a free parameter in the curve-fitting process (rather than
using the value specified above which results from Eq. (5)). It is an
interesting question for future work, whether a weakly nonlinear theory can
capture more accurately the correction to the leading exponential dependence;
however, for our present purposes, the current prediction capturing adequately
the leading order exponential tail behavior will suffice.

Figure 1: The top three panels show a single kink, $\phi_{K}(x)$, for each of
the cases $\alpha=1$, $\alpha=4$, $\alpha=5$ and $\beta=1$ (respectively, left
to right). The bottom three panels illustrate the behavior of the right-side
tails of these kinks by graphing $1-\phi_{K}(x)$ multiplied be $\exp(kx)$,
where $k$ (positive) is Re$(\lambda)$ in the complex case, $\lambda$ in the
critical case, and the smaller (in absolute value) $\lambda$ in the real case.
Superimposed are the fitted curves, also multiplied by $\exp(kx)$. In all
cases the red solid curve is $\exp(kx)(1-\phi_{K})$. The blue dash-dot curve
is the fitted equation multiplied by $\exp(kx)$; the specific equations for
each case are $y=1.205\cos(0.8660(x-0.9890))$, $y=3.274(x-0.8606)$ and
$y=3.308-14.9e^{-x}$ (right to left respectively)

.

##  III THE FORCE BETWEEN KINK AND ANTIKINK

In order to find the force or acceleration between kink and antikink, we use
the approach of Manton as in Manton_nuclear ; manton-decker ; see also manton-
decker  for details of the calculation in the case where $\alpha=0$ and
$\beta=1$. We now briefly review some of the details of this force
calculation. Consider the momentum $P=-\int_{x_{1}}^{x_{2}}u_{t}u_{x}dx$ on
the interval $[x_{1},x_{2}]$ ($P$ is conserved when the integral is over the
entire real line). Differentiating under the integral and using Eq. (3) we
find that the force is given by

$\frac{dP}{dt}=F=\left[-\frac{1}{2}u_{t}^{2}-\alpha\frac{1}{2}u_{x}^{2}+\beta
u_{x}u_{xxx}-\beta\frac{1}{2}u_{xx}^{2}+V(u)\right]_{x_{1}}^{x_{2}}.$ (9)

For a field configuration that is static or almost so, we can ignore the first
term in the right hand side bracket of Eq. (9). We consider a configuration

$u(t,x)=\phi(x)=\phi_{K}(x+X(t))+\phi_{AK}(x-X(t))-1,$ (10)

where $-X(t)$ is the position of the kink and $X(t)$ is the position of the
antikink, which represents a kink-antikink pair approaching each other (as $t$
gets larger). Then, set $\eta=1-\phi$, $\eta_{K}=1-\phi_{K}$ and
$\eta_{AK}=1-\phi_{AK}$ (where $\phi_{AK}$ is an antikink solution to Eq.
(4)). Evaluating Eq. (9) from $x_{1}=x$ to $x_{2}\rightarrow\infty$ results in

$F=\alpha(\eta_{K})_{x}(\eta_{AK})_{x}-\beta(\eta_{K})_{x}(\eta_{AK})_{xxx}-\beta(\eta_{K})_{xxx}(\eta_{AK})_{x}+\beta(\eta_{K})_{xx}(\eta_{AK})_{xx}-4(\eta_{K})(\eta_{AK}),$
(11)

where we assume that the kink-antikink separation $2X$ is large. We have also
assumed that $|x|<<X$ and have used the approximation
$V(\phi)=V(1-\eta)\approx 2\eta^{2}$. Note that only the cross terms are left
at this point.

For the complex case the model for the tail behavior is
$\eta_{K}=be^{-ax}\cos(c(x-d)$, appropriate for $x$ sufficiently large.

Carrying out the derivatives in Eq. (11) and using $a={\rm
Re}(\lambda)=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}$ and $c={\rm
Im}(\lambda)=\frac{1}{2}\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}$ we get the
following expression for the force:

$F=-2\sqrt{\frac{16\beta-\alpha^{2}}{\beta}}b^{2}e^{-X\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}}\cos\left(\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}\left(X-d\right)+\theta\right).$
(12)

Here, we have that $\theta\in[0,\frac{\pi}{2}]$ such that
$\tan\theta=\frac{\alpha}{\sqrt{16\beta-\alpha^{2}}}$.

For the critical case the model for the tail behavior is $be^{-ax}(x-d)$.
This, upon substituting $a=\sqrt{\frac{\alpha}{2\beta}}$ and
$\alpha=4\sqrt{\beta}$ results in the force formula:

$F=-8b^{2}\sqrt{2\alpha}e^{-2aX}\left(X-\frac{\sqrt{2\alpha}}{4}-d\right),$
(13)

once again featuring a functional form reminiscent of that of the kink tail.

For the real case the tail behavior is $be^{-ax}+de^{-cx}$ and the force
becomes

$F=-b^{2}e^{-2Xa}\left(a^{2}\alpha-3a^{4}\beta+4\right)-d^{2}e^{-2Xc}\left(c^{2}\alpha-3c^{4}\beta+4\right),\vskip
12.0pt plus 4.0pt minus 4.0pt$

where:

$a=\sqrt{\frac{\alpha-\sqrt{\alpha^{2}-16\beta}}{2\beta}}\;\;\;\textnormal{and}\;\;\;c=\sqrt{\frac{\alpha+\sqrt{\alpha^{2}-16\beta}}{2\beta}}.$
(14)

Using these values the force formula can be written as

$F=\frac{\alpha^{2}-16\beta-\alpha\sqrt{\alpha^{2}-16\beta}}{\beta}b^{2}e^{-2Xa}+\frac{\alpha^{2}-16\beta+\alpha\sqrt{\alpha^{2}-16\beta}}{\beta}d^{2}e^{-2Xc}.$
(15)

Notice that the coefficient of the slow term is always negative while the
coefficient of the fast term is always positive. Dividing the above formulae
for the force by the mass gives the results in the acceleration column of
Table 2. The values for $b$ and $d$ are determined by curve fitting the tail
of a single kink, and are shown in Table 1.

$\alpha$ | $\beta$ | Mass | Acceleration
---|---|---|---
0 | 1 | 1.1852 | $6.351e^{-2x}\cos(2x-0.8166)$
1 | 1 | 0.9540 | $11.79e^{-2.236x}\cos(1.732x-1.464)$
2 | 1 | 0.8052 | $27.66e^{-2.4494x}\cos(1.4142x-2.0559)$
3 | 1 | 0.7031 | $98.30e^{-2.6458x}\cos(x-2.451)$
3.5 | 1 | 0.6633 | $259.1e^{-2.7386x}\cos(0.7072x-2.429)$
4 | 1 | 0.6290 | $407.0e^{-2.828x}(x-1.686)$
4.5 | 1 | 0.5991 | $166.2e^{-2.208x}-3769e^{-3.623x}$
5 | 1 | 0.5728 | $117.8e^{-2x}-27545e^{-4x}$
6 | 1 | 0.5287 | $92.75e^{-1.748x}-2194577e^{-4.576x}$
1 | 0 | 4/3 | $24e^{-4x}$
Table 2: Mass and acceleration as a function of the half-separation distance
$x$ of the kink and antikink.

Next, we integrate the expressions for the force on a kink to get the
potential energy for each case, and also find all fixed points and their
stability type. For the potential function in the complex case we have

$U=-b^{2}\sqrt{\frac{16\beta-\alpha^{2}}{2\sqrt{\beta}}}e^{-X\sqrt{\frac{4\sqrt{\beta}+\alpha}{\beta}}}\cos\left(\sqrt{\frac{4\sqrt{\beta}-\alpha}{\beta}}\left(X-d\right)+\vartheta\right)$
(16)

for $\vartheta\in[0,\frac{\pi}{2}]$ such that
$\tan\vartheta=\frac{\sqrt{4\sqrt{\beta}+\alpha}}{\sqrt{4\sqrt{\beta}-\alpha}}$.
For saddle points we have

$X^{*}(k)=d+\left(\frac{\pi}{2}+2k\pi-\theta\right)\sqrt{\frac{\beta}{4\sqrt{\beta}-\alpha}},$

and for centers we get

$X^{*}(k)=d+\left(\frac{\pi}{2}+(2k+1)\pi-\theta\right)\sqrt{\frac{\beta}{4\sqrt{\beta}-\alpha}}.$

For the critical case (repeated $\lambda$) the potential function is

$U=-2b^{2}\alpha
e^{\frac{-4\sqrt{2}}{\sqrt{\alpha}}X}\left(X-d-\frac{\sqrt{2\alpha}}{8}\right)$
(17)

and the center is given by

$X^{*}=d+\frac{\sqrt{2\alpha}}{4}.$

For the real case we have the potential function

$U=-b^{2}a\sqrt{\alpha^{2}-16\beta}e^{-2aX}+d^{2}c\sqrt{\alpha^{2}-16\beta}e^{-2cX},\
$ (18)

and the center:

$\frac{1}{2(c-a)}\left(\log{\frac{c^{2}d^{2}}{a^{2}b^{2}}}\right),$

with $a$ and $c$ as given in Eq. (14).

##  IV Comparison of ODE and PDE models

Using the expressions for the force we can now write an ODE for the time
evolution of the kink and antikink position for any given ($\alpha$, $\beta$)
combination. In Table 2 we have divided the force by the numerically
calculated mass and used the curve-fitted values for $b$ and $d$ to get an
acceleration expression for specific values of $\alpha$ and $\beta$. The
corresponding ODE is then $\ddot{X}=-dU/dX$, where the acceleration of the
right hand side is provided in Table 2. This ODE for the position of the one
coherent structure (while the other one is symmetrically located) is amenable
to a phase portrait analysis, as shown in Figure 2 and a comparison with the
corresponding PDE results of Eq. (3) can be obtained both at that level and at
the spatio-temporal evolution one as shown in Fig. 3.

Figure 2: $\alpha=1$, $\beta=1$, Phase portrait of the ODE Equation (ODE) in
comparison with Eq. (3) (PDE). The blue solid curve corresponds to
$X_{AK}(0)=8$, $\dot{X}_{AK}(0)=-0.02$. The red dash-dotted curve:
$X_{AK}(0)=8$, $\dot{X}_{AK}(0)=-0.00555$. The light blue closed orbit:
$X_{AK}(0)=3.3$, $\dot{X}_{AK}(0)=0$. The green curve corresponds to
$X_{AK}(0)=8$, $\dot{X}_{AK}(0)=-0.00356$. The pink solid closed orbit:
$X_{AK}(0)=7.4$, $\dot{X}_{AK}(0)=-0.0002$.

In Figure 2 we show trajectories for the case $\alpha=1$ and $\beta=1$
(complex case) that illustrate behavior near the steady states of the PDE (the
fixed points of the ODE). For these cases, there is clearly excellent
agreement, between the ODE and PDE phase planes which validates our force
calculations. The calculated force laws work well as long as the separation
between kink and antikink is sufficiently large. In the case of the PDE, we
identify the motion of the coherent structure by using the intersection of the
kink or antikink with the horizontal axis as the position, and also find the
corresponding speed, and thus extract an effective phase portrait to be
compared with the ODE results. In Figure 3 we show two of the trajectories
from Figure 2 as contour plots of the PDE with the ODE trajectory superimposed
on top (in blue). The left one among them is a robust oscillation around a
stable fixed point in the form of a center (the light blue curve in Figure 3).
The other is a trajectory that is scattered from the innermost potential
energy barrier –due to the presence of the innermost saddle point,
corresponding to a maximum– of the effective energy landscape and is thus
reflected. In this case, we see that the kinks do not make it to a collision
but are rather reflected due to their interaction landscape before the
collision.

Figure 3: $\alpha=1$, $\beta=1$. Comparisons of the PDE contour plot of the
displacement field $u(x;t)$ and the ODE trajectory. (blue solid curve). Left:
$x_{0}$ = 3.3, $v_{in}=0$. This corresponds to the closed orbit (light blue)
in Fig. 2. Right: $x_{0}$ = 8, $v_{in}=-0.02$. This corresponds to blue solid
curve in Fig. 2.

For another perspective on the quantification of the agreement between PDE and
ODE results, see Figure 4. In the left panel of this figure, we show the
potential energy plot of the ODE, again for $\alpha=1$, $\beta=1$, in blue
(using Eq. 16). The data points represent the potential energies of the steady
states of the PDE calculated using the PDE energy $E=\int{\frac{1}{2}\alpha
u_{x}^{2}+\frac{1}{2}\beta u_{xx}^{2}+\frac{1}{2}(u^{2}-1)^{2}dx}$ of the
associated steady state configurations. The corresponding steady-states
themselves are shown in the same figure, right panel. Note that since the
calculated potential energy curve of the ODE approaches zero as the separation
distance increases, the potential energies of the steady states of the PDE
must also be normalized (i.e., calibrated) so that the limiting value is zero
(by subtracting the potential energy of a steady state with very large
separation). Again, clearly, the local maximum and minimum values of the ODE
energy landscape coincide with the potential energies of the corresponding
steady states of the PDE. The local minima correspond to stable steady states
of the PDE (centers of the ODE) and the local maxima correspond to unstable
steady states of the PDE (saddle points of the ODE). Importantly, aside from
the center-most potential energy maximum where the kink structures are so
close that we cannot identify them as independent entities (and thus we do not
expect the collective coordinate characterization to be as accurate), we
observe that the agreement is very good. We remind the reader that the
presence of this oscillatory energy landscape, its associated minima (centers)
and maxima (saddles), and the respective stationary PDE configurations are
distinctive features of the prevalence of the biharmonic term and are
genuinely absent in the harmonic case (and more generally for
$\alpha>4\sqrt{\beta}$, when the harmonic contribution is dominant).

Figure 4: The left panel shows the energy vs $x_{0}$ for $\alpha=1$,
$\beta=1$. Blue curve is twice the potential function of the ODE for the
complex case (given in Eq. (16)). The data points are the (normalized)
potential energies of the steady states of the PDE at
$x_{0}=1.825,3.56,5.38,7.19,9.01$ which are shown in the right panel. The need
to multiply the potential function of the ODE by two when comparing ODE and
PDE stems from the fact that the energy calculation using a steady state of
the PDE involves two solitons - kink and antikink. The right panel presents
the static, equilibrium solutions corresponding to $x_{0}\approx 1.825$
(orange dashed-dot curve), $x_{0}\approx 3.56$ (blue solid curve),
$x_{0}=5.38$ (red dashed curve), $x_{0}\approx 7.19$ (green dashed-dot curve)
and $x_{0}\approx 9.01$ (purple dotted curve). Note that the steady state for
the PDE occurs at $x_{0}\approx 1.825$ but the fixed point of the ODE is at
$x_{0}\approx 1.75$.

While in Figure 2 we showed example phase portraits that resulted in very
proximal correspondence between PDE and ODE (for the complex case), in Figure
5 we show phase trajectories for both the real and complex cases that
illustrate at what point the PDE and ODE solution curves may depart from each
other (recall that the dash-dotted green line represents the position of the
antikink as measured by its intersection with the horizontal axis). In these
cases, the kink and antikink get too close for the force law to remain valid.
One can see that in the real case of the left panel this occurs at about
$X_{AK}=3$, while in the complex case of the right panel at about $X_{AK}=2$.
Note also that the green curve (position of the antikink) indicates the
formation of a bound state which is losing energy (in a way somewhat akin to a
stable spiral but keeping in mind that in a bound state there are no longer an
identifiable kink and antikink). Here, the important differences of the PDE
dynamics from the conservative ODE of 1 degree-of-freedom (dof) become
evident. The latter being energy conserving can only reflect in some way,
while the former can transfer energy from the kink translational motion to
other degrees of freedom (internal ones or radiation ones ourp4 ), thus
leading to the effective translational energy dispersion and thus the apparent
trapping of the kinks into a so-called bion state. Successive “breathings” of
this bion state at the PDE level are mirrored in the progressively inward
green curves (carrying less and less energy). At the ODE level, we make two
more minor (in terms of the bigger picture of our story), yet technically
relevant observations. Given the absence of ODE-PDE correspondence in the
right panel we stop the ODE evolution once the kink-antikink pair directly
collides (i.e., at $X=0$). On the other hand, the left panel has another
intriguing but non-physical trait: the double exponential force (of opposite
signs between the two exponentials) results in a landscape with a local
minimum very close to $X=0$. While this would result in a turning point of
$X\neq 0$, we have found this feature to be an artifact of the theory and its
lack of accuracy in the immediate vicinity of $X=0$. Let us reiterate, also in
light of the above remarks, that the ODE models are based only on the behavior
of the tails of the kinks and antikinks. When a kink and antikink are involved
in an interaction, the model makes sense only when the structures are well-
separated. Therefore the ODE model should not be expected to reflect the
actual behavior of the system beyond the point where the waves are at a
distance comparable to or smaller than their width, at which time they
essentially forego their individual character.

Figure 5:  The left panel shows phase plots for the real $\lambda$ case of
$\alpha=5$ and $\beta=1$ using initial conditions $X_{AK}(0)=8$ and
$\dot{X}_{AK}(0)=-0.003593$. The right panel illustrates phase plots for the
complex $\lambda$ case of $\alpha=1$ and $\beta=1$ using initial conditions
$X_{AK}(0)=8$ and $\dot{X}_{AK}(0)=-0.35$. In both cases, the ODE is shown by
the blue solid curve and the PDE by the green dash-dotted curve. Insets show
at what points the ODE model diverges from the PDE model. The ODE trajectory
in the right panel is stopped at the point when $X_{AK}=0$ because it becomes
physically unrealistic beyond that point.

In Figure 6 we show contour plots for the same values of $\alpha$ and $\beta$
and the same initial conditions as in Figure 5, again with the ODE solution
curves superimposed. I.e., these panels represent the spatio-temporal contour
plot representation of the failure of the ODE theory to capture the PDE
dynamics, as explained in the above discussion. As in Figure 5 we can see that
the ODE tracks the PDE simulation until around the time that the collision
occurs. As a result of the latter, at the PDE level, a bound state emerges,
while in the case of the ODE, the conservative nature of the 1 dof system does
not allow any scenario other than reflection.

Figure 6: Contour plots of the PDE corresponding to the same parameter values
and initial conditions as in the corresponding panels of Figure 5. The ODE
trajectory is superimposed in blue.

##  V Velocity in versus velocity out curves and Soliton Collisions

We now investigate kink-antikink collisions in the context of escape velocity
($v_{out}$) and multi-bounce windows as a function of incoming velocity
($v_{in}$), in line with the extensive literature on the subject discussed in
the Introduction (for a relatively recent summary in the $\phi^{4}$ case, see,
e.g., ourp4 ). For the latter model, it has been shown that there exists a
critical $v_{in}$ value, which we label $v_{crit}$, such that for
$v_{in}>v_{crit}$, the kink and antikink interact once and then separate
forever. For $v_{in}<v_{crit}$ the kink and antikink can form a bound state,
or can interact (bounce) any number of times, depending on $v_{in}$, before
separating forever. Furthermore, it is well-established since the work of Ann
for the $\phi^{4}$ model, that the bounce windows corresponding to different
numbers of bounces are nested in a fractal pattern. For example, three bounce
windows occur at the edges of two-bounce windows, four bounce windows occur at
the edges of three-bounce windows, and so on. Also, the $v_{in}-v_{out}$ graph
for a given window has the appearance of an inverted parabola, with the
$v_{out}$ values going to zero at the edges of the window.

For a model with only a biharmonic term ($\alpha=0$ in this paper) it was
shown in beam1  that two critical $v_{in}$ values exist, with
$v_{1,crit}<v_{2,crit}$. $v_{2,crit}$ is similar to $v_{crit}$ for the
$\phi^{4}$ model in that for $v_{in}>v_{2,crit}$ the kink and antikink
interact once and then separate. For $v_{in}<v_{1,crit}$ the kink and antikink
repel elastically before interacting. For $v_{1,crit}<v_{in}<v_{2,crit}$ the
kink and antikink form a bound state. Near both critical values, we see
oscillations in the $v_{in}-v_{out}$ graph, where the frequency of the
oscillations rapidly increases as the critical values are approached.

In the model described in this paper, which is a mixture of the two cases just
described, we find that by fixing $\beta$ at $\beta=1$ and letting $\alpha$
increase from zero to six, we see a transition from one case to the other. As
$\alpha$ is increased, more and more bounce windows begin to populate the
region between $v_{1,crit}$ and $v_{2,crit}$. At first these new bounce
windows display oscillations in the $v_{in}-v_{out}$ graphs near the edges of
the windows, similar to what is seen in the bound-state region of the pure
biharmonic case. With increasing $\alpha$ the oscillations diminish and the
$v_{out}$ values at the edges of each window begin to approach zero as in the
pure $\phi^{4}$ case. We will showcase these features qualitatively in the
results that follow. Nevertheless, the delicate nature of the associated
computations renders especially difficult the identification of effective
“critical points” where the behavior changes from the one reminiscent of the
pure biharmonic problem to that reminiscent of the pure harmonic one. In the
case of the critical velocities, by rescaling we are able to relate the
solutions to Eq. (3) for general $\alpha$ and $\beta=1$ to other combinations
of $\alpha$ and $\beta$, as is now shown.

Let $u^{1,\beta}$ be a solution to

$u_{tt}=u_{xx}-\beta u_{xxxx}+2u-2u^{3}.$

and consider the coordinate transformation

$x\mapsto\xi=\frac{x}{a}.$

In the new coordinate system the solution can be rewritten as
$u^{1,\beta}(x,t)=\tilde{u}(\xi,t)$. Of course,
$u^{1,\beta}_{tt}=\tilde{u}_{tt}$ and
$u^{1,\beta}_{x}=\frac{1}{a}\tilde{u}_{\xi}$, therefore $\tilde{u}$ obeys the
equation

$\tilde{u}_{tt}=\frac{1}{a^{2}}\tilde{u}_{\xi\xi}-\frac{\beta}{a^{4}}\tilde{u}_{\xi\xi\xi\xi}+2\tilde{u}-2\tilde{u}^{3}.$

For $a^{4}=\beta$ we get

$\tilde{u}_{tt}=\frac{1}{\sqrt{\beta}}\tilde{u}_{\xi\xi}-\tilde{u}_{\xi\xi\xi\xi}+2\tilde{u}-2\tilde{u}^{3}$

so, $\tilde{u}(\xi,t)\equiv u^{\alpha,1}(\xi,t)$ or,

$u^{1,\beta}(x,t)=u^{\alpha,1}\left(\frac{x}{\beta^{1/4}}\right)$

for $\alpha=\frac{1}{\sqrt{\beta}}$.

Therefore we can obtain solutions to the model for the parameters $\alpha=1$
and $\beta$, using the solution for parameters $\alpha$ and $\beta=1$. Then,
using this coordinate transformation, we find that the critical velocity
$v_{1,crit}^{1,\beta}$ of the solution $u^{1,\beta}$ can be expressed in terms
of the corresponding critical velocity $v_{1,crit}^{\alpha,1}$ of the
$u^{\alpha,1}$ solution as

$v_{1,crit}^{1,\beta}=\frac{dx}{dt}=a\frac{d\xi}{dt}=\beta^{1/4}v_{1,crit}^{\alpha,1},$
(19)

where $\alpha=\frac{1}{\sqrt{\beta}}$.

Notice that when $\beta$ becomes large enough, we get
$v_{1,crit}^{\alpha,1}\approx v_{1,crit}^{0,1}$, so

$v_{1,crit}^{1,\beta}\sim\beta^{1/4}v_{1,crit}^{0,1}.$

Similarly

$v_{1,crit}^{\alpha,1}\sim\sqrt{\alpha}v_{1,crit}^{1,0}.$

Furthermore, the above equations hold when $v_{1,crit}^{\alpha,1}$ is replaced
by $v_{2,crit}^{\alpha,1}$.

In Figure 7 we show graphs of $v_{1,crit}$ versus $\alpha$ for $\beta=1$ (left
panel) and $v_{1,crit}$ versus $\beta$ for $\alpha=1$ (right panel). Figure 8
is similar, but for $v_{2,crit}$. The blue circles on all panels are obtained
by the numerical simulation of Eq. (3) where the left panels represent
$v_{1,crit}$ vs $\alpha$ when $\beta=1$ and the right panels represent
$v_{1,crit}$ vs $\beta$ when $\alpha=1$. In the panels of both figures, the
red curves are obtained from the transformation given by Eq. (19) (plotted
without markers for the transformed points and with connecting lines in order
to make the graph more readable). The red curves are included to demonstrate
the validity of Eq. (19) in comparison with direct PDE simulations.

Figure 7: The blue circles on both panels are obtained by the numerical
simulation of Eq. (3) where left panel represents $v_{1,crit}$ vs $\alpha$
when $\beta=1$ and the right panel represents $v_{1,crit}$ vs $\beta$ when
$\alpha=1$. The red solid curve on the left panel is obtained by applying the
formula $v_{1,crit}^{\alpha,1}=\sqrt{\alpha}v_{1,crit}^{1,\beta}$ where
$\beta=\frac{1}{\alpha^{2}}$ to the numerically obtained data (blue circles)
on the right. The red solid curve on the right panel is obtained by applying
the formula $v_{1,crit}^{1,\beta}=\beta^{1/4}v_{1,crit}^{\alpha,1}$ to the
numerically obtained data (blue circles) on the left.

Figure 8: The left panel shows $v_{2,crit}$ vs $\alpha$ when $\beta=1$ and the
right panel shows $v_{2,crit}$ vs $\beta$ when $\alpha=1$. The blue circles
and the red solid curves were obtained as described in Figure 7.

.

Figure 9: The top left panel shows $v_{\rm out}$ vs $v_{\rm in}$ when
$\alpha=1$ and $\beta=1$ with $v_{2,crit}\approx 0.5902$. The top right panel
is the zoom-in about the first two-bounce curve. The bottom left panel is the
zoom-in about the two three-bounce windows right before the critical velocity
$v_{2,crit}$. The bottom right panel is the zoom in about the leftmost three-
bounce window on the bottom left panel. In both top right and bottom panels,
the tails and their oscillatory behaviors are shown. One-bounce windows in the
figures are in solid black. Two-bounce windows are in blue and three bounce
windows are in green. The gray solid line on the top left panel is when the
kink-antikink repel each other elastically.

Having identified the critical point scaling relations, we now turn to a
direct examination of the collision features and associated multi-bounce
windows. In Figure 9 we show $v_{in}-v_{out}$ curves for $\alpha=1$ and
$\beta=1$. For fixed $\beta=1$ we know that the force law changes from the
complex $\lambda$ case to the real $\lambda$ case at $\alpha=4$, so we expect
the case $\alpha=1$ to be somewhat similar to the pure biharmonic case of
$\alpha=0$. In the upper left panel of Figure 9 we see that the elastic
collision region corresponds to $0<v_{in}<v_{1,crit}\approx 0.30805$ and the
one-bounce region corresponds to $v_{2,crit}\approx 0.5902<v_{in}<1$. Note
that we have chosen not to include $v_{in}$ values greater than one. The
region $v_{1,crit}<v_{in}<v_{2,crit}$, which corresponds to a bound state when
$\alpha=0$, is beginning to be populated by two and three bounce windows, a
byproduct of the inclusion of the quadratic dispersion. The top right panel
shows the first two-bounce window. The bottom left panel shows the next two-
bounce window, with three-bounce windows appearing just to the left. The
bottom right panel shows the first three-bounce window in more detail. All
windows display the characteristic oscillations at the edges. Notice the
important features of this case: on the one hand, the multi-bounce windows
(which did not appear in the pure biharmonic case) are now present. On the
other hand, they do not terminate as, e.g., in the case of the standard
$\phi^{4}$ model ourp4 ; lizunova , but rather have the oscillatory
terminations (with progressively shorter periodicity) encountered previously
in beam1  for the pure biharmonic case, Moreover, we have encountered a
feature also absent in the standard (pure) $\phi^{4}$ case, namely higher-
bounce windows appear only on one side (to the left) of the two bounce
windows, while it is well-known Ann  that they appear on both sides in the
pure harmonic $\phi^{4}$ problem.

Figure 10: $v_{\mathrm{out}}$ vs $v_{\mathrm{in}}$ when $\alpha=5$ and
$\beta=1$ with $v_{c}\approx 0.7295$. The one-bounce window is in solid black.
Two-bounce windows are in blue and three-bounce windows are in green.

In Figure 10 we show $v_{in}-v_{out}$ curves for $\alpha=5$ and $\beta=1$. For
this case, since $\alpha>4$ we have $\lambda$ real, and expect some similarity
with the case of the pure $\phi^{4}$ model ($\beta=0$). Indeed, the structure
is similar to the fractal pattern we see in the $\phi^{4}$ case, with three-
bounce windows at the edges of the two bounce windows. However, we were not
able to find three-bounce windows to the right of the two-bounce windows.
These should, presumably, emerge as $\alpha$ gets larger, or as $\beta$ gets
smaller. However, it is an open question requiring further systematic
investigation how the self-similar (on both sides) picture of the pure
$\phi^{4}$ model arises.

Figure 11: Transition from dominant quartic progressively closer to dominant
harmonic behavior, by changing $\alpha$ from $2$ (left) to $3$ (middle) and
finally the critical case of $\alpha=4$ (right panel).

We can begin to see how the system transitions from the $\alpha=0$, $\beta=1$
(pure quartic dispersion) case to the $\alpha=5$, $\beta=1$ (harmonic term
dominant) case through some additional (less detailed) $v_{in}$-$v_{out}$
graphs; see, e.g., Figure 11. The first two-bounce window for $\alpha=2$ (left
panel, main figure, and see also Figure 9 top left panel for $\alpha=1$)
demonstrates a curious behavior. It appears (out of nowhere) at about
$\alpha=0.711$, persisting to about $\alpha=2.8$ where it disappears. This is
why it arises in the left panel, but not the middle one. Similarly, notice how
the transition shrinks progressively the size of the gray line interval of “no
collision” for $0<v_{in}<v_{1,crit}$. It can be seen that this interval
eventually disappears for $\alpha=4$ in the right panel of the figure, again
showcasing how the transition between the two regimes (biharmonic vs.
harmonic) emerges.

Figure 12: Contour plots when $\alpha=2.05$ and $\beta=1$ with $v_{c}\approx
0.6222$ for $X(0)=10$. Left panel is when $v_{\mathrm{in}}=0.622$ and right
panel is when $v_{\mathrm{in}}=0.621$.

While in the discussion above, we have focused on features that transition the
phenomenology between the two limits, it is important to realize that the
wealth of the model considered here transcends that of solely the limit cases.
For instance, Figure 12 illustrates a phenomenon not seen in either the pure
biharmonic or the pure harmonic $\phi^{4}$ case, with $\alpha=2.05$ and
$\beta=1$. What we see here is an initial interaction between kink and
antikink, followed by separation of the solitons for a period of time, and
then another approach of the pair. At this point one or more elastic
collisions can occur, resulting in the appearance of multiple bounces. In the
first panel of Figure 12 we see a “pseudo” two-bounce result, and in the
second panel a pseudo three-bounce result. This can occur when the speed at
which the kink and antikink approach each other for the second (or third) time
is very small and therefore when we find ourselves in the small-speed
reflection window of the complex eigenvalue case. In short, this is an
unprecedented type of two-bounce since two-bounces cannot happen in the pure
biharmonic case (where there is only bion formation and single bounce events
beam1 ), but it can also not happen for pure harmonic $\phi^{4}$ where the
small speed reflection scenario is absent. This is yet another manifestation
of the rich phenomenology of the model combining harmonic and biharmonic
dispersion.

##  VI Conclusions & Future Challenges

In this work we have explored a model featuring the competition of a harmonic
and biharmonic linear operator in a quadratic-quartic $\phi^{4}$ model. We
have argued that this model is of intrinsic mathematical interest due to the
distinct implications of the different linear operators and also the unique
features created by their interplay that neither of the “pure” (quadratic or
quartic dispersion) models possesses. The harmonic part creates a saddle point
in the spatial dynamics and hence leads to exponentially decaying waveforms.
On the other hand, the biharmonic operator leads to complex eigenvalues and a
spiral point in the corresponding spatial dynamics. Here, we have seen the
interplay of these two possibilities creating an effective competition between
the two tendencies. We have observed that this competition leads to a critical
point (with an intriguing behavior in its own right, i.e., a linearly
modulated exponential) and on the two sides of this criticality either the
biharmonic oscillatory effect or the harmonic exponential decay effect prevail
respectively. This crucially affects the interactions between the kinks which
we have also explicitly identified and corroborated by means of detailed
comparison of both the single wave tails and of the two-coherent-structure
interactions. We have also elucidated the extent to which this collective
coordinate approach can be reliably used and illustrated its failure when the
kinks get too close to each other. Additionally, we have examined the
collisions, bounce and multi-bounce windows stemming from the kink
interactions and have shown how the critical velocities and corresponding
windows are modified as a function of the quadratic-quartic model parameters.

These findings will clearly have a significant bearing on the corresponding
quadratic-quartic NLS model which is a natural possibility for future work,
given the recent developments in pure quartic solitons pqs ; pqs3  and on
optical media featuring quadratic and quartic dispersion pqs2 . This natural
extension will have similar existence properties to the case considered
herein, however the stability of the dark solitons of the latter problem is an
interesting open question. Aside from this, both at the Klein-Gordon (real
field-theoretic) level and at the NLS one, it will be interesting to
generalize considerations to higher dimensional settings, where it is natural
to expect, e.g., vortical and more generally topologically charged
excitations. Corresponding studies are currently in progress and will be
reported in future publications.

## Acknowledgments

This material is based upon work supported by the US National Science
Foundation under Grants PHY-1602994 and DMS-1809074 (P.G.K.).

## References

  * (1) M.J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform, SIAM (Philadelphia, 1981).
  * (2) M.J. Ablowitz and P.A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering, Cambridge University Press (Cambridge, 1991).
  * (3) P.G. Drazin, R.S. Johnson, Solitons: An introduction, Cambridge University Press (Cambridge, 1989).
  * (4) R. K. Dodd, J. C. Eilbeck, J. D. Gibbon, H. C. Morris, Solitons and Nonlinear Wave Equations (Academic Press, London, 1982).
  * (5) J. Cuevas, P.G. Kevrekidis, F.L. Williams (Eds.), The sine-Gordon Model and its Applications: From Pendula and Josephson Junctions to Gravity and High Energy Physics, Springer-Verlag, (Heidelberg, 2014).
  * (6) P.G. Kevrekidis and J. Cuevas-Maraver (Eds.), A Dynamical Perspective on the $\phi^{4}$ model, Springer Nature (Heidelberg, 2019).
  * (7) A. E. Kudryavtsev, Solitonlike solutions for a Higgs scalar Field, Pis’ma Zh. Eksp. JETP Lett. 22, 82 (1975).
  * (8) S. Aubry, Unified approach to interpretation of displacive and order-disorder systems 2\. Displacive systems J. Chem. Phys. 64, 3392 (1976).
  * (9) B.S. Getmanov, Bound States of Soliton in $\phi_{2}^{4}$ Field-Theory Model, JETP Lett. 24, 291 (1976).
  * (10) M. J. Ablowitz, M. D. Kruskal, and J. F. Ladik, Solitary Wave Collision, SIAM J. Appl. Math. 36, 428 (1979).
  * (11) T. Sugiyama, Kink-antikink collisions in the two-dimensional $\phi^{4}$ model, Prog. Theor. Phys., 61, 1550 (1979).
  * (12) D.K. Campbell, J.S. Schonfeld, and C.A. Wingate, Resonance structure in kink-antikink interactions n $\phi^{4}$ theory, Physica D, 9, 1 (1983).
  * (13) T.I. Belova and A.E. Kudryavtsev, Solitons and their interactions in classical field theory, Phys. Usp., 40, 359 (1997).
  * (14) P. Anninos, S. Oliveira, and R.A. Matzner, Fractal structure in the scalar $\lambda(\phi^{2}-1)^{2}$ theory, Phys. Rev. D, 44, 1147 (1991).
  * (15) R.H. Goodman and R. Haberman, Kink-antikink collisions in the $\phi^{4}$ equation: The $n$-bounce resonance and the separatrix map, SIAM J. Appl. Dyn. Sys., 4, 1105 (2005).
  * (16) R.H. Goodman, Chaotic scattering in solitary wave interactions: A singular iterated-map description, Chaos, 18, 023113 (2008).
  * (17) H. Weigel, Kink–antikink scattering in $\phi^{4}$ and $\phi^{6}$ models, J. Phys. Conf. Ser., 482, 012045 (2014).
  * (18) I. Takyi and H. Weigel, Collective coordinates in one-dimensional soliton models revisited, Phys. Rev. D, 94, 085008 (2016).
  * (19) N. S. Manton, K. Oleś, T. Roma?czukiewicz, A. Wereszczyński, Kink Moduli Spaces – Collective Coordinates Reconsidered, arXiv:2008.01026.
  * (20) C.F.S. Pereira, G. Luchini, T. Tassis, and C.P. Constantinidis, Some novel considerations about the collective coordinates approximation for the scattering of $\phi^{4}$ kinks, arXiv:2004.00571.
  * (21) D. Bazeia, Adalto R. Gomes, Fabiano C. Simas, Semi-compactness and multiple oscillating pulses in kink scattering, arXiv:2011.11157.
  * (22) M. Lizunova, J. van Wezel, An introduction to kinks in $\phi^{4}$-theory, arxiv:2009.00355.
  * (23) S. Levandosky,  Stability and instability of fourth order solitary waves,  J. Dynam. Diff. Eqs., 10, 151 (1998).
  * (24) A.R. Champneys, P.J. McKenna, and P.A. Zegeling,  Solitary waves in nonlinear beam equations: stability, fission and fusion,  Nonlinear Dynamics, 21, 31 (2000).
  * (25) Y. Chen and P.J. McKenna, Traveling waves in a nonlinearly suspended beam: theoretical results and numerical observations, J. Diff. Eqs., 136, 325 (1997).
  * (26) P. Karageorgis and P.J. McKenna, The existence of ground states for fourth-order wave equations, Nonlinear Anal., 73, 367 (2010).
  * (27) A. Blanco-Redondo, C. Martijn de Sterke, J.E. Sipe, T.F. Krauss, B.J. Eggleton, and C. Husko, Pure-quartic solitons, Nature Comms., 7, 10427 (2016).
  * (28) A.F.J. Runge, D.D. Hudson, K.K.K. Tam, C.M. de Sterke, A. Blanco-Redondo, The pure quartic soliton laser Nat. Photonics 14, 492 (2020).
  * (29) A. Demirkaya and M. Stanislavova, Numerical results on existence and stability of standing and traveling waves for the fourth order beam equation, Discrete Contin. Dyn. Syst. B, 24, 197 (2019).
  * (30) R. Decker, A. Demirkaya, P.G. Kevrekidis, D. Iglesias, J. Severino, and Y. Shavit, Kink dynamics in a nonlinear beam model, arXiv:2001.06973.
  * (31) I. Posukhovskyi and A. Stefanov,  On the normalized ground states for the Kawahara equation and a fourth order NLS, Discr. Cont. Dyn. Sys. A 40. 4131-4162 (2020).
  * (32) K.K.K. Tam, T.J. Alexander, A. Blanco-Redondo, and C.M. de Sterke, Generalized dispersion Kerr solitons, Phys. Rev. A , 101, 043822 (2020).
  * (33) R. J. Decker, A. Demirkaya, N. S. Manton and P. G. Kevrekidis, Kink-antikink interaction forces and bound states in a biharmonic $\phi^{4}$ model, J.Phys. A, 53, 37 (2020).
  * (34) N.S. Manton, An effective Lagrangian for solitons, Nucl. Phys. B, 150, 397 (1979).
  * (35) I.C. Christov, R. Decker, A. Demirkaya, P.G. Kevrekidis, and V.A. Gani, Long range interactions of kinks, Phys. Rev. D, 99, 016010, (2019).


Closed universes, factorization, and ensemble averaging

Mykhaylo Usatyuk, Ying Zhao

Kavli Institute for Theoretical Physics, Santa Barbara, CA 93106, USA

<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study closed universes in holographic theories of quantum gravity. We argue
that within any fixed theory, factorization implies there is one unique closed
universe state. The wave function of any state that can be prepared by the
path integral is proportional to the Hartle-Hawking wave function. This unique
wave function depends on the properties of the underlying holographic theory
such as the energy spectrum. We show these properties explicitly in JT
gravity, which is known to be dual to an ensemble of random Hamiltonians. For
each member of the ensemble, the corresponding wave function is erratic as a
result of the spectrum being chaotic. After ensemble averaging, we obtain
smooth semi-classical wave functions as well as different closed universe
states.

###### Contents

  1. 1 Introduction
  2. 2 Factorization and the uniqueness of closed universe state
    1. 2.1 Preparing states of closed universes
    2. 2.2 Hartle-Hawking wave function in a single theory
    3. 2.3 Ensemble averaging
  3. 3 From matrix integrals to closed universes
    1. 3.1 Review of JT gravity, matrix integral, and a factorizing theory
    2. 3.2 Closed universes in JT gravity
    3. 3.3 The unique closed universe wave function in a fixed theory
    4. 3.4 Ensemble averaging: Going back to JT gravity with smooth wave functions
  4. 4 Discussion
  5. A Review of factorization in deformed JT gravity dual to a fixed Hamiltonian

##  1 Introduction

Gauge/gravity duality is a powerful framework to study quantum gravity [1, 2,
3]. A lot of progress has been made in understanding quantum gravity in
asymptotically AdS spacetimes through the AdS/CFT correspondence. However, the
universe we live in has a positive cosmological constant and may not have an
asymptotic boundary. The understanding of closed universes without boundaries
has been limited.

In [4] we studied simple models of closed universes in two dimensional gravity
and pointed out various puzzling features. In particular, semi-classically
there appear to be many quantum gravity states encoding rich physics, but once
non-perturbative effects are included there is only a single closed universe
state [5, 6, 7, 8, 9, 10]. In this work we continue to study these issues. The
approach we will take is to begin with a well-defined holographic theory dual
to a certain bulk theory, and we will study closed universes in this theory.
In section 2 we point out the following general features.

  * •

With a single holographic theory:

    * –

Factorization implies there is only one closed universe state. Any state that
can be prepared by the path integral has a wave function proportional to the
Hartle-Hawking wave function.

    * –

The data of the boundary theory, such as the spectrum, is intricately encoded
in the unique wave function. For chaotic holographic theories, we expect the
wave function will be erratic.

    * –

A more refined argument using factorization shows that the unique wave
function, in a certain sense, contains trivial information. All states we can
input into the wave function are equivalent.

  * •

With an ensemble average over theories:

    * –

Factorization is no longer true. Distinct closed universe states can be
prepared by specifying different boundary conditions.

    * –

The wave functions contain non-trivial information. There are a large number
of orthogonal states that can be used as inputs for the wave functions.

    * –

The wave functions are smooth functions of the inputs, and are well
approximated by a sum over semi-classical geometries.

We will be able to show some of the features listed above hold generally,
regardless of the specific model being studied. For other features, we can
only rigorously show them in simple models. Nevertheless, we expect these are
general lessons that hold more broadly for holographic theories.

In section 3 we consider the concrete model of 2d JT gravity. Pure JT gravity
is known to be dual to an ensemble average of Hamiltonians [11]. We will
restrict to a single Hamiltonian by considering a specific deformation of the
theory[12]. We will demonstrate all of the points listed above in this model,
and find an explicit expression for the unique wavefunction in terms of the
boundary data. One surprising feature we will see is that even though we have
argued the closed universe wave function within a fixed single theory is
trivial, it nevertheless shares some qualitative similarities with the wave
function of the averaged theory.

In section 4 we point out unanswered questions and future directions.

##  2 Factorization and the uniqueness of closed universe state

###  2.1 Preparing states of closed universes

In this section we will assume we are starting with a single holographic
boundary theory dual to a bulk gravitational theory, and that the bulk theory
is specified by a bulk gravitational path integral. We are interested in
understanding the bulk theory on spatial slices that do not have any
asymptotic boundaries. In such cases, the connection between the bulk and
boundary theory is not well understood.

We begin by constructing quantum states of closed universes on a compact
Cauchy slice $\Sigma$. We allow for $\Sigma$ to be the union of multiple
disconnected compact slices. States are specified by a wave function of the
induced metric $h$ on the slice, along with other potential fields. The most
obvious state we can construct is the Hartle-Hawking state
$\Psi_{\text{HH}}[\Sigma]$, where we perform a path integral over all
manifolds bounded by a slice $\Sigma$ with data $h$ on the slice [13]

$\Psi_{\text{HH}}[\Sigma]=~{}\raisebox{-17.07182pt}{\includegraphics[width=56.9055pt]{Figure_new/sec2HH.png}}~{}~{}+\ldots$
(2.1)

This is not the only choice. We can construct other states by specifying
additional boundaries with boundary conditions denoted by $\Sigma_{B}$. Then
$\Psi_{\Sigma_{B}}[\Sigma]$ is given by the gravitational path integral over
all manifolds bounded by $\Sigma_{B}$ and the spatial slice $\Sigma$ with
induced metric $h$. In this notation the Hartle-Hawking state is specified by
the empty set $\Sigma_{B}=\varnothing$. Naively, by specifying different
boundary conditions we obtain different wave functions $\Psi_{\Sigma_{B}}$,
and hence different closed universe states. Let’s examine this more closely.
Consider the case where $\Sigma_{B}$ is an asymptotic boundary circle with
length $\beta$, which computes $Z(\beta)=\operatorname{Tr}(e^{-\beta H})$ in
the corresponding boundary theory. From the bulk path integral perspective it
appears we have produced a new wave function since there exist bulk manifolds
that connect $\Sigma_{B}$ to $\Sigma$, pictorially

$\Psi_{\beta}[\Sigma]=\raisebox{-48.36958pt}{\includegraphics[width=170.71652pt]{Figure_new/sec2HHdisk.png}}\,.$
(2.2)

However, since the boundary theory is a well-defined quantum mechanical
theory, $Z(\beta)$ is simply a number. This implies that for arbitrary spatial
slices $\Sigma$ we must have

$\displaystyle\Psi_{\beta}[\Sigma]=\Psi_{\text{HH}}[\Sigma]\times Z(\beta)\,.$
(2.3)

Another way to state the above is that the bulk path integral that computes
the cosmological state with additional asymptotic boundary conditions must
factorize, which is a variant of the standard factorization problem in
holography[6] applied to the case of computing the closed universe wave
function. One may worry that factorization may not hold here because $\Sigma$
is a finite boundary. However, from the perspective of the boundary quantum
theory, the inclusion of asymptotic boundaries $\Sigma_{B}$ computes an
independently well defined quantity, in this case $Z(\beta)$. As a result, we
must have that asymptotic boundaries factorize from the finite boundary
$\Sigma$.

This argument applies to arbitrary asymptotic boundary conditions
$\Sigma_{B}$. The boundary theory on asymptotic boundaries $\Sigma_{B}$ simply
computes an independently well defined boundary quantity $Z[\Sigma_{B}]$. By
factorization within a single theory we have

$\displaystyle\Psi_{\Sigma_{B}}[\Sigma]=\Psi_{\text{HH}}[\Sigma]\times
Z[\Sigma_{B}]\,.$ (2.4)

We conclude that for any fixed holographic theory there exists one closed
universe state: the Hartle-Hawking state. All other states prepared through
the bulk path integral have wave functions proportional to the Hartle-Hawking
wave function. In [4] it was argued that in a fixed theory the states of
closed universes are numbers. From the above argument these numbers are the
prefactors in front of the Hartle-Hawking wave function, given by
$Z[\Sigma_{B}]$ if the state is prepared with asymptotic boundary conditions
$\Sigma_{B}$.

We again emphasize that $\Sigma$ can be a union of disconnected compact
manifolds. $\Psi_{\text{HH}}[\Sigma]$ gives the unique closed universe wave
function with support on any number of disconnected universes.

###  2.2 Hartle-Hawking wave function in a single theory

A natural question to ask is how to compute the unique Hartle-Hawking wave
function given a well defined holographic theory. There are two approaches one
can take.

From the bulk perspective, if we know the exact bulk theory we can directly
perform the gravitational path integral following the procedure given by
Hartle and Hawking [13]. Additionally, given the complete description of the
bulk path integral we should find that factorization (2.4) automatically
holds.

From the boundary perspective, we would need to understand how the closed
universe is encoded in the boundary data. As the spatial slice $\Sigma$ does
not have a boundary, one may worry that the boundary theory knows nothing
about the closed universe. We find simple theories where this is not the case.
We will show in the next section that boundary data, such as the energy
spectrum of the Hamiltonian, is intricately encoded in the exact Hartle-
Hawking wave function. Going in reverse, given boundary data we can
reconstruct the Hartle-Hawking wave function for closed universes.

A natural guess is that these properties extend more generally beyond the
models we consider. If this is the case we can infer some general properties
of the wave function in a single theory. Holographic theories with bulk duals
that resemble Einstein gravity are chaotic [14]. Assuming the wave function
depends on the details of the spectrum, we should expect to find that
$\Psi_{\text{HH}}[\Sigma]$ is erratic for a single theory. However, performing
an ensemble average can smooth it out and give us wave functions that look
smooth and semi-classical.

#### Trivial wave function

We have concluded that given a single holographic theory, factorization
implies there is a unique wave function of any states prepared by path
integral. Naively this wave function carries interesting physical information,
for example it can tell us the relative amplitude to find a closed universe
with a particular field configuration on a Cauchy slice. This turns out to be
incorrect.

We can take two states specified by different field configurations
$\left(\Sigma_{1},h_{1}\right)$, $\left(\Sigma_{2},h_{2}\right)$ and consider
their inner product. The inner product is defined by performing a bulk path
integral over all manifolds bounded by $\Sigma_{1}$ and $\Sigma_{2}$ with
their specified data. Since the bulk path integral factorizes for a single
theory, we immediately have

$\displaystyle\bra{\Sigma_{1}}\ket{\Sigma_{2}}=\Psi_{\text{HH}}[\Sigma_{1}\cup\Sigma_{2}]=\Psi_{\text{HH}}[\Sigma_{1}]\Psi_{\text{HH}}[\Sigma_{2}]\,.$
(2.5)

Here we should emphasize that we are assuming the stronger condition that
finite boundaries also factorize in a single theory, see section 4 for
discussions. To check whether states are linearly independent and correspond
to distinct states in the Hilbert space, we can construct their Gram matrix

$\displaystyle
G_{ij}\equiv\bra{\Sigma_{i}}\ket{\Sigma_{j}}=\Psi_{\text{HH}}[\Sigma_{i}]\Psi_{\text{HH}}[\Sigma_{j}]\,.$
(2.6)

The rank of the Gram matrix tells us the number of linearly independent
states. From equation (2.6) one immediately sees that the gram matrix has rank
one. This implies that slices with different metrics/field content correspond
to the same physical state, and in this sense the unique wave function in a
single theory appears to be trivial.

###  2.3 Ensemble averaging

Let’s briefly summarize how ensemble averaging over boundary theories can
change the preceding story. It has been speculated that ensemble averaging may
play an important role in cosmology, see [15, 16, 17, 18, 10]. Our argument
for a single closed universe state crucially relies on the property of
factorization. With an ensemble average of theories, factorization no longer
holds, and we have the following consequences:

  * •

After ensemble averaging the gram matrix defined in (2.6) no longer has rank
one, which implies different field configurations can correspond to physically
distinct closed universe states.

  * •

Distinct states with different wave functions can be prepared by specifying
different asymptotic boundary conditions.

  * •

The closed universe wave functions can be smooth and behave semi-classically.

##  3 From matrix integrals to closed universes

In this section we will work with JT gravity which is dual to a matrix
integral. We will explicitly illustrate the statements made in the previous
section.

###  3.1 Review of JT gravity, matrix integral, and a factorizing theory

The JT gravity action in Euclidean signature is given by [19, 20,
21]111Throughout this section we set $8\pi G_{N}=1$.

$I_{\text{JT}}[g,\Phi]=-S_{0}\chi-\frac{1}{2}\int\sqrt{g}\Phi(R+2)-\int_{\text{bdy}}\sqrt{h}\Phi_{b}(K-1)\,,$
(3.1)

in terms of the dilaton $\Phi$ and metric $g$. The term $S_{0}\chi$ suppresses
higher topologies by the Euler characteristic. We have included a boundary
term so that good asymptotic boundary conditions consist of fixing the dilaton
$\Phi|_{\partial M}$ and the metric along the boundary.

Pure JT gravity is dual to an integral over random matrices [11]. More
precisely, the duality is

$\int
d\mu(H)\operatorname{Tr}\left(e^{-\beta_{1}H}\right)\ldots\operatorname{Tr}\left(e^{-\beta_{n}H}\right)=\int_{\text{b.c.}}\mathcal{D}g\mathcal{D}\Phi
e^{-I_{\text{JT}}[g,\Phi]}\,.$ (3.2)

On the left-hand side we have an integral over Hermitian matrices $H$ with a
particular measure $\mu(H)$.222The double scaling limit must be taken where
the size of the matrices is sent to infinity while we zoom into the end of the
spectrum. The measure is suitably tuned so that amplitudes computed in the
matrix integral match order by order in the topological expansion to JT
gravity amplitudes. The matrix $H$ is to be interpreted as the Hamiltonian of
a boundary quantum mechanical theory. Operator insertions in the matrix
integral, such as $\operatorname{Tr}\left(e^{-\beta H}\right)$, map to
gravitational boundary conditions such as an asymptotic boundary of length
$\beta$ [11, 22].

#### A factorizing model

While pure JT gravity is dual to a matrix integral, we are interested in
studying a bulk theory dual to a single quantum mechanical theory with fixed
Hamiltonian. Following [12, 23], we can consider a deformed bulk theory with
the following action:

$I[g,\Phi]=I_{\text{JT}}[g,\Phi]-\underbrace{\int_{0}^{\infty}dbb\mathcal{O}(b)Z_{H_{0}}(b)}_{\text{gives
discrete
spectrum}}+\underbrace{\frac{1}{2}\int_{0}^{\infty}bdb\mathcal{O}(b)\mathcal{O}(b)}_{\begin{subarray}{c}\text{gives
bulk factorization}\end{subarray}}\,.$ (3.3)

The operator $\mathcal{O}(b)$ when inserted in the path integral creates a
geodesic boundary of length $b$ on which surfaces can end.333This operator has
an expression in terms of dilaton gravity variables given by
$\mathcal{O}(b)\Leftrightarrow
e^{-S_{0}}\int_{\mathcal{M}}d^{2}x\sqrt{g(x)}e^{-2\pi\Phi(x)}\cos\left(b\Phi(x)\right)$.
However, it is more convenient to express it simply as $\mathcal{O}(b)$ when
performing the series expansion of the deformation. The function
$Z_{H_{0}}(b)$ is given by

$Z_{H_{0}}(b)=\sum_{i=1}^{N}\frac{2}{b}\cos\left(b\sqrt{E_{i}}\right)-\int_{0}^{\infty}dE\rho_{0}(E)\frac{2}{b}\cos\left(b\sqrt{E}\right)\,,$
(3.4)

where the energies $E_{i}$ are the eigenvalues of a Hamiltonian $H_{0}$ and
$\rho_{0}(E)=\sinh\left(2\pi\sqrt{E}\right)$ is the disk density of states of
pure JT. The second deformation term is non-local in spacetime, and introduces
a condensate of correlated geodesic boundaries. The role of this deformation
is to factorize the bulk theory by canceling all connected wormhole
contributions order by order in the genus expansion. It was argued that this
deformed bulk theory is dual to a single fixed boundary Hamiltonian $H_{0}$
whose discrete spectrum is encoded in the function $Z_{H_{0}}(b)$.

Amplitudes in the theory are defined by first defining a set of gravitational
boundary conditions, expanding the deformations in a perturbative series, and
evaluating the amplitude at each order in the expansion in the pure JT theory.
The gravitational boundary conditions can now be connected to the geodesic
boundaries generated by $\mathcal{O}(b)$ at each order in the expansion. It’s
simplest to represent the deformation terms through pictorial rules for the
bulk path integral444The interaction terms put into the action are
renormalized interaction vertices where all self-interactions have already
been defined to be re-summed. Thus the interaction vertices that are brought
down in the series expansion cannot connect to each other to create degenerate
surfaces.

$\displaystyle\text{Discrete spectrum
interaction:}\qquad\qquad\int_{0}^{\infty}dbb\mathcal{O}(b)Z_{H_{0}}(b)\Leftrightarrow~{}~{}\raisebox{-21.33955pt}{\includegraphics[width=31.2982pt]{Figure_new/ZH0.png}}$
(3.5) $\displaystyle\text{Factorization
branes:}\qquad\qquad-\frac{1}{2}\int_{0}^{\infty}dbb\mathcal{O}(b)\mathcal{O}(b)\Leftrightarrow~{}~{}\raisebox{-21.33955pt}{\includegraphics[width=62.59596pt]{Figure_new/Zfact.png}}$
(3.6)

The factorizing branes are pictured linked since their geodesic lengths $b$
are the same.555For the rest of the paper, we will use orange and red circles
to represent the geodesic boundaries $\mathcal{O}(b)$ generated by the
respective interaction terms. For the convenience of the readers, we briefly
explain how factorization works in this model in appendix A. Readers who are
interested in more details can look at the original paper [12].

###  3.2 Closed universes in JT gravity

We first look at closed universes in pure JT gravity [24, 25, 26, 4]. For
related recent works, see [8, 9, 27, 16, 17, 18, 28, 29, 30, 31, 32, 33]. The
unique classical solution in Lorentzian signature with compact spatial slices
$S^{1}$ gives a big bang/big crunch cosmology

$ds^{2}=-dt^{2}+b^{2}\cos^{2}(t)d\sigma^{2},\qquad\Phi(t)=\phi_{c}\sin(t)\,.$
(3.7)

The spatial coordinate is periodic $\sigma\sim\sigma+1$ and time runs from
$t\in(-\frac{\pi}{2},\frac{\pi}{2})$. The universe exists for a finite amount
of time, and reaches a maximal spatial size $b$ at time $t=0$ with the slice
being a geodesic.

To discuss the quantum theory of these geometries we can compute the Hartle-
Hawking wave function [13]. We integrate over all Euclidean geometries ending
on a slice with extrinsic curvature $K=0$ and geodesic size $b$ (blue circle
in (3.8)) . In JT gravity this wave function has a simple
expression[11]666Integrating the dilaton over a complex contour restricts the
integral to $R=-2$ geometries, so we must sum over constant negative curvature
geometries.

$\Psi_{\text{HH}}(b)=~{}~{}~{}\raisebox{-42.67912pt}{\includegraphics[width=56.9055pt]{Figure_new/Vg1.png}}+\ldots~{}~{}=\sum_{g=1}^{\infty}e^{S_{0}\chi}V_{g}(b)\,.$
(3.8)

given by a sum over hyperbolic manifolds bounded by a geodesic boundary $b$,
with each genus $g$ geometry contributing a Weil-Petersson volume weighed by
the Euler-characteristic.777The knowledge of Weil-Petersson volumes is not
necessary for this work, but closed form expressions are known. The final
result gives an amplitude to find a universe of size $b$.

The Hartle-Hawking prescription is one choice of asymptotic boundary
condition, namely no boundaries, from which we obtain a state. We can obtain
other states by demanding the insertion of additional boundaries $\Sigma_{B}$
in the Euclidean path integral. This would prepare a different closed universe
state depending on $\Sigma_{B}$. For example, we can insert an asymptotic
boundary of length $\beta$, and integrate over all bulk geometries ending on
the geodesic slice $b$ and asymptotic boundary888Explicit expressions for all
these amplitudes can be readily found in [11].

$\displaystyle\Psi_{\Sigma_{B}}(b)=~{}~{}~{}\raisebox{-42.67912pt}{\includegraphics[width=128.0374pt]{Figure_new/Vg1disk.png}}~{}~{}=V_{1}(b)Z_{\text{disk}}(\beta)+Z_{\text{trumpet}}(\beta,b)+\ldots\,.$
(3.9)

The bulk manifold connecting the $\Sigma_{B}$ boundary and the $b$ boundary is
a version of the factorization problem. In this case it is not a paradox since
we are explicitly averaging over boundary theories. The lack of factorization
allows us to construct many closed universe states $\Psi_{\Sigma_{B}}(b)$ by
varying the boundary conditions $\Sigma_{B}$. We shall now explicitly see this
is no longer the case in a theory that displays bulk factorization.

###  3.3 The unique closed universe wave function in a fixed theory

Let’s us now consider the wave function of closed universes in the single
fixed theory. To obtain the wave function we insert a geodesic boundary of
size $b$ and perform the path integral in a perturbative series in genus and
in the interactions. The Hartle-Hawking wave function is given by

$\displaystyle\Psi_{\text{HH}}(b)=~{}~{}\raisebox{-34.14322pt}{\includegraphics[width=170.71652pt]{Figure_new/HHfact.png}}$
(3.10)

It is not difficult to see see that all higher genus geometries will be
cancelled by the factorizing branes. Furthermore, the only geometric
contribution that survives is the one where the $Z_{H_{0}}(b)\mathcal{O}(b)$
vertex connects to our geodesic boundary. Such a geometry is a degenerate
cylinder of zero size sandwiched between geodesic boundaries of lengths $b$
and $b^{\prime}$. The amplitude for this degenerate cylinder is given by
$\langle b|b^{\prime}\rangle=\frac{1}{b}\delta(b-b^{\prime})$. The Hartle-
Hawking wave function is thus given by

$\displaystyle\Psi_{\text{HH}}(b)$
$\displaystyle=\int_{0}^{\infty}db^{\prime}b^{\prime}Z_{H_{0}}(b^{\prime})\times\frac{1}{b}\delta(b-b^{\prime})=Z_{H_{0}}(b)\,.$
(3.11)

We can attempt to modify this closed universe wave function by inserting
additional asymptotic boundaries. Let’s see how the condensate of
factorization branes ensures this does not work. The simplest geometric
contribution is given by connecting the asymptotic boundary to the geodesic
boundary,

$\Psi_{\Sigma_{B}}(b)\supset~{}~{}\raisebox{-34.14322pt}{\includegraphics[width=85.35826pt]{Figure_new/hhcancellation.png}}=0\,.$
(3.12)

However, as in (3.12), the factorizing branes give an identical geometry with
a relative sign, canceling this contribution. All geometries connecting the
asymptotic boundary $\Sigma_{B}$ and the geodesic circle $b$ cancel for
similar reasons. The only geometries that end up surviving are the following:

$\Psi_{\Sigma_{B}}(b)=~{}\raisebox{-19.91684pt}{\includegraphics[width=170.71652pt]{Figure_new/hhfinal.png}}~{}=\Psi_{\text{HH}}(b)\times
Z(\beta)\,.$ (3.13)

This is precisely the answer we should expect from a factorizing theory as
discussed in section 2. Thus we see we cannot change the state of the closed
universe by changing the boundary conditions. This same conclusion continues
to hold if we insert any other type of boundary conditions into the theory,
since the factorizing branes would kill off any kind of wormhole
contributions.

So far we have performed a bulk computation and found that in the factorizing
theory all closed universe states have wave functions proportional to the
Hartle-Hawking wave function

$\Psi_{\Sigma_{B}}(b)=\Psi_{\text{HH}}(b)\times(\text{constant})\,.$ (3.14)

Combining equations (3.4) and (3.11), we have

$\Psi_{\text{HH}}(b)=\sum_{i=1}^{N}\frac{2}{b}\cos\left(bE_{i}^{1/2}\right)-\int_{0}^{\infty}dE\rho_{0}(E)\frac{2}{b}\cos\left(bE^{1/2}\right)\,.$
(3.15)

This equation (3.15) directly relates the unique wave function of the closed
universe to the properties of the dual boundary theory. We point out a few
aspects of this relation.

  * •

The Hartle-Hawking wave function is fixed by the data of the boundary dual. In
this case, the boundary data consist of the energy eigenvalues $E_{i}$ of the
boundary Hamiltonian.

  * •

The wave function is erratic for a single member of the ensemble. For a
particular draw of the random matrix $H_{0}$ we can plot the wave
function:999The numerics were obtained by keeping the square root edge of the
spectrum of a random matrix, and cutting off the integral on the disk
subtraction at the largest eigenvalue of the random matrix. The matrix
ensemble was chosen to be the gaussian unitary ensemble instead of the JT
ensemble for simplicity.

Figure 1: The Hartle-Hawking wave function $\Psi_{\text{HH}}(b)$ for a single
draw of a random matrix $H_{0}$. Since the spectrum is chaotic the resulting
wave function is highly oscillatory and sensitive to the eigenvalues. In this
plot we take the first 100 eigenvalues out of a $10^{4}$ sized Hermitian
matrix, with typical level spacing $\sim 10^{-3}$ near the edge.

The erratic behavior comes from the fact that the spectrum for a random
Hamiltonian $H_{0}$ is chaotic.

  * •

We can extract boundary data if given the exact Hartle-Hawking wave function.

A priori it is unclear what (if any) information closed universes carry about
the dual theory on the boundary. After all, closed universes themselves do not
have an asymptotic boundary. We see in this theory that the Hartle-Hawking
wave function encodes boundary data in a particularly straightforward way.
Given the above wave function (3.15), we can extract the eigenvalues of the
boundary Hamiltonian by performing integral transforms. For example, we have

$\frac{1}{\pi}\int_{0}^{\infty}db\hskip
0.85355ptb\Psi_{\text{HH}}(b)\cos(b\sqrt{E})=\sum_{i=1}^{N}\quantity[\delta(\sqrt{E}-\sqrt{E_{i}})+\delta(\sqrt{E}+\sqrt{E_{i}})]-\rho_{0}(E)\,.$
(3.16)

In section 2 we also argued that this unique wave function is trivial in the
sense that all inputs of the wave function are the same, i.e., different $b$’s
correspond to the same state. In the current model this is obvious, as from
(3.15), with a fixed Hamiltonian $H_{0}$ the overlap between two geodesic
states is $\langle
b|b^{\prime}\rangle=\Psi_{\text{HH}}(b)\Psi_{\text{HH}}(b^{\prime})$ and (2.6)
manifestly holds. Said another way, $b,b^{\prime}$ are not distinct states.

###  3.4 Ensemble averaging: Going back to JT gravity with smooth wave
functions

We now briefly comment on the boundary matrix integral interpretation of these
wave functions. The matrix operator that corresponds to inserting a geodesic
boundary in gravity is given by[22, 12]101010The subtraction of the second
term gets rid of the disk level density of states contribution when the
operator is integrated over all random matrices. This is needed because in
gravity a hyperbolic disk with geodesic boundary does not exist, and the first
geometry contributing to the amplitude exists at genus one.

$\mathcal{O}(b)\Longleftrightarrow\frac{2}{b}\operatorname{Tr}\cos\left(b\sqrt{H}\right)-\int_{0}^{\infty}dE\rho_{0}(E)\frac{2}{b}\cos\left(bE^{1/2}\right)\,.$
(3.17)

The expectation value of this operator in the boundary theory is thus the
Hartle-Hawking wave function in the bulk. The deformations added to the bulk
action in (3.3) map to a deformation of the matrix integral potential (3.2),
and localize the integral onto a single Hamiltonian $H_{0}$ with eigenvalues
$\left(E_{1},\ldots,E_{N}\right)$ set by the parameter $Z_{H_{0}}$ [12]. This
reproduces the Hartle-Hawking wave function as computed from the bulk earlier.

The unique wave function with a fixed Hamiltonian has some undesirable
features: it is erratic, does not appear to carry non-trivial information, and
it is unclear how to do quantum mechanics with one state. These issues can be
solved by ensemble averaging over Hamiltonians to get back to pure JT. The
expectation of $\mathcal{O}(b)$ for each Hamiltonian will be erratic, but if
we average over the Hamiltonians we recover a smooth answer

$\int
d\mu(H)\mathcal{O}(b)\Longleftrightarrow\Psi_{\text{HH}}(b)=\sum_{g=1}^{\infty}e^{S_{0}\chi}V_{g}(b)\,.$
(3.18)

Furthermore, after averaging we can now prepare different cosmological states,
since the bulk path integral no longer factorizes. For example, we have that
the overlap between different $b$ basis states is now
$\bra{b}\ket{b^{\prime}}=\frac{1}{b}\delta(b-b^{\prime})$, whereas before it
was the product of two random numbers given by (3.15).

One natural question is whether the wave function of a single theory (3.15)
knows anything about the averaged wave function (3.18). The answer is yes but
in a subtle way. The full answer in (3.18) will grow as a function of $b$. One
may wonder if we can see this growth from the wave function of a single
theory. At least qualitatively, the answer is yes:

(a) Typical draw of $H_{0}$

(b) Atypical draw of $H_{0}$

Figure 2: We plot the Hartle-Hawking wave function $\Psi_{\text{HH}}(b)$ for a
typical and an atypical draw of the ensemble of Hamiltonians. At large $b$ the
atypical draw better approximates the bulk gravity answer. See footnote 12.

Figure 2 contains plots of the wave function (3.15). In 2(2(a)) we plot the
typical draw of the ensemble, while in 2(2(b)) we plot an atypical
draw.111111Plot 1 and plot 2(2(a)) are from the same draw of $H_{0}$. We
zoomed into figure 1 to illustrate the erratic nature of the wave function.
Here, by “atypical draw”, we mean it has at least one negative eigenvalue,
which causes growth in the expectation value of $\mathcal{O}(b)$ (3.17). The
probability to get a Hamiltonian with a negative eigenvalue depends on the
details of the theory. It is atypical in the sense that its spectrum contains
an eigenvalue that is negative after shifting the edge of the spectrum to take
the double scaling limit. This gives rise to growth as a function of $b$ which
is the expected bulk gravity answer.121212From the boundary perspective, the
exact wavefunction (3.17) will grow with $b$ if we have a negative eigenvalue.
A typical eigenvalue of size $\sim-e^{-S_{0}}$ would give growth
$\sim\frac{1}{b}\cosh\left(\\#be^{-S_{0}/2}\right)$. From the bulk
perspective, the Weil-Petersson volumes in (3.18) also grow as functions of
$b$. After doing a full resummation of the genus expansion along with the
inclusion of other doubly non-perturbative effects, the theory should see that
there are negative eigenvalues in the spectrum. We thus see a qualitative
match between the numerical and bulk computations, but these answers do not
quantitatively match. This is because the exact answer will be sensitive to
the entire probability distribution of negative eigenvalues.  One surprising
lesson is that certain features of the averaged wave functions actually come
from behaviors of atypical members of the ensemble.

We would like to emphasize that even though the wave function in a single
theory is unique, and as we have argued trivial, the coarse grained features
of the averaged wave function are captured by the wave function in a single
theory. In this sense even though there is a unique state it seems to capture
features of a sum over geometries in the averaged theory.

##  4 Discussion

In this paper we have studied closed universes in holographic theories of
quantum gravity. We have shown that in a single theory, factorization of
boundary partition functions implies that there is unique closed universe
state given by the Hartle-Hawking wave function. We have argued that for
holographic theories that are chaotic this wave function will be highly
erratic.

We have explicitly shown these properties in a simple two dimensional model of
quantum gravity given by JT gravity. One surprising conclusion is that the
wavefunction of closed universes is encoded in a highly non-trivial way in the
boundary data of the dual theory. We end with some open questions and
comments.

##### Ensemble averaging.

We see that ensemble averaging plays a crucial role in obtaining smooth and
rich semi-classical physics. Is it necessary? If it is, it will be quite
different from the standard gauge-gravity duality for spacetimes with
asymptotic boundaries. While some low dimensional bulk theories are known to
be dual to an ensemble average of boundary theories, how should one do
ensemble averaging in higher dimensional quantum mechanical theories like
super Yang-Mills? Another way to ask this is, how do smooth and semi-classical
wave functions emerge for theories for which there is no obvious ensemble
averaging procedure?

##### Observables in closed universes.

It is an interesting open problem to find clear observables for closed
universes. Since the wave function in a single theory is unique, and all input
states are secretly equivalent, it is unclear in what sense
$\Psi_{\text{HH}}(b)$ contains information. Nevertheless, we showed that if
you naively examine $\Psi_{\text{HH}}(b)$ for a single theory, it seems to
qualitatively reproduce on a coarse grained level some features of the semi-
classical wave function in the averaged theory. It thus seems tempting to
speculate that one may be able to extract non-trivial physics from
$\Psi_{\text{HH}}(b)$ even within a single theory.

Ensemble averaging significantly improves the situation. We obtain a large
number of orthogonal states and we can define some set of observable with
these states. However, naively defined observables seem to have large
fluctuations due to effects of spacetime wormholes [4] when the variance of
the observable is calculated. It would be interesting if there was a
mechanisms to suppress these fluctuations.

##### Expanding closed universes.

In this work we have taken the approach of starting from a well defined
holographic theory, and studied closed universes within the theory. The
simplest examples of closed universes we find with this approach are AdS big
bang/big crunch cosmologies [6, 18]. It would be interesting if there were
examples of holographic theories where the closed universe had other behavior,
such as expanding eternally like de Sitter space [34].

##### Factorization for finite boundary conditions.

In this work it was important that the bulk path integral factorizes when
evaluated with multiple boundary conditions. This includes both asymptotic and
finite boundaries. In the case of asymptotic boundaries, this must be true
since we are evaluating boundary CFT partition functions which are well
defined quantities that factorize. However, for boundary conditions that are
not asymptotic, such as finite boundaries used for the evaluation of the
Hartle-Hawking wavefunction, the argument for factorization is less clear. It
would be desirable if the path integral for finite slices was also non-
perturbatively defined in terms of some independent quantum mechanical theory.
In the case of JT gravity we found this is precisely the case; the geodesic
boundary is independently defined by the data of the boundary quantum
mechanics (3.17), and so must factorize.

Does this property continue to hold in higher dimensional theories such as the
bulk dual to $\mathcal{N}=4$ super Yang-Mills? One appealing approach is that
of Cauchy Slice Holography [35], where the quantum gravity wave function on a
non-asymptotic slice can be defined in terms of a $T\overline{T}$ flowed
partition function of the boundary theory evaluated on the slice.131313This is
currently being investigated in the case of JT gravity, where the quantum
theory at finite cutoff is a deformation of the Schwarzian theory[36]. This
would provide an independent quantum mechanical definition of the wave
function on non-asymptotic slices, and would explain the uniqueness and
factorization properties of such boundary conditions.

## Acknowledgment

We are grateful to Douglas Stanford for explaining the bulk origin of negative
eigenvalues, and Don Marolf for helpful comments on the draft. M.U. was
supported in part by grant NSF PHY-2309135 to the Kavli Institute for
Theoretical Physics (KITP), and by a grant from the Simons Foundation (Grant
Number 994312, DG). Y.Z. was supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958 and by a grant from the Simons
Foundation (815727, LB).

##  Appendix A Review of factorization in deformed JT gravity dual to a fixed
Hamiltonian

Let us first briefly explain how the factorization interactions in (3.3) leads
to a cancellation of all wormholes in the gravity path integral, and results
in factorization between disconnected boundaries. Let’s first compute the
leading contribution to the connected component of the product of two thermal
partition functions

$Z_{\text{conn}}(\beta_{1},\beta_{2})=\raisebox{-21.33955pt}{\includegraphics[width=227.62204pt]{Figure_new/trumpetfact.png}}$
(A.1)

The double trumpet contribution is cancelled out by a single factorizing
interaction. We see that Order by order the factorizing interaction will give
contributions that subtract off the wormhole contributions.141414The basic
idea behind the interaction is that all hyperbolic surfaces can be decomposed
into geometries where the asymptotic boundary is separated from the rest of
the geometry by a geodesic $b$. The introduction of
$-\frac{1}{2}\int_{0}^{\infty}dbb\mathcal{O}(b)\mathcal{O}(b)$ introduces
additional linked geodesic boundaries that cancel out any corresponding
wormhole. In [12] it was shown that these cancellations work to all orders in
the genus expansion. The only geometries that are not cancelled are those
which do not have a geodesic circle somewhere on the geometry, with the only
such geometry being the disk.

We now discuss the effect of the interaction vertex that gives a discrete
spectrum. Let us calculate the disk amplitude

$Z(\beta)=\raisebox{-21.33955pt}{\includegraphics[width=142.26378pt]{Figure_new/Zdisk.png}}$
(A.2)

The only surviving geometries to all orders in the genus expansion are the
disk and the trumpet weighed by $Z_{H_{0}}(b)$.151515As a reminder, self-
interactions for the vertices have already been re-summed so we do not have
have a factorizing vertex connect to the $Z_{H_{0}}(b)$ vertex to cancel out
this trumpet. We can obtain a discrete spectrum corresponding to a boundary
Hamiltonian $H_{0}=\text{diag}\left(E_{1},\ldots,E_{N}\right)$ by choosing

$Z_{H_{0}}(b)=\sum_{i=1}^{N}\frac{2}{b}\cos\left(b\sqrt{E_{i}}\right)-\int_{0}^{\infty}dE\rho_{0}(E)\frac{2}{b}\cos\left(b\sqrt{E}\right)\,,$
(A.3)

where $\rho_{0}(E)=\sinh\left(2\pi\sqrt{E}\right)$ is the disk density of
states of JT gravity, and the subtracted term has the effect of getting rid of
the disk contribution in (A.2). With this choice, we find

$Z(\beta)=Z_{\text{disk}}(\beta)+\int_{0}^{\infty}dbbZ_{\text{trumpet}}(\beta,b)Z_{H_{0}}(b)=\operatorname{Tr}e^{-\beta
H_{0}}\,,$ (A.4)

where we have used that
$Z_{\text{trumpet}}(\beta,b)=\frac{1}{\sqrt{\beta}}e^{-\frac{b^{2}}{4\beta}}$
and $Z_{\text{disk}}(\beta)=\int dE\rho_{0}(E)e^{-\beta E}$. We thus see that
the effect of the $Z_{H_{0}}(b)$ interaction is to get the bulk path integral
to reproduce a theory with a discrete spectrum given by a Hamiltonian $H_{0}$.
Formally, we should send the number of eigenvalues $N\to\infty$.

## References

  * [1] J. M. Maldacena, “The Large N limit of superconformal field theories and supergravity,” Adv. Theor. Math. Phys. 2 (1998) 231–252, arXiv:hep-th/9711200.
  * [2] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, “Gauge theory correlators from noncritical string theory,” Phys. Lett. B 428 (1998) 105–114, arXiv:hep-th/9802109.
  * [3] E. Witten, “Anti-de Sitter space and holography,” Adv. Theor. Math. Phys. 2 (1998) 253–291, arXiv:hep-th/9802150.
  * [4] M. Usatyuk, Z.-Y. Wang, and Y. Zhao, “Closed universes in two dimensional gravity,” arXiv:2402.00098 [hep-th].
  * [5] S. R. Coleman, “Black Holes as Red Herrings: Topological Fluctuations and the Loss of Quantum Coherence,” Nucl. Phys. B 307 (1988) 867–882.
  * [6] J. M. Maldacena and L. Maoz, “Wormholes in AdS,” JHEP 02 (2004) 053, arXiv:hep-th/0401024.
  * [7] G. Penington, S. H. Shenker, D. Stanford, and Z. Yang, “Replica wormholes and the black hole interior,” arXiv:1911.11977 [hep-th].
  * [8] D. Marolf and H. Maxfield, “Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information,” arXiv:2002.08950 [hep-th].
  * [9] D. Marolf and H. Maxfield, “Observations of Hawking radiation: the Page curve and baby universes,” JHEP 04 (2021) 272, arXiv:2010.06602 [hep-th].
  * [10] L. Susskind, “De Sitter Holography: Fluctuations, Anomalous Symmetry, and Wormholes,” Universe 7 no. 12, (2021) 464, arXiv:2106.03964 [hep-th].
  * [11] P. Saad, S. H. Shenker, and D. Stanford, “JT gravity as a matrix integral,” arXiv:1903.11115 [hep-th].
  * [12] A. Blommaert, L. V. Iliesiu, and J. Kruthoff, “Gravity factorized,” arXiv:2111.07863 [hep-th].
  * [13] J. B. Hartle and S. W. Hawking, “Wave function of the universe,” Phys. Rev. D 28 (Dec, 1983) 2960–2975. https://link.aps.org/doi/10.1103/PhysRevD.28.2960.
  * [14] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, “Black Holes and Random Matrices,” JHEP 05 (2017) 118, arXiv:1611.04650 [hep-th]. [Erratum: JHEP 09, 002 (2018)].
  * [15] T. Banks, W. Fischler, and S. Paban, “Recurrent nightmares? Measurement theory in de Sitter space,” JHEP 12 (2002) 062, arXiv:hep-th/0210160.
  * [16] J. Maldacena, G. J. Turiaci, and Z. Yang, “Two dimensional Nearly de Sitter gravity,” arXiv:1904.01911 [hep-th].
  * [17] Y. Chen, V. Gorbenko, and J. Maldacena, “Bra-ket wormholes in gravitationally prepared states,” JHEP 02 (2021) 009, arXiv:2007.16091 [hep-th].
  * [18] M. Van Raamsdonk, “Comments on wormholes, ensembles, and cosmology,” JHEP 12 (2021) 156, arXiv:2008.02259 [hep-th].
  * [19] C. Teitelboim, “Gravitation and Hamiltonian Structure in Two Space-Time Dimensions,” Phys. Lett. 126B (1983) 41–45.
  * [20] R. Jackiw, “Lower Dimensional Gravity,” Nucl. Phys. B252 (1985) 343–356.
  * [21] J. Maldacena, D. Stanford, and Z. Yang, “Conformal symmetry and its breaking in two dimensional Nearly Anti-de-Sitter space,” PTEP 2016 no. 12, (2016) 12C104, arXiv:1606.01857 [hep-th].
  * [22] A. Goel, L. V. Iliesiu, J. Kruthoff, and Z. Yang, “Classifying boundary conditions in JT gravity: from energy-branes to $\alpha$-branes,” arXiv:2010.12592 [hep-th].
  * [23] A. Blommaert, L. V. Iliesiu, and J. Kruthoff, “Alpha states demystified — towards microscopic models of AdS2 holography,” JHEP 08 (2022) 071, arXiv:2203.07384 [hep-th].
  * [24] D. Louis-Martinez, J. Gegenberg, and G. Kunstatter, “Exact Dirac quantization of all 2-D dilaton gravity theories,” Phys. Lett. B 321 (1994) 193–198, arXiv:gr-qc/9309018.
  * [25] J. Navarro-Salas, M. Navarro, and V. Aldaya, “Covariant phase space quantization of the Jackiw-Teitelboim model of 2-D gravity,” Phys. Lett. B 292 (1992) 19–24.
  * [26] M. Henneaux, “QUANTUM GRAVITY IN TWO-DIMENSIONS: EXACT SOLUTION OF THE JACKIW MODEL,” Phys. Rev. Lett. 54 (1985) 959–962.
  * [27] L. V. Iliesiu, J. Kruthoff, G. J. Turiaci, and H. Verlinde, “JT gravity at finite cutoff,” arXiv:2004.07242 [hep-th].
  * [28] S. Cooper, M. Rozali, B. Swingle, M. Van Raamsdonk, C. Waddell, and D. Wakeham, “Black hole microstate cosmology,” JHEP 07 (2019) 065, arXiv:1810.10601 [hep-th].
  * [29] M. Van Raamsdonk, “Cosmology from confinement?,” JHEP 03 (2022) 039, arXiv:2102.05057 [hep-th].
  * [30] S. Antonini, P. Simidzija, B. Swingle, and M. Van Raamsdonk, “Cosmology from the vacuum,” arXiv:2203.11220 [hep-th].
  * [31] S. Antonini, P. Simidzija, B. Swingle, and M. Van Raamsdonk, “Accelerating Cosmology from a Holographic Wormhole,” Phys. Rev. Lett. 130 no. 22, (2023) 221601, arXiv:2206.14821 [hep-th].
  * [32] A. Sahu, P. Simidzija, and M. Van Raamsdonk, “Bubbles of cosmology in AdS/CFT,” JHEP 11 (2023) 010, arXiv:2306.13143 [hep-th].
  * [33] S. Antonini, M. Sasieta, and B. Swingle, “Cosmology from random entanglement,” JHEP 11 (2023) 188, arXiv:2307.14416 [hep-th].
  * [34] S. Kachru, R. Kallosh, A. D. Linde, and S. P. Trivedi, “De Sitter vacua in string theory,” Phys. Rev. D 68 (2003) 046005, arXiv:hep-th/0301240.
  * [35] G. Araujo-Regado, R. Khan, and A. C. Wall, “Cauchy slice holography: a new AdS/CFT dictionary,” JHEP 03 (2023) 026, arXiv:2204.00591 [hep-th].
  * [36] R. Khan, M. Usatyuk, and W. Weng, “Work in progress,”.