Simulating the universe(s): from cosmic bubble collisions to cosmological observables with numerical relativity

The theory of eternal inflation in an inflaton potential with multiple vacua predicts that our universe is one of many bubble universes nucleating and growing inside an ever-expanding false vacuum. The collision of our bubble with another could provide an important observational signature to test this scenario. We develop and implement an algorithm for accurately computing the cosmological observables arising from bubble collisions directly from the Lagrangian of a single scalar field. We first simulate the collision spacetime by solving Einstein's equations, starting from nucleation and ending at reheating. Taking advantage of the collision's hyperbolic symmetry, simulations are performed with a 1+1-dimensional fully relativistic code that uses adaptive mesh refinement. We then calculate the comoving curvature perturbation in an open Friedmann-Robertson-Walker universe, which is used to determine the temperature anisotropies of the cosmic microwave background radiation. For a fiducial Lagrangian, the anisotropies are well described by a power law in the cosine of the angular distance from the center of the collision signature. For a given form of the Lagrangian, the resulting observational predictions are inherently statistical due to stochastic elements of the bubble nucleation process. Further uncertainties arise due to our imperfect knowledge about inflationary and pre-recombination physics. We characterize observational predictions by computing the probability distributions over four phenomenological parameters which capture these intrinsic and model uncertainties. This represents the first fully-relativistic set of predictions from an ensemble of scalar field models giving rise to eternal inflation, yielding significant differences from previous non-relativistic approximations. Thus, our results provide a basis for a rigorous confrontation of these theories with cosmological data.


Introduction
Do we live in a bubble? In a picture of eternal inflation driven by an inflaton field with multiple potential minima, our universe is predicted to lie inside one of many bubble "universes" nucleated out of an inflating "parent" vacuum (see Refs. [1,2] for a review of eternal inflation). An active area of inquiry is to search for observable signatures of eternal inflation, which might allow us to confirm this radical picture of the universe on the largest of scales. One unambiguous signature arises from the collision between bubble universes [3,4], which leaves distinctive cosmological signatures in the cosmic microwave background (CMB) radiation [5][6][7][8][9][10][11] and the distribution of large scale structure [12].
There has been enormous progress in understanding whether or not bubble collisions are a feasible observable consequence of the eternal inflation scenario. Previous work exploring the detailed outcome of bubble collisions [13][14][15][16][17][18][19][20][21][21][22][23][24][25][26][27][28] has firmly established that collisions can be a minor perturbation on standard cosmology, and also addressed the probability with which we can expect to see bubble collisions under various model assumptions [3,4,29,30]. A review of much of this previous work can be found in Refs. [18,31]. However, to date there has been no quantitative connection between a particular scalar field model giving rise to eternal inflation and the detailed signature of bubble collisions imprinted on the CMB. This paper aims to fill this gap, constructing for the first time a quantitative prediction for the signatures of bubble collisions in an ensemble of scalar field models. To this end, fully relativistic numerical simulations are required because the collision between bubbles is a highly non-linear process. We employ a powerful new set of numerical relativity codes, 1 significantly improving upon previous work [27] (see also Ref. [39]). Using a combination of a specialized gauge, adaptive mesh refinement, and a strategic choice of which regions to simulate, we are able to accurately simulate the collision itself, as well as the evolution of induced perturbations over the entire duration of slow-roll inflation inside each bubble. By tracing geodesics through the simulation, we directly extract the perturbed cosmological metric inside each bubble, allowing us to determine the observational signatures. We show that the signature in the CMB is well described by a set of four phenomenological parameters. The values of these parameters are only probabilistically determined, a consequence of the intrinsic variation in the relative position of the collision and the observer, the center of mass energy of the collision, and uncertainties in the underlying model. To provide a complete case study, we compute the probability distribution over these phenomenological parameters.
Observational searches for the signature of bubble collisions -and the interpretation of their findings -are only as good as the assumptions going into the theoretical predictions. Searches for the predicted signatures of bubble collisions have been performed [10,11,[40][41][42][43] on CMB data from the Wilkinson Microwave Anisotropy Probe (WMAP) [44]. Feeney et. al. [40] performed a Bayesian analysis, which constrained the number of observable collisions on the sky to be N < 4.0 at 95% confidence. Osbourne et. al. [42] took a different approach, constraining a combination of the amplitude a of the initial perturbation and the angular size θ crit of a bubble collision on the sky to lie in the range −4.66 × 10 −8 < a (sin θ crit ) 4/3 < 4.73 × 10 −8 [Mpc −1 ] at 95% confidence. Both analyses made assumptions about the form of the signature and the prior probability distribution over parameter values. The methods outlined in this paper will help guide future observational searches, as well as laying the foundations for translating those searches into constraints on fundamental physics.

Summary of background, methods, and results
Before launching into the technical details, let us pause to summarize the necessary background on bubble collisions in eternal inflation, outline the most important features of our methods, and highlight some of the most important results. The reader can use this section as a guide to the rest of the paper.
To accurately predict the observational signature of bubble collisions, one must in effect simulate the entire history of the early universe. The steps along the way are sketched in Fig. 1. One begins by specifying a model that drives eternal inflation. In this paper, we will consider a single scalar field with a potential containing a false vacuum and a single true vacuum. Given an initial condition where the field is at rest in the false vacuum over One begins by specifying a scalar-field potential with multiple vacua. The false vacuum is unstable to bubble formation. Each bubble has a wall determined by region 1 of the potential, and contains an open FRW universe, the evolution of which is determined by regions 2 and 3 of the potential. Surfaces of constant density are depicted as solid lines inside the bubble. Bubbles collide, and we can conveniently simulate them in a reference frame, the Collision Frame, where they nucleate at the same cosmological time. Because the field excursion caused by the collision is bounded, the dynamics of the collision are determined entirely by regions 1 and 2 of the potential. The subsequent cosmological evolution is determined by region 3. We reconstruct the perturbed FRW metric inside the observation bubble by evolving a set of geodesics through the simulation. A gauge transformation then allows us to extract the comoving curvature perturbation late in the inflationary epoch. Since this quantity is frozen in, and largely insensitive to the details of reheating, it is possible to calculate cosmological signatures such as the CMB temperature anisotropies. a region of size greater than H −1 F (the false vacuum Hubble radius), bubble nucleation will occur via the Coleman-de Luccia (CDL) instanton [45,46]. If the rate of bubble formation is somewhat less than one bubble nucleated per Hubble volume per Hubble time, then the phase transition will not complete, and inflation becomes eternal.
The rate of bubble formation and the profile of the bubble wall is entirely determined by the CDL instanton solution pertaining to the region of the potential labeled "1". The reader can find a summary of these instantons and the methods for constructing them in Sec. 3.1. The structure of the instanton also implies that the exterior of the null cone emanating from the bubble center contains values of the field between the instanton endpoints, labeled by the red dots in Fig. 1. The symmetry group of the instanton implies that each bubble contains an open Friedmann-Robertson-Walker (FRW) universe. In order to dilute the curvature and to produce phenomenologically-viable density fluctuations, we must invoke an epoch of slowroll inflation inside the bubble. These models of "open inflation" [47,48] have been widely discussed in the literature (see e.g., Ref. [49] for a review). We allow for an initial transient era described by evolution on the section of the potential labeled "2," before slow-roll occurs in region "3". This transient region is motivated by the hierarchy in scales necessary to support both slow-roll inflation and tunneling on the same potential 2 [50]. Inside a single bubble, evolution in regions "2" and "3" occurs along surfaces of constant density in the open FRW universe, as depicted in Fig. 1. It is important to note that throughout this paper we treat the field evolution in the post-nucleation region entirely classically, neglecting ambiguities about the transition between quantum mechanical and classical behavior after tunneling, 3 as well as fluctuations during slow-roll inflation and in the bubble wall.
Any given bubble will undergo collisions with an infinite number of other bubbles [4]. The subset of potentially observable collisions depends on the cosmology inside the bubble and the effect of collisions on it. A discussion can be found in Ref. [18]. For widely-separated collisions, we can consider each collision with the "observation bubble" (i.e., the one whose cosmology we are interested in) separately. The potential depicted in Fig. 1 allows for only one type of bubble, so all collisions are between identical bubbles. Extensions of this simplest model would allow for the collision between different bubbles, connecting the false vacuum to different basins of attraction. We restrict ourselves in this paper to this simplest class of models, leaving an assessment of the full range of possibilities to subsequent work.
The full collision spacetime for two bubbles possesses an SO(2,1) symmetry, allowing us to extract the spacetime (with metric Eq. 3.4) using a 1 + 1-dimensional simulation. This symmetry arises from the hyperbolic SO(3,1) symmetry of the individual bubbles (rotations and radial boosts): the intersection of two hyperboloids is a hyperboloid of one lower dimension. Rotations about the line separating the bubble centers is retained, as well as boosts transverse to this axis, yielding a residual SO(2,1) symmetry [14,15]. The equations of motion for the scalar field and metric (see Eq. 3.8) respect this symmetry, and unless additional effects are considered, there is no source for breaking it. One such source that violates our symmetry assumption is the effect of fluctuations on the bubble wall induced by tensor modes quantum-mechanically produced in the false vacuum de Sitter space (see e.g., Refs. [49,[52][53][54][55]). These fluctuations are generally small in amplitude, unless the background false vacuum has an energy density of order M Pl . More generally, even though instabilities on domain walls can arise, they become increasingly length-contracted, and hence unimportant, on expanding bubble walls as they reach relativistic velocities [56,57].
In a previous paper [27], three of the authors presented the first fully relativistic simulation of bubble collisions. In this paper, we use a significantly improved code which employs adaptive mesh refinement (AMR) and a new gauge which allows us to simulate the entire epoch of inflation inside each bubble. The results (see Fig. 7 and Fig. 8) are solutions for the metric functions and the scalar field, as functions of the coordinate x separating the bubble centers and a time variable N measuring the elapsed number of e-folds of expansion in the false vacuum. The code is convergent (see Fig. 9) and respects the constraints imposed by Einstein's equations (see Fig. 10).
Given the structure of the individual bubbles, regions 1 and 2 on the potential are most relevant for determining the immediate outcome of a collision. The immediate aftermath of the collision is quite simple when the bubble walls are relativistic: the field profiles making up the individual bubbles simply superpose [21,27]. Therefore, only a finite range in field space is relevant for predicting the immediate outcome of the collision. To evolve further, the full non-linearities of General Relativity become important. Because the bubbles collide in an ambient false vacuum de Sitter space, the center of mass energy for a collision is bounded, and the strength of the disturbance cannot grow arbitrarily large. Since the perturbation can be made small by tailoring the potential, we can ensure the the future of the collision is not a radical deviation from standard cosmological evolution. Region 3 of the potential will be important for determining the evolution of the background inflationary cosmology, and the evolution of the perturbation on that background sourced by the collision.
To simulate the collision and to extract the appropriate information requires multiple coordinate systems, summarized in Table 1. For the collision itself, we use coordinates that cover the interior of both bubbles, as well as the false vacuum they evolve in. To extract the signature of bubble collisions in the observation bubble, we must first go from the simulation coordinates to a set of perturbed FRW coordinates appropriate for describing the cosmology inside the observation bubble. We accomplish this by tracking a set of geodesics emanating from the nucleation center. Geodesics label coordinate positions {ξ, ρ, ϕ}, and evolve in proper time τ ; this yields a map 4 {N (ξ, τ ), x(ξ, τ )} between the simulation coordinates and the anisotropic slicing of the open FRW spacetime labeled by the geodesics. Geodesics are evolved numerically; an example of the map is shown in Fig. 12. We then find the metric in terms of coordinates {τ, ξ, ρ, ϕ} using the standard tensor transformation law, Eq. 4.31.

Name
Variables Metric False-vacuum Hyperbolic Coordinates Ψ, Υ, θ, φ Eq. 7.12 Used to label the nucleation sites of collision bubbles in the false vacuum. The simulations are performed in a reference frame where the colliding bubbles nucleate at the same time, which we refer to as the Collision Frame. To compute observational signatures, a convenient frame will be one in which the observer is at the center of the bubble. We refer to this as the Observation Frame. The two reference frames are connected by a diffeomorphism that can be enacted by a boost in an embedding-space picture (see Ref. [18]), which corresponds to a spatial translation in the cosmological spacetime contained inside the observation bubble, and to a time-translation in the spacetime outside the observation bubble.
Because the collision can only affect regions of the observation bubble in its causal future, the universe is split into regions affected and not affected by the collision. Near the causal boundary, the metric must, by continuity, be close to the empty open FRW metric. In this region, we can explicitly find the difference between the metric inside a bubble with a collision and a bubble without a collision, in order to map onto the standard metric fluctuation variables of cosmological perturbation theory. The resulting metric is explicitly in the synchronous gauge, and satisfies to a very good approximation the perturbative expansion of Einstein's equations (see Fig. 14). The synchronous gauge metric functions can be used to extract cosmological observables. However, a much cleaner variable to calculate is the comoving curvature perturbation R. For the adiabatic perturbations we consider, R remains constant on superhorizon scales and is insensitive to the details of reheating and subsequent evolution. Fig. 16 shows the comoving curvature perturbation near the end of inflation as a function of the Cartesian coordinates X, Y, Z in an open universe. As can be seen in this figure, R varies mainly along the direction connecting the bubble centers, with deviations from this near-planar symmetry only becoming apparent on scales close to the radius of curvature of the open FRW spacetime. Near the collision boundary, we find that R is well described by a power law in the coordinate ξ, as in Eq. 5.27. This allows us to develop a set of phenomenological parameters that describe the collision. The detailed profile depends on the initial separation ∆x; the free parameters in the template Eq. 5.27 vary with ∆x as shown in Fig. 19.
The comoving curvature perturbation can be converted into the observed temperature and polarization signatures in the CMB. In this paper, we work in the Sachs-Wolfe [58] approximation, where the temperature anisotropy is proportional to the comoving curvature on the surface of last-scattering, evaluated on the past light-cone of an observer, δT /T = R(x ls )/5. By symmetry, collisions project onto a disc on the CMB sky. From the detailed form of R, the temperature profile is a power law in the cosine of the angular distance from the center of the disc, Eq. 6.6. This fitting function agrees with the full projection of R over most scales, as shown in Fig. 20 and Fig. 21.
We have mentioned two sources of variation in the observed signature after the underlying inflationary potential is fixed. In the Collision Frame, as depicted in Fig. 2, variations in the observed signature arise from variations in the position of the observer relative to the causal boundary of the collision ξ obs and from variations in the initial separation of the bubbles ∆x. To these, we add two parameters describing variations in the underlying scalar-field potential itself. First, we can perform an overall rescaling of region 1 of the potential by a factor β; second, we can alter region 3 so as to effect a change in the observable portion of the surface of last scattering R ls . This yields four fundamental parameters {ξ obs , ∆x, β, R ls }, which we translate into four phenomenological parameters describing the characteristics that determine the observational signatures of this ensemble of scalar-field models: the observed angular scale θ c , number of collisions N coll , power-law index κ, and R ls . These parameters are analogous to, for example, the amplitude and spectral index predicted by an inflationary model -complicated scalar-field models map onto a few phenomenological parameters. A main result of this paper is a derivation of the probability distribution over these phenomenological parameters, Eq. 7.24. The marginalized probability distribution for a fiducial class of models is plotted in Fig. 23.
The probability distribution, along with the phenomenological parameters and fitting function describing the signal, are the complete prediction of this set of models. Similar distributions exist for other models -or classes of models -and are the central input for any observational search. In future work, following the procedure outlined above will allow Observer Position Kinematics Figure 2. Observer position and kinematics: the two intrinsic sources of variation in the collision signature in the Collision Frame. An observer's position today, as well as the intersection of their past light-cone with the surface of last scattering, are depicted using filled dots. Moving the observer along this constant-time hypersurface (left panel to right panel), more of the region affected by the collision is seen. The causal boundary (dashed light-cone emanating from the collision) and perturbed region projects onto a disc of variable size on the observer's CMB sky. Alternatively, fixing the observer's position in its bubble, and decreasing the separation of the colliding bubbles, more of the region affected by the collision is seen, and the strength of the underlying perturbation will change. Both of these effects project onto a variable size, amplitude, and perturbation profile in the affected region on the sky. us to identify how features of the underlying model interface with observational constraints.

CDL Instantons
To simulate the collision of two bubbles in an inflating background we first need the initial conditions: a field and metric configuration for each bubble that interpolates between the metastable false vacuum state and a region with lower vacuum energy separated from the false vacuum by a potential barrier. If one of the bubbles is to represent our observable universe, its interior must undergo a period of slow-roll inflation.
Tunneling between vacuum states was first studied in detail by Coleman and de Luccia [45,46]. Their basic conclusion was that quantum tunneling occurs via bubble nucleation -a bubble with a lower energy vacuum state fluctuates into existence in a sea of false vacuum -with an exponentially-suppressed tunneling rate. The rate, as well as the field configuration -and with gravity, the metric -of the bubble is found by solving Euclidean field equations to find an "instanton" solution that extremizes the Euclidean action.
Once nucleated, the bubble's internal pressure will accelerate its bubble wall outwards, quickly approaching the speed of light. The instanton with the smallest action, and thus the dominant contribution to the tunneling rate, comes from the bubble which is O(4)-symmetric in the Euclidean metric.
The general CDL instanton with O(4) symmetry has a Euclidean metric of the form where r is the bubble's (Euclidean) radius and dΩ 3 is the metric on a three-sphere. The function ρ(r) determines the curvature. For example, an instanton in pure de Sitter space has ρ(r) = R dS sin(r/R dS ), where R dS is the de Sitter radius. Such an instanton has a maximum radius of πR dS , corresponding to half the circumference of the de Sitter space. Any instanton that interpolates between regions in a positive-definite potential will likewise be compact with some maximum radius r max . The instanton equations take the form: 3) The field φ should be in different basins of attraction (near different minima in the potential) at r = 0 and r = r max in order to achieve a non-trivial solution to the equations of motion. Conventionally, the field is in the lower energy vacuum at r = 0. Non-singular solutions must also satisfy dφ/dr| r=0 = dφ/dr| r=rmax = 0 and ρ(r = 0) = ρ(r = r max ) = 0. We also specify dρ/dr| r=0 = 1, which essentially just sets ρ and r to have the same units. Note that an overall rescaling of V just rescales the radius r (and ρ); it does not otherwise affect the instanton's shape. We solve the instanton equations using an iterative overshoot/undershoot method. First, we make an initial guess for the value of the field φ at r = 0. The guess should have V (φ guess ) < V (φ metastable ). The equations of motion dictate that the field will roll up the potential barrier as r increases, gaining momentum. If the field has enough momentum (|dφ/dr| is large enough), it will continue down the other side of the barrier and, if the initial guess was correct, stop with dφ/dr = 0 exactly at r = r max near the metastable vacuum. If the field does not have enough momentum, it will stop rolling before r = r max and reverse course towards the top of the potential barrier. This is an undershoot, and the initial guess must be modified to have a lower vacuum energy and sit closer to the stable minimum. If the initial guess is very close to a minimum (or an inflection point) such that dV /dφ is very small, the field will not move much over a large radius and the friction term 3 ρ dρ dr dφ dr will become small or negative. Therefore, the field can always be made to have enough momentum to cross the barrier. If the momentum is too large, the field will either cross the metastable minimum or not stop (dφ/dr = 0) by r = r max . These are overshoots, and the initial guess must be modified to have a higher vacuum energy and sit closer to the potential barrier. One can converge upon the correct initial condition by making a sequence of initial guesses, integrating the equations of motion, and then modifying the guesses in response to overshoots and undershoots.
We implement the overshoot/undershoot method using the CosmoTransitions software package [59] with several modifications to account for gravity. Including the dynamics for ρ (whereas ρ(r) = r in flat space) is relatively straightforward; the criterion for an undershoot V (φ) Figure 3. The inflationary potential used throughout this study, rescaled such that V (φ F ) = 3 8π (its original scale only affects the conversion between dimensionless simulation coordinates and physical scales). In order to achieve ∼ 60 e-foldings of inflation, the potential requires a large slow-roll region (left), while the requirement that the instanton's size be much smaller than the de Sitter radius forces the potential barrier to be very narrow (magnified on the right).  remains unchanged, but we had to add an additional criterion for an overshoot. The field overshoots the desired end point if dφ/dr = 0 when ρ(r) = 0, even if the field has not passed the metastable vacuum. Note that if the bubble wall radius is of similar size to, or larger than, the de Sitter radius, then the only solution may be the Hawking-Moss instanton [60]. In this case, the field is constant and resides at the top of the potential barrier separating the two vacua. The modified CosmoTransitions code will converge to this solution when appropriate.
For the purposes of this paper, we consider just a single potential with both bubbles coming from the same type of instanton. The potential is identical to the 'L2' potential in Ref. [27], replotted here in Fig. 3. The metastable false vacuum is at φ F = 0; there is an inflection point at φ I = 0.0074 M Pl ; and the true vacuum (with zero cosmological constant) is at φ T = 2.75 M Pl . Fig. 4 shows the shape of the associated instanton. A broad survey of different scalar field potentials will appear in future work.

Simulation equations and initial conditions
Each bubble by itself is invariant under rotations and boosts about its origin. When there are two bubbles, the SO(3, 1) symmetry is broken down to SO(2, 1): rotation about, and boosts transverse to, the vector separating the two bubbles are still symmetries of the field configuration. Therefore, we can fully model the collision with a (1 + 1)-dimensional simulation. We choose coordinates with metric where H F is the false vacuum Hubble constant (and H −1 F is the false vacuum de Sitter radius), and α(N, x) and a(N, x) are unknown metric functions. In the pure de Sitter case, α = a = 1 and N measures the number of false vacuum e-foldings. Inside the bubbles, where α, a = 1, the variable N does not correspond to the number of slow-roll e-foldings in a straightforward way. Instead, the number of e-foldings can be found by taking N in concert with the metric functions α and a. Distances and times are measured in units of H −1 F ; this is a natural set of units to describe the collision, since all dimensionless quantities will be order one. For the potential we simulate in this paper, the inflationary energy scale is comparable to the false vacuum Hubble scale, and we expect a and α to be order one until the simulation reaches the end of slow-roll inflation inside the bubbles. However, this will not necessarily be true for other potentials which have a large hierarchy between the energy scale of the false vacuum and slow-roll inflationary epoch inside the bubbles. The coordinate x is measured in radians such that the point at (N = 0, x = 0) is antipodal to the point at (N = 0, x = π). The SO(2, 1) symmetry is manifest in the coordinates χ and ϕ: a boost with rapidity η transverse to the x-axis will transform a point at χ = 0 to χ = η, and a rotation about the x-axis simply shifts the coordinate ϕ.
These coordinates are much better suited to the simulation of collisions than those used in Ref. [27]. The relations between the coordinates and metric functions used by Ref. [27] and in this work are given by: N = sinh −1z , α new = α old cosh N , and a new = a old / cosh N . It can be seen that we avoid exponential changes in the lapse and shift, and can follow the simulation through the end of inflation (N ∼ 60).
The equation of motion for the field φ is given by φ = ∂ φ V , and the equations for the metric functions α and a come from the Einstein equations. These were worked out in detail in Ref. [27]. We restate the results here in terms of our new gauge choice. First, let us define where φ is measured in units of the Planck mass, φ = dφ/dx, and V (φ) has been rescaled such that V (φ F ) = 3 8π . Then the field and metric functions evolve as The solutions to these equations should satisfy the constraint The initial field configuration comes almost directly from the instanton calculation. In the limit where N → 0, the field-dependent terms in the equations of motion Eq. 3.8 are sub-dominant. At N = 0, the metric functions are α(N = 0, x) = a(N = 0, x) = 1. Comparing metric equations 3.1 and 3.4 near N = 0, we can see that for a single bubble, after appropriately rescaling φ and r, we simply have φ(x, N = 0) = φ(r = |x|). Since the instanton is symmetric in time, Π(N = 0) = 0. A second instanton can be added by translating the first instanton along the x-axis and then summing the fields. Technically, a single instanton modifies the field and the metric across the entire compact Euclidean space. Its effect upon a second instanton's field configuration is generally non-negligible. However, we shall assume that for cases of interest, the two instantons can be considered to have finite size, and as long as their separation significantly exceeds the sum of their radii, they can be considered independently.
We start the simulation at some small value of N , where the metric Eq. 3.4 becomes arbitrarily close to a slicing of Minkowski space with SO(2,1) symmetry. The spatial dependence of the metric functions will be negligible in this limit, where anisotropic curvature dominates the Einstein equations. For further discussion of this point, see Ref. [27]. Performing a Taylor expansion about N = 0, the lowest order surviving terms are Inserting these into the equations of motion Eq. 3.8, we obtain where a prime denotes differentiation with respect to x. We then use this as the input for the simulation.

Integration and adaptive mesh refinement
We integrate the above equations of motion using the method of lines (see e.g., Ref. [61] for a general discussion of the numerical methods incorporated below). In the simplest possible implementation, one would discretize the field in the x-direction on a uniform grid and find the x-derivatives using finite differences. One would then integrate the field in time using any standard numerical integrator, with the field's value at each discrete point being a separate integration variable. This was the technique used in Ref. [27], and with sufficient resolution it should work here as well. Unfortunately, as the bubbles expand, their walls quickly approach the speed of light and become severely length-contracted. The shape of the wall is extremely important for accurately reproducing the perturbed metric, as we will show later, so we need extremely high resolution in its surrounding region. This is not feasible with a uniform grid. So instead, we switch to a method involving adaptive mesh refinement (AMR). AMR describes computational algorithms which introduce resolution dynamically where and when such resolution is needed. A number of variations appear in the literature, and the one developed here can be considered a variant of that introduced by Berger and Oliger [62]. In common with that work, this scheme adapts both the spatial resolution ∆x and the time-step ∆N . However, it is important to note the following points, which are unique to this approach and appropriate for just a single spatial dimension: (i) spatial resolution is continuously adaptable instead of discretely, (ii) only a single set of points covers any particular region of the spatial domain, (iii) AMR boundaries are handled in the spirit of the 'tapered boundary' approach of Ref. [63], and (iv) regridding occurs only globally when all regions are time-aligned. Our algorithm proceeds as follows: 1. The first step, after setting the initial conditions for the simulation, is to create the initial grid. To do this, the algorithm calculates a desired grid spacing at each of the input points. The grid density should be high near the bubble walls where the field and metric functions change rapidly, and relatively low elsewhere. We set the grid density m ≡ dn grid /dx to be directly proportional to dφ/dx, such that the number of grid points along the bubble wall is n wall = wall m dx. We also enforce a minimum grid density m min . In this paper, we use n wall = 70 and m min = 500. It is helpful if the grid density does not vary too quickly: otherwise, sharp features in the simulation could migrate from high resolution regions to low resolution regions before the algorithm recomputes the grid. We set an additional minimum of for all points i and j. This guarantees that there will be at least n 2 points along the grid before the density drops by a factor of 2. We take n 2 = 40. With the grid density set, the code integrates m to create the new grid points.
2. Next, the values of the field and metric functions need to be computed along the new grid. The algorithm creates a cubic spline using the SciPy [64] routine interpolate.splrep(), which matches each of the input points along the old grid and is C 2 continuous. The new values are then interpolated from the spline at the new grid points.
3. The equations of motion depend on first and second spatial derivatives, which one can calculate using the differences between neighboring grid values. For a uniform grid, the differentiation stencil is the same for all grid points (excluding the boundaries). For example, to fourth order in ∆x, For a non-uniform grid, a different stencil must be calculated at each grid point. The code finds these stencils to 4th order in ∆x for first derivatives and to 3rd order for second derivatives using the recurrence relations defined in Ref. [65]. At the simulation boundaries the stencils are the same order, but they are no longer centered on the grid points whose derivatives they calculate.
4. The entire computational domain naturally separates into regions based on their local resolution. These regions will share a common time-step ∆N and will be integrated in time together (as per point 6 below), and thus their grid spacing will all be within roughly a factor of two of each other. We additionally require that each region contain at least 40 grid points. 5 The boundaries of these regions are purely an artifact of the AMR itself and are called AMR boundaries. To enable proper time integration at such boundaries, the regions are extended spatially to overlap with neighboring regions. 6 This overlap is similar in spirit to the 'tapered boundaries' introduced in Ref. [63] and avoids the generally lower-order temporal interpolation of the original Berger and Oliger scheme [62].
5. Each region evolves independently for at least one time-step using a fourth-order Runge-Kutta integration. The time-step of the most finely-spaced region is set by the Courant-Friedrichs-Lewy (CFL) condition [66], which demands that the time-step be smaller than the amount of time it takes information to pass between neighboring grid points. That is, c∆N ∆x (we use c∆N ≈ 0.2∆x). Time-steps in more coarsely-spaced regions are set to be powers of 2 larger than that in the most finely-spaced region (so that such regions will be periodically time-aligned). Because of the overlap, we need not worry about boundary effects at the edge of each region: no information can propagate from outside the overlap into a region's interior before a single time-step has completed. At the end of each time-step in each region, the overlapping points are replaced with their corresponding values in the interior of each regions' neighbors, as long as the neighbors have been integrated up to the same time. Since the time-steps come in powers of two, a fine region must integrate two steps for every one step in a coarser neighboring region before the two will match. We set an upper bound on the time-step of the coarsest region to be no more than 2 6 times greater than the time-step in the finest region. This, along with the above-defined minimum region width and time-step size, guarantees that features in high resolution regions will not evolve into low density regions before the coarsest-gridded region has taken a single time-step. 6. After each region has made at least one time-step (the finest-spaced region may make as many as 2 6 steps so that all regions can be matched), the code loops back to step 1 and recreates the non-uniform grid. At specified values of N , the entire grid can be saved to file for later retrieval.  Figure 5. Bubble wall profile at N = 6 at the full simulation resolution. Each dot is one of the discrete points in the simulation itself. The wall is highly length-contracted, so the density of points needs to be extremely large in this region. The region shown here is roughly one hundred thousand times smaller than the size of the full simulation.
In addition to the CFL condition, the integration time-step is limited by the frequency of oscillations about minima in the potential. This will always become the limiting factor near the end of inflation, since α increases roughly exponentially when N is large and V (φ) = 0, We require that the timestep be no larger than 1 30 of the oscillatory period about the inflationary minimum. The entire simulation is written in C with python wrapper functions for easy interactive analysis and for communication with the CosmoTransitions package.

Universes collide
Here we present the results for a single simulation of two bubbles colliding in the 'L2' potential with an initial separation of ∆x sep = 1.0. We run the simulation in two steps. In the first step, the boundaries of the simulation are outside the bubble walls, and at each time-step we enforce the boundary conditions φ = Π = 0; α = a = 1. The simulation region grows along with the bubbles 7 so that it does not needlessly integrate the field in the false vacuum (see Fig. 7). As the bubbles grow, their walls get length-contracted, and more points are added by the adaptive grid algorithm to resolve the wall structure (see Fig. 5). The adaptive grid follows the bubble walls very well, but eventually the extra points noticeably slow down the computation.
In the second step, which we choose to start at N = 7, the simulation region is trimmed to exclude the bubble walls and thereafter held constant. We calculate the spatial derivatives near the boundary using non-centered finite differences. Of course, we should not trust the results of the simulation near the boundary, or anywhere in the boundary's causal future, but luckily this region is very small. For constant α and a, null geodesics obey Since α/a ∼ 1 near the region's boundary at N = 7, the maximum width of the boundary's future light cone is δx ≈ 0.002, which is insignificant in comparison with the total size of the simulation region. By simulating all the way to end of the inflation, we can explicitly check the assumption that the perturbations are frozen in.
The entire simulation takes about 7 minutes to complete on a single core of a modern 2.66 GHz desktop computer, about 6 minutes of which is spent evolving the bubble walls up to N = 7. Fig. 6 shows the initial field configuration and the immediate aftermath of the collision. At N = 0, α = a = 1 and the two instantons are stationary. They then begin to grow, colliding around N ≈ 0.5. The centers of each bubble evolve very little during this timethey are very close to the inflection point φ I , so they do not move much along the potential. After the collision, the two bubbles essentially superimpose -the 'free passage' approximation [21][22][23][24]. The sharp edges of the superimposed bubble walls relax somewhat since there is no longer a strong pressure (potential energy) gradient supporting them. Since the field in the collision region is farther away from φ I , it begins to slowly roll along the inflationary potential. In contrast, the field outside the collision region does not move significantly away from φ I for several more e-folds. Fig. 7 shows contours of the field and the metric functions for N < 5.5. Again, we see that the field in the collision region starts rolling down the inflationary potential much sooner than the field outside the collision region. The metric function a quickly approaches a stable configuration, while α continues to increase with increasing N . This matches the results from Ref. [27]. Fig. 8 shows the behavior near the end of inflation. At N ≈ 50, the field φ reaches the inflationary minimum in the collision region and begins oscillating coherently inside each bubble. The metric function α starts to grow exponentially, while a decreases (although a cosh N continues to increase, so the universe continues to expand). As noted above, the frequency of oscillations about the minimum becomes very large in this region, while the oscillations themselves die down in amplitude. Since lines of constant φ do not correspond to constant N , high frequency temporal oscillations also imply high frequency spatial oscillations. The density of grid points does not significantly increase during this period (the density of points is set to be proportional to |dφ/dx|, which remains relatively small), so eventually the frequency of oscillations becomes smaller than the grid spacing and information about the oscillations is lost. The quality of the solution degrades at this point. Nevertheless, for our purposes we can reconstruct the perturbed metric well within the regime of accurate solutions.
We explicitly check the convergence of the simulation by varying the grid spacing parameters and, implicitly, the integration time-step. To uniformly double the density of grid points across the simulation, we need to double the input parameters n wall , m min , n 2 , the minimum number of grid points per region, and the minimum number of time-steps per oscillation about the inflationary minimum. The Runge-Kutta integration has an error of order ∆N 4 step ∝ ∆x 4 , and the first spatial derivative has an error of order ∆x 4 . The second derivative formally has an error of order ∆x 3 , but by symmetry this leading term is zero for uniform grids and the error becomes order ∆x 4 . Since the grid is increasingly uniform for smaller grid spacing, the error is effectively ∆x 4 . Fig. 9 shows the error associated with the   field φ and the metric functions α and a as a function of time N . The error at any given grid point is estimated as the value at that point minus the value at the same point in a simulation with half the grid spacing. Solid lines show the root-mean-squared error across individual time slices, and dashed lines show the median absolute error. For our chosen input parameters, the error tends to be smaller inside the collision region, which is better represented by the median error estimate. The error in both φ and α converges as ∆x 4 (that is, the error increases by a factor of ∼ 16 when the grid spacing doubles), as expected. The error in the metric function a only converges as ∆x 2 , but with a seemingly subleading contribution. As expected, the error is large after the end of inflation when the spatial oscillation frequency is comparable to the grid spacing. We also check that the constraint equation (Eq. 3.11) holds (see Fig. 10). The constraint holds extremely well (to about 1 part in 10 7 in relative terms) until inflation ends and the field begins rapid oscillation, after which the simulation (expectedly) performs increasingly poorly in satisfying the constraint. Together with the convergence tests described above, this indicates that our code very accurately models the collision spacetime. Both the convergence properties and the behavior of the constraint residuals in this code represent very significant improvements over the code used in Ref. [27].

Finding the Observer Metric
To simulate the collision between two bubbles, we choose coordinates that cover the entire collision spacetime. These are clearly not the coordinates appropriate for describing a nearly homogeneous and isotropic cosmology inside the colliding bubbles. It is possible to describe cosmology in terms of manifestly geometric and covariant quantities [67]. This was the approach in Ref. [38]. For extracting signatures of bubble collisions, it is most convenient to describe the bubble interior in terms of perturbed open FRW coordinates. The challenge in this case is to define a map between such a set of coordinates and the coordinates used in the simulation. In general, this is an ill-defined problem -it is not easy to find a coordinate transformation connecting two generic metrics.
The key insight is to note that as the proper time goes to zero, an open FRW universe can be very well approximated by the Milne slicing of Minkowski space. In the Milne patch, position is labeled by the initial velocity of geodesics emanating from the origin of spherical coordinates. This suggests a prescription for reconstructing the synchronous gauge metric inside bubbles. First, follow a set of geodesics from the bubble center, each labeled by their initial velocity, and track the proper time along each by evolving the geodesic equations through the simulation. Then identify a new set of coordinates associated with the geodesic labels and the proper time, and explicitly construct the map to the simulation coordinates. For an unperturbed bubble, transforming the metric to the geodesic coordinates yields the metric on an open FRW universe. For a perturbed bubble, this procedure yields a perturbed open FRW metric in the synchronous gauge. The synchronous gauge is in general a continuous family of coordinate systems, leading to an ambiguity. However, in our case this ambiguity is fixed by the geometrical prescription used to construct the metric. In effect, the gauge is anchored by the requirement that geodesics probing the unperturbed region of the bubble must describe precisely an unperturbed open FRW universe.
In the following subsections, we build the metric for arbitrary observers in the observation bubble. In Sec. 4.1, we describe two coordinate systems in terms of geodesic trajectories: a Cartesian coordinate system (τ, X, Y, Z), and an anisotropic hyperbolic coordinate system (τ, ξ, ρ, ϕ). The former is useful for describing individual observers and perturbations, while the latter is convenient for integrating a grid of geodesics. In both cases, the time coordinate τ labels the proper time along the geodesics. In Sec. 4.2 we describe the geodesic integration, which supplies the map from simulation coordinates (N, x, χ, ϕ) to anisotropic hyperbolic coordinates, and we show the integrated results. Importantly, N and x are functions of τ and ξ only, so we only need to vary the initial trajectory ξ to recover all features of all geodesics in the (3 + 1)-dimensional space. We would then like to find observer-centered Cartesian coordinates -that is, coordinates in which the observer resides at X = Y = Z = 0 -and to do this we must perform a change in frame which takes the observer to the origin, which we describe in Sec. 4.3. Since the coordinate systems describe an open universe, the framechanging coordinate transformation is more complicated than the Euclidean translation that one would find in flat space. In Sec. 4.4, we describe how to transform from the simulation metric to the observer metric using this multi-step map between the simulation coordinates and the observer-centered Cartesian coordinates. Table 1 summarizes the coordinate systems discussed below.

Coordinate Systems for the Observation Bubble
An important first step is to identify a convenient coordinate system on the bubble interior in the absence of a collision. An unperturbed single bubble contains an infinite open universe, a fact dictated by the SO(3, 1) symmetry of the instanton. We utilize two coordinate systems on the open FRW universe, each of which manifests different properties of the spacetime.
is the (dimensionful) scale factor, and the coordinates take values in the ranges 0 < τ < ∞, 0 < R < 2. For ease of comparison with the simulation coordinates, the proper time τ and the scale factor a(τ ) are both measured in terms of the false-vacuum Hubble scale. The coordinates X, Y and Z are all unitless. If tunneling occurs directly to an inflationary plateau inside the bubble, the scale factor is given by a(τ ) = H F H I sinh H I H F τ , and the comoving curvature radius is R 0 = H F H I . For R 1, the Cartesian coordinates approach those of a flat universe. The Cartesian foliation of an open universe is sometimes referred to as the Poincaré disk.
As τ → 0, for any non-singular configuration we have a(τ ) → τ , and the metric approaches the Milne slicing of Minkowski space: In terms of the coordinates on Minkowski space t, x, y, z, we have 8 and X = R cos θ, Y = R sin θ cos φ, Z = R sin θ sin φ, we obtain t = τ cosh η, x = τ sinh η cos θ, y = τ sinh η sin θ cos φ, z = τ sinh η sin θ sin φ. (4.5) Therefore, we see that the Cartesian coordinates in an open universe can be assigned a very physical interpretation. The coordinates can be thought of as labels attached to a congruence of geodesics emanating from the origin of Minkowski space: τ is the proper time along each geodesic, η is the initial rapidity (the radial velocity and rapidity are related by v = tanh η) of each geodesic, and θ and ϕ parameterize the initial direction of each geodesic. Therefore, we could build up the metric inside an open FRW universe by tracking a set of labeled geodesics. This is the basic method we employ below; however, there is a slightly more convenient coordinate system, for which the symmetries of the collision spacetime are manifest.

Anisotropic Hyperbolic Coordinates.
To make direct connection with the symmetries of the collision spacetime, we define a new set of coordinates. In terms of the open universe Cartesian coordinates, these are defined by 8) or equivalently, These coordinates take values in the ranges −∞ < ξ < ∞, −∞ < ρ < ∞, 0 < ϕ < 2π, and cover the entirety of the bubble interior. The metric is As τ → 0, again we have a(τ ) → τ , and we can see that this metric is an anisotropic slicing of Minkowski space, In terms of the Cartesian coordinates on Minkowski space, we have t = τ cosh ξ cosh ρ, x = τ sinh ξ, y = τ cosh ξ sinh ρ cos ϕ, z = τ cosh ξ sinh ρ sin ϕ.

Geodesic Evolution
As discussed, the coordinate systems for the observation bubble defined above have a clear geometric interpretation as a congruence of labeled geodesics emanating from the center of the bubble at τ = 0. Since the collision will never affect the point where the geodesics begin, we can straightforwardly extend any of the coordinate systems defined above to describe the interior of the perturbed observation bubble. The map between the simulation coordinates and the coordinates foliating the observation bubble can then be obtained by evolving the geodesic equation for each of the labeled geodesics. Because they explicitly manifest the symmetry of the collision spacetime, we employ the anisotropic hyperbolic coordinates. The simulation metric functions are independent of ρ and ϕ, and therefore geodesic motion in these directions is trivial. Evolving the geodesics produces two maps: N (τ, ξ) and x(τ, ξ).
The metric for the simulation coordinates (Eq. 3.4) in the limit where N → 0 is given by Comparing this with Eq. 4.11, the coordinate map between the two metrics is This provides a map at early times, allowing one to label geodesics by their initial conditions. From Eq. 4.14, the initial conditions for geodesics are The relevant geodesic equations are where the Christoffel symbols in the simulation coordinates are The geodesic equations for ρ(τ ) and ϕ(τ ) are always satisfied for ρ(τ ) = const. and ϕ(τ ) = const. in both the unperturbed and perturbed regions -the collision is symmetric in these coordinates. Integrating the geodesics is fairly simple given α(N, x) and a(N, x) -they are just ordinary differential equations -but finding α(N, x) and a(N, x) requires some extra work. The geodesic integration and simulation integration are logically separate procedures, so we perform them in separate steps. 9 We save the Christoffel symbols to a table during the simulation, and then interpolate from the table to get each symbol Γ at arbitrary N and x. We additionally save ∂ N Γ, ∂ x Γ and ∂ N ∂ x Γ to file so that we can use bicubic rather than bilinear interpolation. Bicubic interpolation is not well-defined for an arbitrary non-uniform grid, but the Christoffel symbols are saved in constant-N slices. We first interpolate along a single slice N i using cubic interpolation which produces functions Γ(N i , x) and ∂ N Γ(N i , x) that are C 1 -continuous in x. We can then interpolate between adjacent slices for any given value of x. The resultant function is continuous and smooth. For straightforward post-processing, it is important that we choose the resolution of output slices carefully: too few slices will not resolve important features of the collision front, while having too many slices results in impractically large files. We choose the recording such that each geodesic crosses approximately 20 saved slices in the time that it takes each to cross the collision front. We check that the geodesic integration converges for higher output resolutions, and that the error associated with the chosen resolution is much smaller than the integrated result (see Fig. 11). It can be shown analytically that sinh N ∝ cosh ξ for the single bubble simulation, which matches our numerical results.  For each different bubble collision simulation, we make a grid of geodesics, each starting with a different trajectory ξ. By interpolating along this grid (using linear interpolation, for simplicity), we can find N (ξ, τ ), x(ξ, τ ), and their ξ and τ derivatives at arbitrary points.

Frames
In an unperturbed observation bubble, the surfaces of constant τ inside the bubble are homogeneous -there is no preferred position. However, the collision induces anisotropies inside the observation bubble, making different positions physically inequivalent. The simulated collisions described above are performed in a reference frame where the bubbles both nucleate simultaneously at N = 0. This is the so-called Collision Frame. According to the geometric interpretation of the open universe Cartesian or anisotropic hyperbolic coordinates discussed above, at early times observers at different positions can be thought of as living on boosted trajectories in Minkowski space. To make contact with observables, it is convenient to transform to the so-called Observation Frame, where a hypothetical observer is at rest, or equivalently, at the origin of the anisotropic hyperbolic and Cartesian coordinates. This amounts to performing Lorentz transformations on the coordinates at early times, or equivalently, performing a re-labeling of the geodesics running through the simulation.
To see how this works, consider the Minkowski embedding of the anisotropic hyperbolic coordinates, Eq. 4.12. A boost along the y-direction with rapidity ρ obs has the following action cosh ρ = cosh ρ obs cosh ρ − sinh ρ obs sinh ρ cos ϕ, Along ϕ = 0, using the identity cosh ρ obs cosh ρ−sinh ρ obs sinh ρ = cosh (ρ − ρ obs ), we see that this boost corresponds to a translation of a point ρ = ρ obs to ρ = 0. Importantly, the boost only mixes the ρ and ϕ coordinates in such a way that dρ 2 + sinh 2 ρ dϕ 2 → dρ 2 + sinh 2 ρ dϕ 2 . Similarly, a boost in the z-direction only mixes the ρ and ϕ coordinates. In the collision spacetime, since the metric functions a and α are independent of ρ and ϕ, boosts in the y-or z-direction do nothing to the mappings N (τ, ξ) and x(τ, ξ). This is the underlying SO(2,1) symmetry of the collision spacetime. Now, consider a boost along the x-direction in the Minkowski embedding with rapidity ξ obs . This has the action sinh ξ = cosh ξ obs sinh ξ − sinh ξ obs cosh ξ cosh ρ, (4.28) cosh ρ = cosh ξ obs cosh ξ cosh ρ − sinh ξ obs sinh ξ cosh ξ , (4.29) Along ρ = 0, using the identity mentioned above, we see that this boost corresponds to a translation of a point at ξ = ξ obs to ξ = 0. Boosts in the x-direction alter the mapping between the simulation and cosmological coordinates, taking N (τ, ξ) → N (τ, ξ , ρ ) and x(τ, ξ) → x(τ, ξ , ρ ).
To make contact with observables, we ultimately want to go from anisotropic hyperbolic coordinates to Cartesian coordinates, which will allow us to calculate the comoving curvature perturbation. This can be accomplished by first going to the Observation Frame, and then using Eq. 4.9 to find the mapping between the simulation coordinates and the Cartesian coordinates in the Observation Frame. In Fig. 13, we show the mapping from the anisotropic hyperbolic coordinates in the Collision Frame to the Cartesian coordinates in the Observation Frame. The field and metric functions are, by symmetry, constant along lines of constant ξ in the Collision Frame. In the Observation Frame, these surfaces of constant field and metric functions get distorted by the coordinate transformation. Schematically, Fig. 13 depicts what the general structure of the collision spacetime will look like in the Observation Frame.
The final result of the steps described above is a map from the simulation coordinates to the Observation Frame in Cartesian coordinates for a specific observer located at {ρ obs , ξ obs } in the Collision Frame. Schematically, we have where we have made explicit the fact that there is no physical difference between observers at different ρ obs . The map {N (X, Y, Z), x(X, Y, Z)} can be created for arbitrary observers, allowing us to describe the collision spacetime from an arbitrary vantage point.

Metric Reconstruction
The metric inside the perturbed observation bubble is explicitly obtained from the coordinate transformations discussed above. Using the standard tensor transformation law, we have where x α represent the simulation coordinates, X µ represent the Cartesian coordinates in the Observation Frame, and g αβ [x(X)] is the metric from the simulation. The resulting metric inside the observation bubble, g µν [X], is by construction in the synchronous gauge: g τ τ = −1 and g τ i = 0. In practice, we perform the various coordinate transformations in sequence. We first find the metric in terms of the anisotropic hyperbolic coordinates, which are the direct products of the geodesic evolution. The accuracy of the metric reconstruction is dependent on how well the geodesics reproduce the coordinate map and its first derivatives. We have explicitly shown convergence in both the simulation and geodesic integration, and thus we expect an accurate reconstruction of the metric. We then transform the metric in the Collision Frame (anisotropic hyperbolic coordinates) to the Observation Frame (described in Cartesian coordinates) for an ensemble of observers at different values of ξ obs .

The Comoving Curvature Perturbation
The metric obtained in the last section describes the interior of a bubble in the presence of a collision. However, a full description of the bubble interior is not necessarily relevant for constructing cosmological observables. Because we observe a nearly homogeneous and isotropic universe, any relevant signatures arise from small (linear) deviations from FRW. By mapping the full metric onto a metric describing a linear perturbation of open FRW, we can take advantage of the extensive results of cosmological perturbation theory (see e.g., Ref. [68]).
In making this map, we must choose a gauge (coordinate system) for the perturbed FRW metric. The metric we obtain from the geodesics is explicitly in the synchronous gauge. This metric can, in principle, be used to extract any observable quantity of interest. However, transforming to the comoving gauge facilitates a more direct connection with observables, which are determined entirely by the comoving curvature perturbation. In regions where the metric is nearly FRW, this is accomplished by a linear gauge transformation, and, in linear theory, the comoving curvature perturbation is conserved on super-horizon scales. This allows us to directly connect the output of our simulations to cosmological observables without considering the details of reheating. 10 The rest of this section describes how to calculate the comoving curvature perturbation R and its behavior for different initial bubble separations ∆x sep . A summary of the necessary steps is presented below. After reading the summary, the more results-oriented reader may wish to skip directly to Sec. 5.4.
We start in Sec. 5.1 by mapping the full metric into scalar perturbations on a background FRW metric in the synchronous gauge. We write the metric perturbation (the fractional difference between metrics in perturbed and unperturbed bubbles) as , where E ij is traceless and symmetric. The D-term is a purely scalar conformal perturbation, whereas E ij contains scalar, vector, and tensor components. To make the extraction of E from E ij tractable, we assume that the tensor and vector components are negligable and that the perturbation is approximately planar (a function of the coordinate X only). We justify these assumptions in Sec. 5.1. For the particular Lagrangian under study, we later find that E ij is small and therefore the exact planar assumption is not necessary to calculate the scalar perturbations. (We have not verified that this is true for other Lagrangians, so we present a way to calculate E here.) Together, D and E represent scalar perturbations on an open FRW metric. To convert to perturbations on a flat metric, we assume that the region of interest satisfies R = √ X 2 + Y 2 + Z 2 1, allowing us to write the overall curvature of the open bubble universe as an additional small scalar perturbation D (curv) . Note that unlike perturbations in flat FRW, the scalar perturbations have coordinatedependent contributions: D (curv) goes to zero at the origin of coordinates and can be large far from the origin. It is most convenient to place the observer at the origin, so picking a point to be the origin amounts to choosing the Observation Frame. The scalar perturbations D (syn) , D (curv) , and E are then defined in that frame.
Collectively the perturbations should satisfy the linearized Einstein equations; we check that this is the case in Sec. 5.2. We then describe how to switch to comoving gauge and retrieve the comoving curvature perturbation R in Sec. 5.3.
We describe our numerical results in Sec. 5.4, verifying that R is time-independent for all initial bubble separations ∆x sep when the linear Einstein equations are satisfied. We find that E = 0 exactly, one can show that D (syn) depends only upon the map N (ξ, τ ) from anisotropic hyperbolic coordinates to simulation coordinates. If we drop the contribution from D (curv) (which is constant on any sphere surrounding the observer, and constant along the surface of last scattering), the perturbation R then becomes explicitly independent of the choice of coordinate system: transforming R from the Collision Frame to an Observation Frame amounts to a straightforward change in the underlying coordinates, and we can treat R as a function of anisotropic hyperbolic coordinates ξ and τ only. When R is time-independent ('frozen out'), it is a function of ξ only. This is an important convenience, allowing us to later determine the CMB signature for arbitrary observers in Sec. 6 in a straightforward manner. Our working assumptions in Sec. 5.4 and onwards are then: E (syn) ij = 0 and δg ij ∝ δ ij (so that R is explicitly coordinate-independent); and R 1 (so that the overall curvature perturbation can be treated linearly). Note that we no longer need to assume that the perturbation is exactly planar, since that was only needed to calculate E.
Finally, in Sec. 5.5, we find that R(ξ) is well described by a power law near the collision boundary. Simulating a number of collisions, we show the power law index κ to vary from κ ≈ 3.5 for small ∆x sep to κ ≈ 2 for large ∆x sep .

Perturbations in the synchronous gauge
With the metric Eq. 4.31 from the geodesic simulations in hand, we can compare with standard results of cosmological perturbation theory to make contact with cosmological observables. The most general perturbed metric in synchronous gauge can be written as When there is no collision, D = E = 0, and we recover the open FRW metric, Eq. 4.1.
We find the components of the perturbed metric Eq. 5.1 as follows. First, we perform a simulation of a single bubble to define a reference background geometry. We then use the geodesic metric reconstruction method described above to find the components of the metric in terms of Cartesian coordinates, g (no coll) ij . This matches with Eq. 4.1 to high accuracy (we quantify this below), providing an important check of our geodesic reconstruction method. Next, we perform a simulation of a collision, and re-run the geodesic reconstruction routine to obtain the metric functions g (coll) ij inside the perturbed observation bubble. Comparing with Eq. 5.1, we can identify Taking the trace, we find The components of E ij can be written as

vector and tensor components). (5.6)
This is traceless and symmetric. The tensor contribution is zero by a hyperbolic version of Birkhoff's theorem [14]. We find the vector component, which decays during inflation, to be small empirically. Therefore, we consider only the scalar component in what follows. We will later need to separate out the various partial derivatives in the scalar part of E (syn) ij given the components of the perturbed metric Eq. 5.3. Equating Eq. 5.5 and Eq. 5.6 (which holds up to the neglected vector contribution), we obtain a set of elliptic equations for E (syn) , This set of equations must be supplemented by boundary conditions. As we will see shortly, an understanding of the correct boundary conditions is not strictly necessary. Two essential simplifications result when we consider cosmologies where curvature is not dominant today, implying R 1 (where R is the Cartesian radius). First, we can treat the spatial curvature as a small perturbation in the region of interest, In linear theory, we can treat this separately from the perturbations sourced by the collision. Second, since the collision is a function of ξ only, and close to R = 0 we have X ≈ ξ − ξ obs , we can assume that E (syn) = E (syn) (X). That is, we assume that the perturbation E is exactly planar. We can then write The elliptical equation for E becomes an ordinary second-order differential equation in X, with the most natural initial condition being E = dE dX = 0 in the unperturbed region. However, we need not solve this equation, since only the derivatives of E are necessary to calculate the comoving curvature perturbation. We will later show that E ij D, so the precise calculation of E and its derivatives has little effect upon our results.
We will find it convenient to define a new variable Ψ, the curvature perturbation in synchronous gauge Making this change, the metric is

Checking the perturbative description
With these definitions, we can verify our perturbative description. First, we check that the components of E ij satisfy Eq. 5.9. These equations imply that E is very nearly planar, an important assumption in our analysis below. Indeed, in the vicinity of R = 0, we find that the planar assumption works well, as described below. Another important test is to check that the metric in synchronous gauge obeys the Einstein equations. To zeroth order in perturbations, the metric functions should satisfy the Friedmann equations where φ 0 is the field with no collision. At first order in perturbations, we get a constraint from the τ -X component of the Einstein equation: where δφ = φ (syn) − φ 0 , and we work along the X-axis such that Y = Z = 0. Going to second order, we have (5.14) Substituting Eq. 5.10 into Eq. 5.14 yields We expect this constraint to hold near the collision boundary where the perturbation is relatively small and the use of perturbation theory is valid. We can explicitly check it by computing the perturbation variables using the procedure outlined above. The result is shown in Fig. 14. The constraint holds quite well at the edge of the collision, but does not match arbitrarily far into the collision region. It can also be seen that the second-order correction yields a result which matches further into the collision region than when including the first order correction only. This is the expected behavior -perturbation theory is working, but higher order terms become more important as the magnitude of R increases deep inside the collision region. The collision can be treated as a small perturbation up to ∆X ≈ 0.15 beyond the collision boundary for ∆x sep = 0.5, and up to ∆X ≈ 0.05 beyond the boundary for ∆x sep = 2.5. This follows a general pattern, which we will explore in more detail below, of smaller kinematic separations between bubbles yielding smoother and wider collision boundaries.
We conclude that -sufficiently close to the causal boundary of the collision -the perturbative description, and the assumption of planarity, are accurate. This is also an excellent test of our code, since numerical errors would also lead to violations of the constraint equation.

Synchronous gauge to comoving gauge
In order to facilitate comparison with observations, we would like to go to the comoving gauge, in which the metric is (5.16) and the energy momentum tensor has the property that T (com) i0 = 0. One can show that metric functions in comoving gauge are related to those in synchronous gauge by where x µ (com) = x µ (syn) + µ (x (syn) ) and µ (x (syn) ) = 0 , ∂ i . See e.g., Ref. [68] for a discussion of linear gauge transformations in cosmological perturbation theory.
The three-curvature on spatial slices in the comoving gauge is given by This defines the comoving curvature perturbation (dotted lines). The evolution of the individual components almost exactly cancel everywhere that the linearized Einstein equations hold, indicating that the perturbation has 'frozen in.' Substituting Eq. 5.10, we can write this as where we have now dropped the curvature perturbation, which is a small constant on any sphere of fixed distance from the origin (relevant for the projection discussed in the next section).

Behavior of the comoving curvature perturbation
With the above machinery we can calculate the comoving perturbation. Fig. 15 shows the shape and evolution of the perturbation for two different collisions as a function of X. In both cases, the perturbation smoothly increases from R = 0 -there is no sharp edge to the collision boundary. Comparing the two collisions, a higher value of ∆x sep yields a faster and sharper rise in the perturbation. Both the metric perturbation D (syn) and the Hubble rate H evolve in time, 11 but near the collision boundary the two contributions almost exactly cancel and R is roughly constant. In perturbation theory, for adiabatic fluctuations, the comoving curvature perturbation should freeze in as exhibited in Fig. 15. Away from the collision boundary, the perturbation exhibits a small but significant time-dependence. The magnitude of the dependence is roughly proportional to the error in the first-order piece of the metric constraint equation (see Fig. 14). The perturbative freeze-in is a result of the linear approximation, so the perturbation 'thaws' to the extent that the linear approximation fails. This is an important consistency check, indicating that we are accurately constructing the curvature perturbation. In Fig. 16, we show the collision in the X − Y plane near the end of inflation. Because the comoving curvature perturbation near the causal boundary of the collision is frozen in, this snapshot of R will be representative of R after inflation. Recall that R is defined in the Observation Frame, where the hypothetical observer is at the origin of coordinates. This observer will have causal access to some portion of this early-time hypersurface, as depicted in Fig. 16. The precise projection of the observer's past light cone onto this early-time surface depends on the details of the cosmology, and is discussed below. However, in the limit where the observed curvature is small, the causally-relevant portion of the collision will be a small region around the origin. As expected, R has an approximate planar symmetry (lines of constant R are very nearly lines of constant X) on scales much smaller than the curvature radius in the open FRW.
We expect the comoving perturbation R to align with the symmetries of the simulation. That is, we expect that at fixed τ , R is a function of ξ only. We can explore the perturbation from the vantage of any observer by performing the mapping described in Sec. 4.3 (see in particular Fig. 13; lines of constant ξ correspond to lines of constant R). The benefit of working with ξ is that, in the limit where E is small, going to the Observation Frame for an arbitrary observer amounts to a translation in ξ. Setting E to zero, there is only a conformal perturbation to the spatial part of the metric, making the coordinate transformation from Cartesian to anisotropic hyperbolic coordinates trivial: where dH 2 2 = dρ 2 + sinh 2 ρ dϕ 2 , and we restored the factor of 1 − R 2 4 −2 (the appropriate non-linear completion of D (curv) ). Moving from ξ obs = 0 to another Observation Frame with ξ obs then amounts to taking That is, we simply change coordinates in R. Lines of constant R in the Observation Frame will appear much as in Fig. 13. If we are interested in a region close to the observer's origin, we can approximate which is a simple translation in ξ.
If E = 0 exactly, we can also re-cast the perturbed metric in terms of variations in N . In this case, g xx = g yy = g zz , and, at ϕ = 0: √ g zz = sinh N 1 + cosh ξ obs cosh ξ cosh ρ − sinh ξ obs sinh ξ 2 cosh ξ (5. 23) or, for ρ = 0, √ g zz = sinh N 1 + cosh(ξ − ξ obs ) 2 cosh ξ . (5.24) Therefore, if we consider the spatial variation in N as a perturbation δN on top of some background mean value N 0 , for large N 0 we have 26) and the entire perturbation is a function of ξ, τ and ∆x sep only. This is also the behavior seen in Fig. 16: lines of constant R are very nearly exactly lines of constant ξ. Note that the coordinate N roughly measures the number of e-folds, so Eq. 5.25 is reminiscent of the δN formalism [69]. Fig. 17 shows the behavior of the perturbation for varying ∆x sep . When the initial separation between bubbles is smaller, the bubbles collide at earlier times and the collision region penetrates farther into the observation bubble (bluer lines have smaller ∆x sep ). The amount of energy in the bubble walls is also smaller since they have had less time to accelerate, so the perturbation is not as steep in the immediate vicinity of the collision boundary. Far Each line shows the perturbations at τ = 60 for a different value of the kinematic separation ∆x sep , ranging from ∆x sep = 0.4 to ∆x sep = 3.0. As can be seen in Fig. 2, the position of the collision boundary depends on the kinematics. We also calculate the perturbation E xx (right plot), finding E xx R everywhere.
inside the collision region, the perturbation is roughly the same for different kinematics. However, far inside the collision region, linear perturbation theory is not reliable, so R does not have a clear physical meaning. We find that the metric perturbation is almost entirely proportional to the identity matrix. That is, D (syn) E ij . Fig. 17 shows that E xx contributes no more than 0.25% to R for ξ < 3, and much less than that for ξ < 1. Interestingly, E xx seems to be roughly proportional to ∂(x (coll) − x (no coll) )/∂ξ (see Fig. 12). At the observer's origin, E yy = E zz = − 1 2 E xx exactly, and off-diagonal components disappear due to symmetry in Y and Z. This is no longer true for Y, Z = 0, but we find that all components of E ij remain small at all coordinates (with R 0.3) for all observers (boosted by ξ obs 3).
As discussed above, when E can be safely neglected the perturbation from the perspective of different observers can be constructed by a straightforward transformation in the underlying coordinates (see Sec. 4.3). For Y = Z = 0 (or equivalently, ρ = ρ = 0), this corresponds to a simple translation in ξ. In the sections that follow, we will generally assume that R is a function of ξ only. That is, we will assume that E ij = 0 exactly and that the perturbation is frozen out so that it has no time dependence. Under these assumptions, the left panel in Fig. 17 shows the complete comoving perturbation for all observers (with ξ obs 3): different observers just transform ξ differently to their local Cartesian coordinates. For R 1, the transformation is X ≈ ξ − ξ obs , effectively translating each R(ξ) curve by a fixed amount.

Fitting the comoving curvature perturbation
To facilitate comparison with observation, we define a template for the shape of the comoving curvature perturbation. We then fit for these coefficients as a function of the initial bubble separation. Very near the causal boundary of the collision, we find that the perturbation is well described by a power law, Θ(ξ − ξ c (∆x)). The coefficients A(∆x), ξ c (∆x), and κ(∆x) are all functions of the initial bubble separation; ξ 0 = 0.05 is a constant characteristic of the collision regions' width and the region over which the linear approximation holds. We find the coefficients by minimizing w i [R fit (ξ i ) − R i ] 2 in some small region surrounding the collision boundary, where the weights w i are chosen to increase the fit's accuracy very close to the boundary where R is small. We use a Nelder-Mead downhill simplex algorithm as implemented in the scipy.optimize package to perform the minimization. Fig. 18 shows some representative fits for different values of ∆x, and Fig. 19 shows the fit parameters as a function of ∆x as well as the range in ξ over which the power law solution is a good approximation.
When the initial separation between the bubbles is very small, the collision region does not have a clear boundary. There is always a well-defined region that is in the causal future of the center of the colliding bubble. However, the exterior of the bubble wall will never actually be in this region, and the wall itself does not have a sharp edge (see Fig. 4). The bubble wall length-contracts as the bubble grows, so the wall sharpens as it asymptotically approaches the well-defined bubble interior. However, for small ∆x sep , the bubbles have not had much time to grow before they collide, and the bubble walls extend significantly beyond the bubble interiors. Because the collision does not have a sharp boundary, a power-law fit cannot, in principle, accurately reproduce the perturbation arbitrarily close to the collision region boundary. This is seen in Fig. 18 at ∆x sep = 0.5, where the perturbation does not go to zero at what one would naively guess to be the collision boundary. The power-law fit does a much better job at larger kinematic separations: there, the non-zero perturbation outside the collision region is primarily due to (small) numerical integration errors. . The lower right plot shows how far into the collision region one can travel before the error between the simulated and fitted perturbations is always greater than 5%.
Most previous work has assumed that the curvature perturbation in the collision region could be considered piecewise continuous, and linear near the boundary. This does not appear to be the case. Near the boundary, we find that a good description of the perturbation requires continuous first derivatives; a power law with an index between 2 and 3 is a good fit. Moreover, the region over which this power law describes the perturbation well is more than an order of magnitude larger than the width of the collision shock as seen in the simulations (making the appropriate transformation from x to ξ). Thus in the thin-wall limit, even if the function could be considered piecewise continuous, it is not clear that it would be linear in ξ (and indeed, as the perturbation in this case would not be an analytic function at the junction point, there is no reason it need be described by a Taylor expansion there.) In this section, we have outlined a mapping procedure between the metric extracted from our simulated bubble collisions (Sec. 4) and the scalar metric fluctuations in comoving gauge. We only considered the case where the causal boundary of the collision is very near the position of a fiducial observer at the origin of coordinates, and focussed on a window of radius R 1 around the origin. Within this small window, we find that the collision is nearly planar, which allows us to calculate E (syn) , a necessary ingredient in the transformation from the synchronous to the comoving gauge. Empirically, we find E (syn) to be small, in which case we can express the perturbation entirely in terms of the comoving curvature perturbation R. In this limit, we can easily translate R to explore the perturbation from different vantage points. (Note that this is a convenience, and not in general necessary.) The comoving curvature perturbation R is well described by the power law (Eq. 5.27, with coefficients shown in Fig. 19) in the vicinity of the causal boundary of the collision. This represents a complete set of predictions given an underlying scalar field Lagrangian, and plays a key role in the following discussion of CMB observables.

CMB Signature
In this section, we compute the observational signatures of bubble collisions in the temperature anisotropies of the CMB. As we saw in the previous section, the comoving curvature perturbation is time-independent in the spatial region over which perturbation theory is valid. For the single-field models we are studying, only adiabatic fluctuations are excited by the collision, and the relevant part of the comoving curvature perturbation will remain timeindependent throughout the duration of inflation. Once reheating occurs, and the standard Big Bang cosmology begins, the fine-grained shape of the template will be processed by the transfer function. The observed temperature anisotropy in the CMB will depend on both the underlying curvature perturbation and the position of the observer.
Centering the collision on the north celestial pole, the temperature anisotropy will depend only on the polar angle θ, and as depicted in Fig. 16 we have Here, X is the position in Cartesian Observer Coordinates of a point on the observer's lightcone, which has radius R ls . The approximation X ξ − ξ obs holds in the small-R ls approximation, which applies when the observed curvature is small (as quantified below), and in which lines of constant ξ correspond to lines of constant X. The location of the causal boundary of the collision, θ c , is then given by where ξ c is the location of the causal boundary in the collision frame and ξ obs is the position of the observer. The comoving curvature perturbation (see Eq. 5.27) is determined entirely by the kinematics, ∆x sep , while the projection is determined entirely by the observer position, ξ obs , and the size of the observer's past light cone on the surface of last-scattering, R ls . This is the minimal set of variables determining the temperature anisotropy caused by a collision. A set of variables which will be more closely tied to the observed signature would include the parameters in our fit for R (as a proxy for ∆x sep ), θ c (as a proxy for ξ obs ), and R ls . Since all the parameters in R are a function of ∆x sep , and κ is monotonic in ∆x sep , we are free to choose κ as a proxy for the template, which then uniquely specifies the other fitting parameters. With this, we have δT 3) The kinematics and observer position are independent of the details of the inner-bubble cosmology; R ls is determined entirely by the cosmology inside an undisturbed bubble. Returning to Fig. 1, the details of R are fixed by region 1 of the potential and the kinematics, while R ls is determined by region 3. For a ΛCDM cosmology [18], we can identify . The full-sky temperature perturbation δT /T as seen by observers with R ls = 0.05 using: R(X, Y ) (solid lines); R(ξ), approximating E = 0 (dashed lines, which are nearly coincident with the solid lines); and R fit (ξ) (dotted lines). We approximate ξ ≈ ξ obs + R ls cos θ, whereas we use the precise coordinate transformations for X = R ls cos θ and Y = R ls sin θ. Different lines show different observer positions ξ obs , or equivalently different bubble sizes θ c . For ∆x sep = 0.5, we subtract off a constant perturbation such that R(θ = π) = 0.
In the limit of instantaneous reheating, the size of the observer's past light cone is directly related to the number of inflationary e-folds: R ls ∝ e −Ne .
To get an idea of the shape of the observable signature in the CMB, we can compute the temperature anisotropy in the Sachs-Wolfe [58] limit, This gives an accurate form for the temperature anisotropies on the largest scales. Using the template for the curvature perturbation, Eq. 5.27, and the Sachs-Wolfe formula, Eq. 6.5, we find The central amplitude of the template (using the terminology of Refs. [10,11,40]) is given by A representative set of profiles from the simulation are shown in Figs. 20 and 21. In these figures, we compare the result for the temperature profile using the full simulation data for R (with and without neglecting the contribution from E) with the fitting function Eq. 5.27. The fitting function works very well for collisions with large θ c , and for large kinematics ∆x. This reflects the increased accuracy of the fit for R for larger values of ∆x.
For the model we have chosen, the signature has a number of noteworthy properties: • For this model, the collision is always a hot spot (δT /T > 0).
• Fixing θ c and R ls , the central amplitude increases with increasing ∆x sep . This makes intuitive sense: larger initial bubble separations give rise to higher energy collisions, yielding a stronger signal. π/2 π/3 π/6 0 ∆x sep = 2.5 Figure 21. The same plots as seen in Fig. 20, but with a logarithmic y-axis. The power-law is not a perfect fit all the way down to the collision boundary, and errors in finding ξ c appear larger when mapped to θ for θ 1 (since dθ/dξ diverges at θ = 0).
• Fixing κ and R ls , the central amplitude and size are not independent parameters. Larger collisions (in terms of θ c ) generally have a larger amplitudes.
• The central amplitude is large in magnitude, and therefore phenomenologically not viable, unless R ls is quite small. Considering large-size collisions (θ c ∼ π/2), for δT /T < 10 −5 (a rough limit from observation), we require R ls < ξ 0 (5 × 10 −5 /A) 1/κ . Kinematics with the largest ∆x sep provide the most stringent constraint, yielding R ls < 10 −3 for the potential we have studied. The energy density in curvature at last-scattering is related to R ls by Ω k = R ls 2 2 , indicating a limit Ω k < 2 × 10 −7 on the curvature in this model. The (in-principle) detectable level of curvature is Ω k ∼ 10 −5 . Therefore, in this case, detectable curvature and detectable (but not-yet-detected) bubble collisions are incompatible. This statement may vary for other models.
Since we have analyzed only one model, it is far from clear that the template parameters found above are representative of a generic single-field potential. Because the form of the curvature perturbation is different from that found in previous studies of the bubble collision signature, our template for the CMB temperature anisotropy is likewise somewhat different. The first theoretical calculation of the CMB signature was performed in Ref. [6], where the collision was treated in the thin-wall limit, and the backreaction of the perturbed inflaton on the geometry was neglected. These authors found, for a particular set of kinematics, a template identical to Eq. 6.6 with κ = 1. In Ref. [70], a power series ansatz was made for the comoving curvature perturbation. With sufficient inflation, and assuming a thin bubble wall, the leading term with κ = 1 was found to dominate. Incorporating a finite wall thickness, it was suggested that the template can change. A mapping from the numerical simulation onto the initial conditions in the open FRW universe would allow a direct comparison between the analytic model of Ref. [70] and the output of the numerical simulations. We leave this to a more systematic study of different models to be performed in future work. However, in Sec. 5.5, we have shown that wall thickness alone cannot account for the shape of the template. These changes have important implications for the detailed polarization signature discussed in Refs. [7,8]; we will return to predictions for polarization in future work. What is clear from our analysis is that the power law index κ can vary significantly with kinematics, meaning there is definitely no universal index.

Collision Prior
What are the detailed predictions of theories with bubble collisions? In the previous section, we saw that any individual collision can be described by the minimal set of parameters θ c , κ, and R ls . However, since bubble nucleation is a random stochastic process, a model of eternal inflation can only make probabilistic predictions for how many collisions, N coll , should appear in the CMB, and for what values the template parameters should take. In addition, there is uncertainty in the model itself, which has implications for observable signatures.
Considering the predictions of all models is beyond the scope of this paper. However, the properties of R are insensitive to: (1) an overall re-scaling β of the scalar field potential in the simulation; and (2) changes in the total number of post-nucleation e-folds of inflation N e . Returning to Fig. 1, this corresponds to considering variations in the overall scale of the potential in regions 1 and 2 with independent variations of the shape of the potential in region 3. Rescaling the potential leads, as we will see below, to changes in the prediction for the total number of observable collisions. Changing the total number of e-folds leads to different values of R ls . We therefore have a restricted class of models in which four fundamental parameters can vary: ∆x sep , ξ obs , N e , and β. We have four observable parameters for which we want predictions: θ c , κ, R ls , and N coll . We hold the details of reheating fixed; allowing this to vary would provide additional sources of variation in R ls and N coll .
Given a set of prior probabilities over the fundamental parameters, the goal of this section is to implement a change of variables, Pr(N coll , R ls , κ, θ c ) = Pr(β, R ls , ∆x sep , ξ obs ) dN coll dR ls dκ dθ c dβ dR ls d∆x sep dξ obs −1 . (7.1) As discussed above, we consider an ensemble of models in which R ls and β are produced by independent variations of the underlying scalar field potential. The kinematics and observer position are independent variables. We assume that these parameters are also independent of the variations in the potential (a good assumption as long as the initial radius of colliding bubbles is a small fraction of the false vacuum cosmological horizon size). The prior over fundamental parameters is therefore separable, Pr(β, R ls , ∆x sep , ξ obs ) = Pr(β)Pr(R ls )Pr(∆x sep )Pr(ξ obs ) , (7.2) allowing us to discuss each piece separately before finally assembling the desired distribution in Eq. 7.1.

Prior over potential re-scalings
First, we discuss the prior over overall re-scalings of the potential, and its relation to the predicted number of collisions in the CMB. The number of collisions whose causal boundary intersects the observable portion of the surface of last-scattering is given by [18,29] where λ is the probability of bubble nucleation per unit four-volume, H F is the Hubble scale in the false vacuum, and H I is the Hubble scale during inflation. This expression is valid in the limit where the collision is a small perturbation on the background cosmology which, as we have shown above, is a valid assumption over some range of parameters. The nucleation probability can be written as In the thin-wall limit, and neglecting gravitational effects, the tunneling exponent is given by where we have parameterized the potential, field, and wall tension by the scales The total number of predicted collisions depends both on the detailed shape of the barrier and the details of the cosmology inside an undisturbed bubble. With this scaling, we also have Note that we have further assumed that the combination AH −4 F does not change under rescalings. The predicted number of collisions is therefore where we have defined which does not change under re-scalings by β. For the models studied in this paper, we have H I ∼ H F , although a large hierarchy between these scales may be essential for obtaining a significant probability for observable bubble collisions. It is not clear what prior probability we should assign to the overall scale of the potential. We can assume a flat prior for illustration, in which case where ∆β is the range in β spanned by our ensemble of models. We consider the fiducial set of parameters: B = 100, 10 < β < 100, and N 0 = 10 5 in the following.

Prior over R ls
The prior over R ls is essentially a prior over the number e-folds of inflation (for fixed reheating temperature), since R ls ∝ e −Ne . A number of proposals have been made for the prior over the number of inflationary e-folds [71][72][73][74][75][76][77]. The results are highly measure-and model-dependent, but most proposals are in agreement that obtaining a large number of e-folds is suppressed.
As an example, consider a power law distribution for the number of e-folds Pr(N e ) ∝ N −m e (such as was found by Refs. [72,73]). Changing variables, we find Pr(R ls ) = Pr(N e ) dN e dR ls ∝ 1 R ls (− log R ls ) m .   Figure 22. A conformal diagram depicting the prior calculation. An observer inside the observation bubble has causal access to a portion of the surface of last-scattering (red line). Bubbles nucleating in the shaded region will produce a collision whose causal boundary projects onto various angular sizes θ c on the CMB sky of the observer. The probability for bubble nucleation is uniform in each parcel of four-volume in the false vacuum.

Prior for kinematics and observer position
The prior over ∆x sep and ξ obs is most easily computed in the Observation Frame, where these variables map onto a nucleation site in the false vacuum. This is depicted in Fig. 22; a detailed treatment of this calculation can be found in Ref. [18]. A particularly convenient foliation of the false vacuum exterior to the observation bubble is given by the metric where −∞ < Υ < ∞ and −∞ < Ψ < ∞. The bubble wall is located at Υ → −∞ in the limit where the wall sits on the light cone, and at finite Υ in the realistic case where the bubble has finite size.
In the Observation Frame (unprimed) a position {Ψ, Υ} in the exterior is mapped to a position {Ψ = 0, Υ = Υ} in the Collision Frame (primed). This will translate the observer in the open FRW universe along x = 0 to the position ξ = −Ψ.
(7.13) Therefore, we see that the Υ coordinate maps onto the kinematics in the Collision Frame, and the Ψ coordinate maps onto the observer's position in the Collision Frame at fixed kinematics. Note that, as described in Sec. 4.3, points not along the x−axis will transform in a non-trivial way under a transformation between the Collision and Observation Frames. The kinematics in the simulation are determined by ∆x sep , which is related to Υ by ∆x sep = 2 arctan tanh Υ 2 + π 2 . (7.14) The differential number of bubbles coming from each volume element in the Observation Frame is given by dN = λH −4 F cosh 2 Ψ cosh 4 Υ sin θ 0 dΨdΥdθ 0 dφ 0 . From Eq. 6.2, this corresponds to collisions that encompass zero angular scale all the way to the full sky. 12 The position of the collision is a function of ∆x sep only. We can estimate it in the thin wall limit as This quantity in practice is directly extracted from simulations: see Fig. 19.

Assembling the prior
Dropping constant factors, the full prior over the parameters describing the ensemble of models and nucleation sites is given by Pr(β, R ls , ∆x sep , ξ obs ) ∝ 1 R 2 ls (− log R ls ) m cosh 2 (ξ obs ) sin 3 (∆x sep ) .

(7.21)
Our final task is to change variables to find the prior over observables, Eq. 7.1. We first evaluate the Jacobian in this expression. From the functional dependence of the observables (N coll (β, R ls ), κ(∆x sep ), and θ c (∆x sep , ξ obs , R ls )) the Jacobian is  In this expression, κ(∆x sep ) and ξ c (∆x sep ) are obtained directly from the simulation (see Fig. 19), and then inverted to find ξ c (κ) and ∆x sep (κ). From Eq. 7.8, the possible values for N coll are constrained by the constants N 0 , B, and β as well as the value of R ls . In Fig. 23 we show the marginalized distribution function Eq. 7.24 for a fiducial set of models with B = 100, 10 < β < 100, N 0 = 10 5 , and m = 3. There are a number notable properties of the marginalized probability distributions. First, looking at the properties of the individual collisions, the distribution over angular scales is very close to Pr ∼ sin θ c , preferring large-scale collisions on the sky. The power-law index in Eq. 6.6 is preferentially near κ = 2.25. Moving to the distribution over N coll , we see that there is significant probability for N coll > 1. The distribution is weakly weighted towards fewer total collisions in the range N coll > 1. For the fiducial model we study, there is equal marginalized probability for the no collision (N coll < 1), few collisions (1 < N coll < 25), and many collisions (N coll > 25) cases. Lowering the lower bound on β gives higher weight to the no-collision case (N coll < 1), and maintains the relative weights between the few-and many-collision cases. Increasing N 0 gives higher weight to the many-collision case, preserving the relative weight of the no-and few-collision cases. Therefore, the allowed ranges for these parameters can play an important role in the nature of the prediction. Finally, the marginalized distribution over R ls gives significant weight to relatively large values of R ls . There is roughly equal weight for R ls in the range that produces amplitudes compatible with the observed CMB (R ls < 10 −3 ) and incompatible with the observed CMB. The observable range might span roughly an order of magnitude in R ls ; there is roughly a 20% chance to find R ls in this interesting region between 10 −4 < R ls < 10 −3 .
The prior probability distributions relevant for observational searches must take into account some of the properties of the experiment itself. For example, a more relevant prior is over the expected number of observable collisions as opposed to the expected total number (see Refs. [10,11,40]). This is a conditional probability, determined by making a cut in the size θ c and amplitude z 0 (see Eq. 6.7) of the signature. These cuts will affect the shape of the marginalized distributions shown in Fig. 23; we leave a detailed discussion of this point for future work.

Conclusions
This paper has provided, for the first time, a direct quantitative set of predictions for the observable signatures of bubble collisions from an ensemble of scalar-field models that give rise to eternal inflation. This is an important proof-of-principle -we have shown that it is feasible to directly compute the detailed predictions of eternal inflation, linking cosmological observations (such as that of the CMB) to fundamental theory. For the ensemble of models we have studied, we have rigorously shown that bubble collisions can be consistent with a mildly inhomogeneous cosmology, and produce interesting, possibly observable signatures. We have also identified a minimal set of phenomenological parameters, and computed the prior probability distribution over them. This constitutes a complete set of predictions for the ensemble of models we have considered.
An important product of this work is a computer code that allowed us to simulate bubble collisions and the entire post-collision inflationary cosmology inside a bubble. We have extracted a map from the simulation to the standard variable for cosmological perturbation theory, the comoving curvature perturbation R, which allows us to make direct contact with observables. For a comprehensive roadmap of the paper, and an overview of how this was accomplished, see Sec. 2.
Our results corroborate most of the qualitative features of previous analytic, numerical, and semi-analytic studies of bubble collisions. There are nevertheless some quantitative differences. The detailed form of the curvature perturbation depends on the kinematics of the collision, and we have been able to simulate the entire kinematic range. The curvature perturbation is well described by a power law (of index between 2 and 3) in the region to the future of the collision for the ensemble of models we have studied. This leads to an azimuthally-symmetric CMB temperature anisotropy that falls off with angular distance from the center. In the Sachs-Wolfe limit, the temperature anisotropy follows a power-law in the cosine of the angle from the center. This fall-off is faster than was found in previous studies, which assumed that the curvature perturbation was linear in the region to the future of the collision.
We have identified four phenomenological parameters describing bubble collisions in this model: the total number of predicted collisions (N coll ), the radius of our past light-cone at last-scattering (R ls ), and the angular size (θ c ) and detailed shape (specified by κ) of each collision. We derived the prior probability over these phenomenological parameters for an ensemble of models by considering the distribution of allowed kinematics and observer positions, and by varying the overall scale of the underlying scalar-field potential as well as the number of e-folds of inflation inside the bubble. This prior probability prefers large-scale collisions, a power-law index of around κ = 2.25, and with some reasonable assumptions about measure, appreciable weight for N coll > 1 and observationally interesting values of R ls .
The methods developed in this paper can be used to study the detailed phenomenology of bubble collisions in many models of eternal inflation. The model studied in this paper is simple, but it is finely tuned and features a somewhat unnatural hierarchy of scales between the inflationary and barrier regions of the potential. There is an extremely large landscape of more complicated, and perhaps more realistic, inflationary models, featuring for example multiple scalar fields, branes, fluxes, and extra dimensions. It should be straightforward to extend our program to cover some of these cases. In particular, the code can simulate multiple scalar fields with virtually no modification. More studies will be necessary to determine whether the results presented here for the template and prior are universal across different models.
The results of this paper will be instrumental in obtaining the best limits possible from searches for bubble collisions in data from CMB experiments and other cosmological probes. A previous Bayesian analysis [10,11,40,41] of WMAP 7-year data relied on assumptions about the prior. The results of this paper support some of the assumptions made about the form of the template and prior, but invalidate others. Future work on the universality of our findings will lead to a clear picture for the template and prior, allowing for the most stringent possible tests based on data from the Planck satellite.
The enormity of the implications for detecting a bubble collision signature as a direct experimental test of eternal inflation is obvious; even just the formulation of an observational test of the resulting Multiverse is a breakthrough. Feedback between theory, simulations, and observational tests will allow us to determine how the fundamental theory is constrained, or if a detection is made, to understand what has been discovered.