Cosmic bubble and domain wall instabilities I: parametric amplification of linear fluctuations

This is the first paper in a series where we study collisions of nucleated bubbles taking into account the effects of small initial (quantum) fluctuations in a fully 3+1-dimensional setting. In this paper, we consider the evolution of linear fluctuations around highly symmetric though inhomogeneous backgrounds. We demonstrate that a large degree of asymmetry develops over time from tiny fluctuations superposed upon planar and SO(2,1) symmetric backgrounds. These fluctuations arise from zero-point vacuum oscillations, so excluding them by enforcing a spatial symmetry is inconsistent in a quantum treatment. We consider the limit of two colliding planar walls, with fluctuation mode functions characterized by the wavenumber transverse to the collision direction and a longitudinal shape along the collision direction $x$, which we solve for. Initially, the fluctuations obey a linear wave equation with a time- and space-dependent mass $m_{eff}(x,t)$. When the walls collide multiple times, $m_{eff}$ oscillates in time. We use Floquet theory to study the fluctuations and generalize techniques familiar from preheating to the case with many coupled degrees of freedom. This inhomogeneous case has bands of unstable transverse wavenumbers $k_\perp$ with exponentially growing mode functions. From the detailed spatial structure of the mode functions in $x$, we identify both broad and narrow parametric resonance generalizations of the homogeneous $m_{eff}(t)$ case of preheating. The unstable $k_\perp$ modes are longitudinally localized, yet can be described as quasiparticles in the Bogoliubov sense. We define an effective occupation number to show they are created in bursts for the case of well-defined collisions in the background. The transverse-longitudinal coupling accompanying nonlinearity radically breaks this localized particle description, with nonseparable 3D modes arising.


Introduction
We study the behaviour of linearized fluctuations around colliding scalar field domain walls in situations where the walls possess a high degree of spatial symmetry. Domain walls arise when discrete symmetries are spontaneously broken and are familiar from the magnetic domains formed in ferromagnets. In the context of the early universe they can form in both high-temperature and vacuum phase transitions, either through self-ordering dynamics following a rapid quench or as the walls of nucleated bubbles during a first-order transition. However, if the walls are stable, their energy density falls only as the square of the expansion factor, strongly constraining the allowed outcome of these phase transitions. Domain walls and similar objects such as D-branes are also a common ingredient in early universe model building. Examples include braneworld cosmologies in which our observable dimensions are confined to either a lower-dimensional brane or a domain wall embedded in a higher dimensional space [1][2][3][4][5], inflationary cosmologies including stacks of D-branes [6,7], and some cyclic cosmologies [8].
In a complete theory the domain walls interact with their own dynamics, either inherited from an underlying scalar field theory or intrinsically in the case of D-branes. When several such walls are present, the evolution may result in collisions. In some cases, such as the selfordering after a quench or a rapid percolating first-order phase transition, the domain walls form a complicated network with interactions and collisions occurring in a wide variety of orientations. However, in other scenarios the collisions possess a large amount of symmetry, such as planar symmetry or SO(2,1) symmetry. These highly symmetric configurations may arise from tuning of the initial conditions as in braneworld cosmologies. In other cases, the dynamics naturally lead to symmetric collisions, although the underlying theory might still require tuning to realize the appropriate limit (e.g., of planarity). An example of this is tuning to ensure a slow first order phase transition where collisions are infrequent and the bubbles expand to several times their initial radius before colliding.
Here our focus is on the particular case of colliding parallel planar walls formed by the condensate of some scalar field φ. The qualitative behaviour of the fluctuations around the planar walls also carries over to the case of collisions with an SO(2,1) symmetry. These two symmetry assumptions -planar and SO(2,1) -are widely invoked to study collisions in braneworld scenarios [5,[9][10][11] and false vacuum decay [12], respectively. In both cases, assuming so much symmetry reduces the underlying field equations to a one-dimensional nonlinear wave equation. This reduction greatly simplifies the problem and has been central to many past studies of domain wall collisions. We use the solutions to these one-dimensional nonlinear wave equations as backgrounds upon which our fluctuations are superposed.
An important difference between the classical and quantum problems, which has been neglected in previous work, is the extent to which the dynamical evolution preserves the initially assumed symmetries. While the classical dynamics may possess exact planar or SO(2,1) symmetry, the quantum fluctuations break the symmetry in any realization, albeit at quite a small level initially. These fluctuations may be subject to unstable growth as the combined background plus fluctuation system evolves, with the early stages of the instability described by the linearized perturbation equations we solve for in this paper. If the fluctuations grow they will have considerable backreaction and rescattering effects on the symmetric part of the field, which can significantly modify the overall dynamics. A proper treatment of these effects requires studying the nonlinear problem and is the subject of two companion papers [13] and [14].
We restrict our considerations to two different single-field scalar theories possessing domain wall solutions. We refer to this field φ as the symmetry breaking scalar field. The background spacetime is assumed to be Minkowski throughout. In addition to the choice of underlying theory, the evolution of the fluctuations depends on the particular background around which we expand the fluctuations. Therefore, we consider a variety of collisions in each potential. We show that nonplanar fluctuations in φ can experience exponential instabilities for a broad class of collisions.
Our analysis uses Floquet theory applied to a non-separable PDE. This generalizes the techniques used in preheating, where the spatial homogeneity of the background results in Floquet theory applied to the ODE for the fluctuations. We find generalizations of broad parametric resonance and narrow parametric resonance to the case of fluctuations around a spatially inhomogeneous background.
Although we focus on two specific scalar field models, the dynamical mechanism that leads to the rapid growth of fluctuations is much more general. As we will explicitly demonstrate, the broad parametric resonance instability is essentially particle production in the Bogoliubov sense for fluctuations bound to the walls. These fluctuations are the transverse generalization of the Goldstone mode arising from the spontaneous breaking of translation invariance by the domain wall. Therefore, these modes exist for any membrane-like structure appearing in a translation invariant theory. Examples of such membrane-like structures include domain walls in other field theories and D-branes in string theory. When the two "branes" are well separated, the fluctuations are trapped by an effective potential well. As long as the shape of these wells is modified by the collision we expect similar instabilities to arise regardless of the underlying theory.
The remainder of this paper explores the rich dynamics of linear fluctuations around colliding domain walls. We first introduce our two models in section 2 and present the domain wall solutions that each potential supports. In section 3 we introduce our decomposition of the field into a background and fluctuations, followed by a review of the background dynamics. The central analysis in contained in section 4, where we use Floquet theory to understand the dynamics of the fluctuations. We provide instability charts for the fluctuations and study the mode functions in detail. We also comment on the applicability of our results to a broader class of theories. Finally, in section 5 we briefly comment on the implications for SO(2,1) bubble collisions. We present our conclusions of this linear study in section 6. Some of the more technical details explaining the construction of approximate background solutions are contained in appendix A. Details of our numerical methods and convergence tests demonstrating their superb accuracy and convergence properties are in appendix B.

Model Lagrangians and Domain Wall Solutions
In this section, we introduce the two potentials we consider, review the domain wall solutions each supports and discuss the types of perturbations that exist around these solutions. Since we will ultimately work in three spatial dimensions, we also discuss the embedding of lowerdimensional domain wall solutions in three dimensions.
Our first choice of potential is the sine − Gordon model : which supports a family of static inhomogeneous solutions -kinks -with profiles given by φ SG kink (x) = 4φ 0 tan −1 (e √ Λφ −1 0 (x−x 0 ) ) + 2πn, n ∈ Z . (2.2) These solutions interpolate between neighbouring minima of the potential (2.1) with φ(∞) = φ(−∞) + 2π, and they are the one dimensional version of domain walls. Here x 0 determines the spatial position of the kink and n is an integer determining which minima the kink interpolates between. There is also a corresponding antikink solution which is obtained by the substitution (x − x 0 ) → −(x − x 0 ). Kinks moving at a constant velocity can be obtained by Lorentz boosting the static solution.
At linear order in one spatial dimension the only normalizable localized perturbation of the kink is the zero mode corresponding to an infinitesimal translation of the centre of mass We will later consider planar kink solutions in two or more spatial dimensions. For planar walls in higher dimensions, these localized perturbations give rise to a spectrum of bound  Figure 1. Plots of the double-well potential for several choices of the parameter δ controlling the difference in potential energies between the two wells. Domain wall solutions interpolate between spatial regions where the field is near the false vacuum at φ f alse ≈ −1 and regions where it is near the true vacuum φ true ≈ 1.
state fluctuations with dispersion relationship ω = |k ⊥ |, where k ⊥ is the 2D wavenumber parallel to the wall. Our second example of a potential supporting domain wall solutions is figure 1. δ is an adjustable parameter controlling the difference between the false and true vacuum energies and V 0 is a constant. 1 As long as δ is not too large, this potential supports spatially dependent field configurations in which the field is localized near each of the minima in different regions of space, with the requisite domain-wall structures interpolating between the different regions. A well known example occurs for δ = 0 in one spatial dimension. In this case, the kink solution located at x 0 is given by with a corresponding antikink solution again obtained by replacing ( Moving kinks are again obtained by Lorentz boosting. Unlike the sine-Gordon model, which posesses a single bound state excitation, the one-dimensional double-well kink has two localized normalizable linear perturbations. There is also a continuum of delocalized radiative modes with frequencies ω 2 ≥ 2. The localized perturbations are often referred to as the translation mode 1 Unless explicitly indicated, for the remainder of the paper we measure all dimensionful quantities in units of mnorm, with the exception of the fields measured in units of φ0 and the potential in units of m 2 norm φ 2 0 . The natural mass scale mnorm is given by mSG = √ Λφ −1 0 for the sine-Gordon potential and m = √ λφ0 for the double-well potential. and the shape mode The translation mode corresponds to a spatial translation of the centre of the kink and is thus analogous to the sine-Gordon zero mode. The shape mode is an internal excitation which can be thought of as an oscillating wall width. Kink-like solutions continue to exist as we deform the potential by increasing δ. However, the potential energy difference between the two minima causes the "wall" to accelerate, so the kink solutions are no longer time-independent in inertial reference frames. In three spatial dimensions domain walls become embedded two-dimensional hypersurfaces with some small but finite width. We consider two cases that possess a high degree of spatial symmetry. Our main focus is planar walls generated by extending the sine-Gordon and double-well kink solutions to the additional two spatial directions. We will also study bubbles of "true vacuum" nucleating within the false vacuum in the double-well potential, restricting to choices of δ in the double-well potential for which these false vacuum bubbles are well described by the Coleman-deLuccia (CdL) instanton [15][16][17]. In Minkowski space at zero temperature, the most likely initial bubble profile possesses an SO(4) symmetry in Euclidean signature [18] with profile determined by where r 2 E = r 2 + τ 2 , τ is the Euclidean time and d is the number of spatial dimensions. The initial bubble profile is obtained by analytically continuing back to real time. In the thin wall limit (valid if the initial radius of the bubble is much greater than the thickness of the wall), the friction term in (2.8) is dropped and we are left with the same equation as for a domain wall in the corresponding 1 + 1-dimensional theory. In this limit, the initial bubble radius is given by where ∆ρ is the difference in energy density between the true and false vacuum and σ is the surface tension of the wall. In the final equalities on each line we made use of the specific form of the potential (2.4).

Dynamics of Planar Symmetric Collisions
We now study (nonplanar) fluctuations around colliding parallel planar domain walls. We treat this case first because the scale associated with the overall radius of the bubbles does not enter the problem, so parallel domain walls constitute a slightly simpler arena in which to illustrate the underlying fluctuation dynamics. In the limit that the bubbles have radii much larger than any other relevant scale in the problem we also expect bubble collisions to be reasonably approximated by two colliding planar walls.

General Formalism
Our setup consists of a kink starting at x = −x init and moving to the right and an antikink starting at x = x init and moving to the left. For ease of nomenclature, we refer to this as a KK pair. We take the collision axis to be the x direction and split the field as where the fluctuations satisfy δφ(x, y, z, t) = 0. The planar symmetry allows us to approximate ensemble averages with averages over the (y, z) plane. Before solving for the fluctuation field we must find the background solution around which to perturb. If backreaction and rescattering effects are ignored the background field φ bg undergoes the same dynamic evolution as in 1+1-dimensions, namely with initial conditions given by the KK pair. Meanwhile, the linearized equation for the fluctuations is where δφ k ⊥ (x, t) = dydze −i(kyy+kzz) δφ is the 2D Fourier transform of δφ in the directions transverse to the collision axis and k 2 ⊥ ≡ k 2 y +k 2 z . We refer to these planar symmetry breaking fluctuations as transverse. At this level of approximation, the fluctuation δφ behaves as a free field with a time and x-dependent effective mass (V (φ bg (x, t)) determined independently by the background evolution.
If the system is discretized in the x direction, it is just a collection of coupled oscillators. In real space this coupling occurs via our choice of discretization of the Laplacian term, while in momentum space (along x) the oscillators couple via the Fourier transform of V (x, t). This makes it evident that, for a given V (x, t), the (time-dependent) 'normal' mode oscillations fully describe the system. Of course, the rate at which the discretized system converges to the continuum result depends upon our choice of spatial discretization. Nonetheless, we will show that this approach can be carried out approximately and it provides useful insights into the behaviour of the fluctuations.
In this paper we use a Fourier pseudospectral approximation to discretize ∂ xx , thus obtaining exponential convergence as we increase the number of grid sites for a fixed box size. Details about our precise numerical procedures as a well as a demonstration of the superb convergence properties of our techniques can be found in appendix B.

Dynamics of the Planar Background
We now briefly review the background dynamics for colliding planar domain walls in our two chosen potentials. These background solutions determine the form of V (φ bg (x, t)) to use as input in the fluctuation equation (3.3). The planar symmetry allows us to reduce the background dynamics to that of a KK pair interacting in 1+1-dimensions.

sine-Gordon Model
This model is particularly useful because in one spatial dimension it is integrable and the kink and antikink solutions are true solitons -they interact with each other while preserving their shapes at early and late times. More importantly, there exist analytically known periodic breather solutions where v > 0 is a free parameter determining the properties of the breather solution and [19]. For v 1 the kink and antikink pair are well separated at t = 0 and have initial positions x init ≈ ± √ 1 + v 2 ln(0.5v), whereas for larger v they are much more tightly bound and the breather is a localized oscillating blob of field with size The edge of the breather is defined to be the point where φ breather (t = 0, r breather ) = φ edge .
Some representative examples are shown in figure 2.

Double Well Model
The background evolution in this potential (2.4) is more complicated, but it has been studied by many other authors who detail a rich phenomenology [20][21][22][23][24]. The kinks in this model are solitary waves rather than true solitons, and thus they emit radiation when they interact. Combined with excitation of the internal mode during collisions, this means that the motion is no longer exactly periodic. Because we do not have exact solutions to (3.2), we must resort to numerical simulations. We use a Gauss Legendre time-integration combined with a Fourier pseudospectral discretization, which allows us to obtain machine precision results for both the spatial discretization and time-evolution. Details are provided in appendix B. We restrict our attention to a kink-antikink pair characterized by an initial separation and specified initial velocities. For more general setups, the interested reader should consult the previously cited works for the wide range of interesting phenomenology that can arise.
To simplify our analysis we work in the centre of mass frame. We take the initial kink and antikink speeds u and separations d sep as free parameters. Some illustrative examples of the dynamics are shown in figure 3. In the symmetric well the attractive force between the kink and antikink decreases exponentially at large separations, and as a result unbound motions with the kink escaping back to infinity are possible. For low initial velocities (u 0.15), the K andK always capture each other after colliding and do not escape back to infinity -see figure 3, top left, for an example with u = 0.05. Rather than immediately annihilating into radiation, the kinks bounce off each other several times and then settle into a long-lived oscillatory blob known as an oscillon (here living in one dimension). During each oscillation, some energy is radiated away, so this localized state eventually decays. However, the rate of energy loss is slow and the oscillons can persist for thousands of oscillations (or more) before finally disappearing.
As the incident velocity increases, the kinks enter a "resonant escape" regime characterized by bands of incident velocities in which the two kinks eventually escape back to infinity. The escape bands are separated by bands in which the KK pair trap each other and form an oscillon as in the low velocity limit (see e.g. [23]). Within each escape band, the number of bounces the walls undergo before escaping back to infinity is a very complicated function of the incident velocity. This seemingly strange behaviour is usually attributed to the shape mode. At each collision, nonlinear interactions transfer some of the kinks' translational kinetic energy into a homogeneous excitation of the shape mode or vice versa. The direction of energy transfer depends on the oscillation phases of the two shape modes (one on the kink and one on the antikink). As a result, the kinetic energy of the kinks can decrease in one bounce as the shape mode is excited causing the KK pair to become temporarily trapped. In the subsequent collision some of this stored energy can then be transferred back into overall translational motion, giving the kink and antikink enough translational kinetic energy to escape back to infinity. An example of this behaviour for u = 0.2 is illustrated in the upper right panel of figure 3. In the first collision, the KK pair loses some energy to radiation and excites internal shape modes. At the second collision, energy is transferred from the shape mode back into translational kinetic energy and the kinks escape each other.
A critical velocity exists above which the KK pair always interacts exactly once before escaping back to infinity. During the collision, some of the energy escapes as radiation. The remaining energy flows between the shape and translational modes of the kink. Provided that the shape modes are not initially excited, the outgoing velocities of the two walls are thus always less than their incident velocities. This sort of interaction appears in the bottom left panel of figure 3.
When we take δ > 0, (i.e. make the double well asymmetric), the difference in vacuum energies across the kink causes an acceleration toward the false vacuum side. The kink and antikink experience an approximately constant attractive force at large separations. As a result, the kinks are no longer able to escape back to infinity and will always undergo multiple collisions while slowly radiating energy. Eventually, an oscillon forms at the location of the original collision. See the bottom right panel of figure 3 for an example of kinks interacting in an asymmetric double well.
In all cases discussed above the shape mode is excited by the collisions. This is visible as an oscillating wall thickness, as seen in figure 3.

Dynamics of Linear Fluctuations
In the previous section we briefly reviewed the dynamics of planar symmetric background field solutions. The key feature for our analysis is the presence of oscillatory behaviour: the walls typically bounce off each other many times or settle into a localized bound state rather than immediately annihilating, and internal vibration modes of the walls are excited in the collisions. We now study the transverse fluctuations described by (3.3) in these types of oscillating backgrounds Our focus is whether small initial (transverse) perturbations to the background dynamics illustrated in figure 2 and figure 3 can be amplified to the extent that they become important to the full field dynamics. The results presented in this section are the main findings of our analysis. We plot (in)stability charts for linear fluctuations around various choices of planar collision backgrounds for both the sine-Gordon and doublewell potentials. In addition, we study several representative choices of the unstable mode functions, which provides insight into the qualitative structure of the instability charts. We find inhomogeneous generalizations of both broad and narrow parametric resonance. Our results demonstrate that the transverse fluctuations in φ can grow rapidly as a result of collisions.
Amplification can happen in several ways. Radiative modes can be excited in the collision region and then subsequently propagate into the bulk, which is the transverse generalization of the outgoing radiation seen in the (1+1)-d simulations above. However, these excitations carry energy away from the collision and cannot experience sustained growth. More interesting is the pumping of fluctuations bound to the kinks (in cases where the walls separate widely) or bound in the oscillating effective mass well created by the late-time onedimensional oscillons. A similar effect for fermionic degrees of freedom has been studied by several authors [25][26][27], and for an additional field coupled to φ in [28]. We only consider fluctuations in the field φ itself, which is required in a consistent quantum treatment, and we do not appeal to additional fields coupled to φ in order to obtain growing fluctuations. As well, unlike the authors of [28] we include the deformation of the spatial structure in V and couplings between the bound modes and radiative modes.
Although strictly true only for the case of the breathers, we approximate the oscillatory evolution as periodic. Because of the periodicity, we can use Floquet theory to quantitatively study instabilities in δφ. The resulting Floquet modes then provide a time-dependent normal mode decomposition in which the evolution of the system is simple. The periodic approximation is very accurate for the late-time bound states and excitations of the shape mode in the planar background, so Floquet analysis is well justified for those two cases. When the walls repeatedly bounce off of each other, we show that the amplification happens in a very short time interval around the time of collision. The short bursts of growth in the fluctuations remain when the bouncing is not periodic, so again the Floquet analysis provides an approximate account of the full evolution.
In early universe cosmology and nonequilibrium field theory, Floquet analysis is familiar from the theory of preheating [29][30][31][32]. For these problems the background is spatially homogeneous, which leads to two simplifications. First, the background evolution is described by a nonlinear ODE instead of a nonlinear PDE, and analytic solutions are known in many interesting cases [33][34][35][36]. Second, the equation for the fluctuations is diagonalized by the 3D Fourier transform. Therefore, instead of involving a large number of coupled oscillators the problem reduces to a collection of decoupled oscillators with periodically changing masses ) depends only on time and k 2 = k 2 x + k 2 y + k 2 z is now the square of the full three-dimensional wavenumber. Solutions to this equation have the well known form δφ k (t) ∼ P (t)e µt with the (possibly complex) exponent µ known as a Floquet exponent. P (t) satisfies P (t + T ) = P (t), where T is the period of the oscillation in m 2 ef f . The µ's depend parametrically on k as well as the functional form of m 2 ef f . Typically bands of stable (max[Re(µ k )] = 0) and unstable (max[Re(µ k )] > 0) modes appear as we vary k while holding the form of m 2 ef f fixed. 2 There are a variety of underlying mechanisms that can be responsible for the instability, such as tachyonic resonance, narrow resonance and broad parametric resonance [30,31,35].
For a spatially dependent effective mass, which is the case we consider, the form of the decoupled solutions generalizes to δφ = F (x, t)e µt , where the profile function F (x, t) satisfies F (x, t + T ) = F (x, t). This is easily seen by discretizing the x direction and placing the system in a finite box. The equations of motion then take the formv = M (t)v, with M (t) a periodic matrix and v = δ φ T , δ˙ φ T T . Such equations fall within the domain of Floquet theory and therefore solutions of the form given above are known to exist (modulo issues of convergence on taking the continuum limit).
To treat the case of spatially inhomogeneous masses, we first discretize the fields δφ k ⊥ on a lattice of N points labelled by index i. We denote the field value at the ith lattice site by φ i and the corresponding field by the vector δ φ. Next, we consider 2N linearly (3.3). Using this set of solutions, we form a 2N x2N fundamental matrix solution F (t). For simplicity, we choose our initial conditions such that the jth row of F (t) is given by the solution with initial condition Of course, we can construct our fundamental matrix from any complete set of initial states, and we verified that choosing an orthonormal basis in Fourier space reproduced the results we present below. Finally, the Floquet exponents are related to the eigenvalues Λ n of F (0) −1 F (T ) via Λ n = e µnT , with the initial conditions for the mode functions given by For each choice of k 2 ⊥ and effective mass there are 2N such exponents, but we focus on the exponent with the largest real part (i.e. the largest Lyapunov exponent). For notational simplicity we refer to this maximal real part of a Floquet exponent as µ max ≡ max [Re(µ)]. This method allows us to find any exponentially growing instabilities given some fixed background evolution, but it is completely blind to other more slowly growing instabilities. In particular, power law growth results when there are degenerate Floquet exponents, and we might expect this to be common in the continuum limit given that we then have infinitely many oscillators. However we expect the exponentially growing modes to be the most important dynamically and thus most interesting for our purposes. After all, we are ultimately interested in the quantum problem where we must integrate over all possible modes.

sine-Gordon Potential
Using the planar symmetric sine-Gordon breather as the background circumvents the problem of finding appropriate approximations to φ bg . Furthermore, as bound states of a kink-antikink pair, breathers resemble the states arising from collisions in the double-well and can be used to gain some qualitative insight into the dynamics of the fluctuations for the double-well as well. The exact equation for linear fluctuations around the breather is lows from the fact that the linear operator M (t) = 0 I ∂xx − V 0 describing the evolution of (δφ, δφ) T has Tr(M ) = 0. We must then have i µi = 0 so that the presence of a single negative Lyapunov exponent necessarily creates a positive Lyapunov exponent. For the homogeneous case this requirement becomes max[Re(µ k )] = 0.
The initial profile of V (x, t) for several representative choices of v can be found in the right panel of figure 4, with the subsequent time evolution for these same three choices in figure 5, figure 8 and figure 10.
The left panel of figure 4 is an instability chart in which µ max T breather appears as a function of k 2 There are several generic features of note. First, µ max = 0 when k 2 ⊥ = 0 for all values of the parameter v considered. Since the k ⊥ = 0 modes are part of the planar symmetric system, this is a validation of our numerical procedure since no exponentially growing modes exist in the sine-Gordon model in 1+1-dimension. 3 For v −1 ≤ 1, V (x, t) has the form of a single well whose depth oscillates with time. There is a single instability band, the width and strength of which increases monotonically as v decreases. Once v −1 > 1 the kink and antikink become less tightly bound and V develops a pair of local minima separated by a small bump. An additional instability band then appears on the chart. As we continue to increase v −1 (decrease v), the kink and antikink begin to separate from each other during the motion, and the two local minima in V develop into two well separated wells as seen for v = 0.01 in the right panel of figure 4. As we continue to decrease v, additional instability bands appear. Each of these bands grows quickly with increasing v −1 and approaches a nearly constant width, although the strength of the instability (per breather period) within each band continues to increase. However, T breather is simultaneously increasing ∼ v −1 √ 1 + v 2 and this results in a decrease of the maximal Floquet exponent if measured in units of t.
The left portion of the chart (v −1 1) corresponds to the small amplitude breathers while the right portion (v −1 1) corresponds to the large amplitude breathers. Right: V (φ breather (x, 0)) for three representative choices of v. For v ≥ 1, there is a localized blob with a single minimum. We have plotted v = ( √ 2 − 1) −1 and v = 1 corresponding to the cases when the middle of the breather just reaches V (φ) = 0 and V (φ) = −1 respectively. For smaller v, the single blob instead develops a pair of minima, with the formation of two distinct wells when v 1 corresponding to the well-separated kink and antikink.
We now study the properties of the mode functions for various regimes of v. This gives us insight into the chart's qualitative features and sheds light on the dynamical mechanism responsible for the amplification.

v 1 : Broad Parametric Resonance
First we consider the v 1 limit. The breather is effectively a weakly bound kink-antikink pair which repeatedly approach and pass through each other. The resulting evolution of V is illustrated in the left panel of figure 5 for the specific case of v = 0.01. Due to the reflection symmetry of the potential about any of its minima, the period For most of the evolution, V has two distinct wells corresponding to the individual kink and antikink. In the 1d field theory, each of these wells has a zero mode (the translation mode) associated with it. When we allow for the transverse fluctuations, the zero mode leads to a continuum of bound excited states with dispersion relation ω bound = k ⊥ > 0. As the breather evolves, the two wells in V periodically come together and "annihilate" each other as seen in the top left panel of figure 5. During these brief moments when the two wells have disappeared, the bound states cease to be approximate eigenstates of the system. The temporary annihilation of the wells allows for the rapid amplification of bound fluctuations. Whether or not a particular k ⊥ receives coherent contributions to its amplitude at successive annihilations depends on the phase ω bound T V ∼ k ⊥ √ 1 + v 2 /v accumulated between collisions. In the stability chart this dependence on the accumulated phase manifests as the repeating structure in k 2 Of course, this process is very similar to the familiar case of broad parametric resonance [31], in which short but large violations of adiabaticity (|Ω/Ω 2 | 1) of a harmonic oscillator with a time-dependent frequency Ω(t) lead to bursts of bound state "particle production" in the Bogoliubov sense. For the familiar homogeneous case the nonadiabaticity is captured by a time-dependent frequency alone, while for the inhomogeneous background we have the additional effect of a strong deformation of the spatial properties of the effective mass.
To illustrate the broad resonance process with a specific example, we plot various aspects of the fastest growing transverse Floquet mode for a breather with v = 0.01 and k 2 ⊥ = 0.05 in figure 5 and figure 6. As evident in the bottom panel of figure 5 the mode function is strongly peaked around the locations of the kink and antikink. The mode function also experiences large jumps around each collision, exactly as expected from our discussion in the previous paragraph. To further illustrate the rapid amplification of bound fluctuations by the collisions, the left panel of figure 6 shows the value of δφ at the leftmost instaneous minimum of V as a function of time, along with the value V min at this minimum. During the collision, V min rapidly changes from −1 to 1 and δφ experiences a nearly instantaneous increase in its amplitude.
To quantify the growth of fluctuations, it is useful to introduce the effective particle number n ω bound ef f which (modulo overall normalization and contributions from the bulk) gives the occupation number of massless transverse fluctuations bound to either the kink or antikink. For a bound fluctuation on an isolated kink with transverse wavenumber k ⊥ this quantity is constant. Therefore, we can use it to track the growth of fluctuations so changes. From the right panel of figure 6, we see that between collisions n ω bound ef f is constant and undergoes nearly discontinuous jumps during the short collisions between the kink and antikink.
This analysis confirms our intuitive expectation that the exponential growth results from the production of bound state fluctuations. However, several additional questions remain γ , is shown as a red line. Bottom: The corresponding mode function, illustrating the constant amplitude oscillations around the kink and antikink and the rapid growth during the short interval when they collide. Radiation moving away from the collision is also visible. Top right: Initial conditions for the mode function, illustrating both the localization near the kink and antikink, and the extended radiating tail.
that we now address. First, since the mode function is symmetric about the origin it cannot distinguish between the following processes: simultaneous amplification of the bound states on both the kink and antikink, amplification of modes bound to the kink (or antikink) only, and growth of fluctuations in only the left (or right) well. 4 Second, it is unclear whether the stability bands arise primarily from dissipation of fluctuations into the bulk, or primarily from phase interference between bound fluctuations. Finally, a somewhat surprising feature of the the mode function is the long radiative tail that extends far from the spatially localized breather into the bulk. To gain insight into these issues, we consider fluctuations with the initial condition . This initial condition corresponds to a bound state fluctuation on the leftmost kink (or antikink) of the breather solution. Figure 7 shows the evolution of this initial state for several choices of k 2 ⊥ . For each run, we show the evolution of n ω bound  5 Finally, we show the profile of δφ 2 + δφ 2 at several times around the first collision of the kink and antikink. This is a useful indication of the growth of fluctuations and in light of (4.3) can be interpreted as a local "particle density".
In each collision some radiation is released from the collision region. As k ⊥ is increased the fluctuations appear to be less tightly bound to the kink, and the amplitude of radiation produced in the collisions tends to be larger. This is likely the origin of the decreasing amplitude of µ max at the centre of the instability bands as we increase k ⊥ while holding v constant. Incidentally, the emitted radiation clarifies the origin of the long radiative tail we found for the mode function in figure 5. Whether or not the mode functions grow, at each collision some radiation is emitted from the collision region. Consider a pair of subsequent collisions for the case when the mode function is growing exponentially. Just prior to both collisions, the phase of the mode function will be the either the same (or differ by π) while the amplitude differs by e µT V . These properties follow from the form P (x, t)e µt of the mode function and the step-like nature of the increase in the amplitude. For simplicity, consider the case when the phases are the same. 6 Since the fluctuations obey a linear equation, the radiation emitted in the second collision is identical to that in the first collision, except it has a larger amplitude by a factor e µT V . Therefore, trains of radiation with the same spatial profile but step-like behaviour in their amplitude will be emitted from the collision. Figure 7 also suggests that whether or not a particular value of k ⊥ will experience an exponential instability is primarily determined by the phase interference between the bound state fluctuations. This can occur either because the fluctuations are not excited in any of the collisions (as in the case k 2 ⊥ = 4.006x10 −4 ), or because the fluctuations are excited in one collision but then de-excited due to phase interference in the subsequent collision (as in the case k 2 ⊥ = 0.83). Finally, as expected from studying the second largest Floquet exponent (see footnote 4), we see that for the unstable modes the excitation tends to occur on both the kink and antikink simultaneously. On the other hand, when the modes are drawn from a stability band the fluctuations on the kink and antikink no longer experience the same degree of excitation or de-excitation at each collision.

v ≥ 1 : Narrow Parametric Resonance
We now turn to the case of larger v's, where we no longer have a well-defined kink and antikink at any point in the breather's motion. We first consider v = 1 and k 2 ⊥ = 0.35, which is located near the centre of the instability band. For this choice of v, the middle of the breather just makes it to a maximum of the potential every half oscillation. Therefore, V (x, t) has the form of a single oscillating well whose middle oscillates between −1 and 1 and asymptotes to 1 at ±∞, as illustrated in the top left panel of figure 8. The kink and antikink are tightly bound so they never have separate identities. Therefore, the notion of particles bound to the kink and particles bound to the antikink is ill-defined, and our previous intuition based on the creation of fluctuations bound to the individual kink and antikink no longer applies. Instead, we expect the pumping to occur more smoothly and be localized in the region of the breather. As seen in the bottom panel of figure 8, this is indeed the case. The mode function looks like a smoothly oscillating function with exponentially growing amplitude. As well, the mode function satisfies δφ(x, t + T V ) = −δφ(x, t)e µT V , so the period of the oscillation is the period of the breather not the period of the effective mass. One way to see the smooth exponential growth is to consider the quantity n (ω breather ) ef f is an effective particle number, like (4.3) for the v 1 breathers, but modified to account for the fact that in this case the oscillation frequency of the breather dominates the oscillation frequency of the mode function. The top right panel of figure 8 demonstrates the smooth exponential growth of n (ω breather ) ef f with some small subleading oscillations. The periodic part of the Floquet modes can be decomposed into an harmonic expansion dominated by the frequency of the breather, with smaller subleading contributions from higher harmonics,  . Evolution of δφ for initial condition δφ init = sech(x + x K ) with x K = − √ 1 + v 2 log(v/2). We used v = 0.01 and four values of k 2 ⊥ to illustrate various types of behaviour. For each k ⊥ we plot: the effective particle number defined in (4.3) as well as the fraction of "particles" in the left and right half of the simulation (left), the value of the field at the instantaneous locations of the two minima in V (middle), and the local "particle density" (δφ 2 + k 2 ⊥ δφ 2 ) for several times around the first collision (right). In the first row k 2 ⊥ = 4.006x10 −4 , which is located in the first stability band; in the second row k 2 ⊥ = 6x10 −4 , which is near the maximum of the second instability band; in the third row k 2 ⊥ = 0.05 which is in one of the higher order instability bands; in the fourth row k 2 ⊥ = 0.83, which is located in one of the higher order stability bands. T b = T breather in the legend for the right panels.
as shown in figure 9 . This dominant frequency justifies our introduction of the quantity n concentrated at the location of the well created by the breather. Once we move away from the breather, the spatial phase in each frequency varies linearly. This linear variation is consistent with the production of radiation in the breather which then travels off to infinity.
Finally, we consider the case v = 1 √ 2−1 and k 2 ⊥ = 0.2. As in the case v = 1, there is a single oscillating well in V . However, the maximum excursion of the field is ± π 2 , so the effective mass for the fluctuations is positive semidefinite. This is illustrated in the top left panel of figure 10. The remaining panels of figure 10 show that the resulting fastest growing Floquet mode is qualitatively similar to the v = 1 case. There is an isolated blob that oscillates with a frequency determined by ω breather whose amplitude grows smoothly as an exponential. Looking at the frequency content in figure 11, we see that the mode function again consists of a large amplitude part concentrated in the potential well of the breather and a smaller amplitude radiative tail propagating away from the breather. The oscillation frequencies in the two regions are more monochromatic than in the v = 1 cases. The dominant frequency of the radiative part differs from the dominant frequency of the bound part. In both the v = 1 and v = 1 to the case of narrow resonance for a single oscillator.

Symmetric Double Well
We now consider the symmetric double-well. As demonstrated in figure 3, for certain choices of the initial kink speed the background solutions can undergo oscillatory motion. For the case of the asymmetric well, the oscillatory motion is generic due to the attractive force experienced by a well separated KK pair. We will not consider the asymmetric well explicitly, but the behaviour of the fluctuations is qualitatively the same as for the symmetric well. For the double-well, oscillatory motion of the background solution comes in three forms: repeated collisions of the walls, oscillations of the internal planar shape mode, and evolution of a late-time planar oscillon. The first of these is analogous to the sine-Gordon breathers with v 1, while the last is qualitatively similar to the breathers with v ≥ 1. For the repeated collisions and late-time oscillons, the only new feature is that we have V (x = 0) > V (x = ∞) for part of the evolution. The oscillations of the planar shape mode are a new feature not present in the sine-Gordon model. We will study each of these background motions independently, although there are several caveats to this approach. First of all, the repeated collisions and shape mode oscillations typically occur together in actual background solutions. In the homogeneous case, the stability chart for a harmonic oscillator with effective mass containing multiple frequencies is not the same as a superposition of the stability charts for each individual frequency [37]. Despite this, from the sine-Gordon analysis above we know that the growth of fluctuations due to the collision will occur in a very short time interval during the actual collision. The collision also excites the (planar) shape mode in the background, which is then free to pump excitations during the much longer intervals while the walls separate from each other. Since the two amplification effects are temporally separated, we can gain good qualitative understanding by considering the processes in isolation. The resulting interference from the two mechanisms could then be done using projections on the appropriate Floquet basis (either for the wall collisions or the shape oscillations) just before and after each collision. Of course, since the exact background is not exactly periodic, the interference effects will tend to get smeared out leading to a smoothing of the instability diagram. This smoothing is analogous to what happens in the homogeneous case when expansion of the universe is included.

Bouncing of the Walls : Broad Parametric Resonance
Keeping the above caveats in mind, we begin with the case where the KK pair separate widely from each other between collisions. Unlike the sine-Gordon breather, we have no analytic solutions for this case. Also, due to the emission of radiation and the excitation of the shape mode, the background motion is no longer periodic. In order to create a periodic  Figure 11. approximation that is amenable to our Floquet analysis, we take the background to have the following form where r(t) is a dynamical variable and γ = (1 −ṙ 2 ) −1/2 . This ansatz ignores the production of radiation, excitation of the planar shape mode and any additional deformation of the kink profile due to interactions. After several further simplifying assumptions, we arrive at the following approximate solution for r(t) where r min is the (minimum) solution to V ef f (r min ) = V ef f (r max ) with V ef f (r) defined in (A.8) of appendix A. Further details of the construction are provided in appendix A. An alternative approach would have been to simply insert a periodic function for r(t) and study the resulting fluctuation behaviour. However, (4.7) is better justified than a completely ad hoc choice for r(t) because it partially incorporates the full field dynamics.  Figure 12 shows the resulting Floquet exponents for several choices of r max . As expected, the bands are distributed evenly in the phase accumulated between subsequent collisions by bound fluctuations, k ⊥ T walls . There is a very strong instability as k ⊥ → 0, which is not unexpected given that our approximation ignores the radiation and planar shape mode that are excited during the collisions. One interesting new feature is the presence of a second growing mode, which we illustrate by plotting the second largest real part of a Floquet exponent in addition to the largest one.

Planar Shape Mode Oscillations : Narrow Parametric Resonance
Next consider the pumping of fluctuations by oscillations of the planar shape mode. This is a new effect that is not present in the sine-Gordon model. Since the planar shape mode is generically excited (or de-excited) during collisions, this source amplification can occur in conjunction with the nonadiabatic production of fluctuations due to the colliding walls described in the previous paragraph. We parameterize the background motion as As in the case of the repeated collisions described above, (4.9) is not an exact solution of the 1D field equations. Nonlinear couplings in the potential will modify the oscillation frequency of the shape mode and also cause it to radiate. Hence the amplitude will gradually decrease in time, leading to a slowly changing oscillation frequency. It is even possible to tune the amplitude so that the subsequent evolution leads to the creation of a kink-antikink pair in addition to the kink that was initially present [24]. However, provided the amplitude is not too large, the time-varying amplitude and frequency can be approximated as an adiabatic tracing of modes on the instability chart. Hence the chart and corresponding mode functions are a good approximation to the fluctuations in the true background.  Figure 13. Instability chart for planar oscillations of the shape mode. Shown is µ max T shape for Plugging (4.9) into the equation for fluctuations gives (4.10) whose Floquet chart is given in figure 13. In the top left panel of figure 14 we plot V (x, t) for A shape = 0.5. As expected from the interpretation of the shape mode as a perturbation to the width of the kink, V looks like a well whose width oscillates in time. There is also some additional oscillating side-lobe structure. The remaining panels in figure 14 show the fastest growing Floquet mode, the effective particle number, and the various frequency components of the solution. These plots are analogous to those for the v ≥ 1 breathers in figure 8, figure 9, figure 10 and figure 11. Because the mode function completes only half an oscillation per period of the shape mode, the frequency used in the definition of n ω ef f is ω = π/T shape . The mode function displays more structure than sine-Gordon breather mode functions, in particular in the higher harmonics. This additional structure is probably due to the additional sidelobe structure of V . However, the solution is still well described by a single oscillation frequency near the core of the kink, with additional harmonics becoming important as we move into the oscillating sidelobes and then into the radiating regime.

1D Oscillon Background : Narrow Parametric Resonance
Consider fluctuations around the late-time (1-dimensional) oscillon state. In order to approximate the background motion, we first expand the solution about the true vacuum minimum as φ bg = φ min +φ(x, t). Up to an additive constant, the potential forφ then takes the form V (φ) =m T shape Figure 14. For an unstable mode around an oscillating planar shape mode excitation, the same series of plots are shown as for the Floquet modes of the v = 1 (figure 8, figure 9) and v = ( √ 2 − 1) −1 (figure 10, figure 11) sine-Gordon breathers. In the definition of n ω ef f , we used the frequency ω = π/T shape . As for the sine-Gordon breathers, this mode function consists of a well-defined core region near the location of the kink as well as a much smaller radiative tail. Due to the additional spatial complexity of V , the mode function displays more spatial and frequency structure than for the breathers.
Expanding the solution asφ(u, w) = ∞ n=1 n φ n (u, w) and solving the resulting equations order by order in , we find φ osc = (P (x) + g( √ 2mx)) cos( √ 2 1 − 2 mt) For the particular form of the double well potential, we have α = 12/φ 2 0 . In the literature, two choices for the function g have been considered. Fodor et al. [38] assumed that no bounded solutions exist and set g = 0, while Amin [39] instead demanded φ 2 (t = 0) = 0 to set g(x) = 3P (x) 2 /φ 0 . We follow Amin when we plot the second order oscillon instability chart. The resulting equation for the transverse fluctuations is In figure 15 we show stability charts for the fastest growing mode around the oscillon background (4.11). We plot the charts for both the leading order and second-order in approximations to the background. 8 The detailed structure of the instability does display some sensitivity to our choice of approximation for the background. For k ⊥ = 0 and 0.2, the higher order approximation removes a weak instability that was present in the leading order approximation, indicating that it is indeed a better approximation to the background at small . However, for larger the higher-order background is actually more unstable than the leading order approximation. This is merely a reminder that our approximation is asymptotic rather than convergent. When we consider k ⊥ = 0, we see that the additional oscillating frequencies in the second-order background lead to several weak instability bands at small . For larger the main instability band extends for a wider range of k ⊥ and has larger Floquet exponents. This is not surprising since the oscillation amplitude is larger in this approximation and we would thus expect it to drive stronger instabilities.

Comments on Collisions of other Membrane-like Objects
Although we focussed on two specific scalar field models, the dynamical mechanism that leads to rapid growth of the fluctuations is much more general. In particular, for well-separated walls the explosive growth of fluctuations relied only on the presence of bound fluctuations around each of the kinks and violation of adiabaticity for these bound states. Recall that these fluctuations are the transverse generalization of the Goldstone mode (i.e. the translation mode) resulting from the spontaneous breaking of translation invariance by the kink in the one-dimensional theory. An equivalent Goldstone mode will occur for kinks in any translation invariant theory, and thus these bound fluctuations are ubiquitous. For kink-antikink type collisions such as those studied here, we typically expect large deformations in the shape of V during each collision. Therefore, the parametric amplification of wall-bound fluctuations we have found here should be very common. We leave for future work the interesting question of kink-kink collisions, where the effective potential wells resulting from each of the kinks do not completely annihilate during collision.
Our results also have relevance to collisions of other membrane-like objects such as Dbranes. When the kink and antikink are well separated, the transverse translation modes are characterized by a spatially dependent location for the centre of the kink, and are thus well-described by an effective action for a membrane. If two such membranes are in close proximity to each other, it is natural to expect them to interact. For the case of D-branes, this interaction is usually described in terms of the excitation and production of string modes. Since string production is a local process from the viewpoint of a field theory on the brane, we expect that the resulting fluctuations will be inhomogeneous and analogous to the excitation of our transverse modes. We show in [13] and [14] that the inhomogeneity of the growing fluctuations dramatically changes the full three-dimensional dynamics compared to the planar approximation. The detailed effects of the onset of nonlinearities amongst the fluctuations is dependent on the details of the high-energy theory (in our case a scalar field theory). However, the linear dynamics described in this paper should be much less sensitive to the UV completion and we expect D-brane collisions to also result in the rapid growth of "transverse" fluctuations that eventually require a full nonlinear higher-dimensional treatment.

Comments on Fluctuations Around Colliding Bubbles
Beginning with the work of Hawking, Moss and Stewart [12], the two-bubble collision problem has been explored by many authors [40][41][42][43][44][45][46][47][48]. A common feature of these analyses is the assumption of SO(2,1) symmetry for the field profiles. This is motivated by the SO(4) symmetry of the minimum action bounce solution for a single bubble [18], which translates into an SO(3,1) symmetry for the nucleation and subsequent expansion of the bubble in real time. The nucleation of a second bubble destroys the boost symmetry along the axis connecting their centres, leaving a residual SO(2,1) symmetry for the 2-bubble solution, making of the problem effectively 1 + 1-dimensional, which greatly eases the numerical challenges and has even led to a general relativistic treatment [49][50][51].
However this is not the full story. Quantum fluctuations are present, in fact are responsible for the bubble nucleation in the first place. We now discuss the linear fluctuation dynamics in the background of the pair of colliding bubbles. For a discussion of perturba-tions around a single bubble see [52][53][54]. Results for the fully nonlinear three-dimensional dynamics of bubble collisions are presented in [14].

Background Dynamics of Highly Symmetric Collisions
As in the planar wall case, we decompose the field into a symmetric background component and a nonsymmetric fluctuation. A convenient variable change is aligned so the centres of the two bubbles both lie on the x-axis. The SO(2,1) symmetry is now manifest, as the background depending only on s and x. We separate the field into a symmetric background and symmetry breaking fluctuations φ = φ bg (s, x) + δφ(s, x, χ, θ).
Ignoring backreaction of the fluctuations, the background solution obeys and the linearized fluctuations evolve according to Here we have factored the perturbation into a sum over eigenfunctions, C ,n labelled by and integer n and eigenvalues : l,n represents an integral over the continuous part of and a sum over n and any discrete normalizable modes C ,n . The curvature of the bubble walls thus influences the fluctuation dynamics in three ways: damping the overall amplitude as s −1 , redshifting the effective transverse wavenumber squared 2 as s −2 , and modifying the dynamics of φ bg and by extension V (φ bg ).
Treatments that assume SO(2,1) symmetry study (5.2) and ignore (5.3). A sample collision between two such bubbles in the asymmetric well with δ = 1 10 is shown in figure 16. As in the case of the kink-antikink collisions, the bubble walls undergo multiple collisions, each time opening up a pocket where the field is localized near the false vacuum minimum. The bouncing behaviour we observe is characteristic of bubble collisions in double-well potentials with mildly broken Z 2 symmetry, and was first noted by Hawking, Moss and Stewart [12].
Considering the implications of this behaviour for the full 3 + 1-dimensional evolution suggests that two instabilities may occur. The first is the generalization of our previous results to the SO(2,1) symmetric rather than the planar symmetric case. Given the background evolution depicted in figure 16, we see that (5.3) again describes a field with an oscillating As in the planar symmetric case, the two walls bounce off of each other multiple times rather than immediately annihilating. During this process, scalar radiation is emitted from the collision region.
time-and space-dependent mass. Figure 16 also reveals the presence of another possible instability. Because of the SO(2,1) symmetry, each pocket with the field near the false vacuum corresponds to a torus with growing radius centred on the initial collision in the full 3-dimensional evolution. Roughly, we can think of this as a torus of false vacuum with a thin membrane separating it from the true vacuum on the outside. The energy difference between the false and true vacuum leads to a pressure acting normal to the local surface of the membrane. Since this pressure wants to push the membrane into the false vacuum, this will tend to cause small initial ripples on the surface of the torus to grow. The surface tension of the membrane and the stretching of the torus as it expands tend to counteract this effect, so that a three-dimensional calculation is required to determine the ultimate fate of these ripples. As we will see in [14], both of these instabilities manifest themselves in the fully 3+1-dimensional problem.

Conclusions
We studied the dynamics of linear asymmetric fluctuations around highly symmetric collisions between planar domain walls and vacuum bubbles. Parallel planar walls are a common ingredient in cosmological braneworld model building, and SO(2,1) bubble collisions are generally believed to be an accurate description of individual bubble collisions during false vacuum decay. The results shown here have implications for such cosmological scenarios, as well as being of intrinsic theoretical interest.
Fluctuations are generically present in the field that forms the domain wall and therefore must be included in a quantum treatment of the problem. However, nearly all past studies of planar wall and vacuum bubble collisions dynamics have used symmetry to reduce the problem to a 1+1-dimensional PDE, thus excluding the fluctuations a priori. Assessing the validity of this drastic reduction in the effective number of dimensions requires a more sophisticated treatment of the problem, and this paper provides such a treatment.
Once we have fixed the symmetric background dynamics for the collision, the fluctuations behave as a free scalar field with a time-and space-dependent effective frequency. Using Floquet theory and extending well-developed methods for ODEs to PDEs, we were able to show that the time-dependence of the effective frequency can lead to exponential growth of the symmetry breaking fluctuations. We also studied the spatial structure of the amplified modes to obtain an understanding of the mechanism responsible for the amplification. We found generalizations of both broad parametric resonance and narrow parametric resonance. Due to the spatial dependence of the effective frequency the amplified modes are localized along the collision direction. For collisions between well defined wall-antiwall pairs, the resulting amplification can be interpreted as Bogoliubov production of particles bound to the walls, i.e., of longitudinally localized modes. The evolution of the effective occupation number shows they are created in bursts when there are well-defined collisions in the background. Once nonlinearity onsets, the strong coupling of the transverse to the longitudinal degrees of freedom breaks this localized particle description. Numerical lattice solutions are needed to describe the subsequent evolution.
Our detailed study of the unstable modes reveals that the dynamical mechanism responsible for the rapid growth of fluctuations is much more general than our two scalar single-field model examples. In particular, for collisions between a pair of well defined walls, the most strongly amplified modes are the transverse generalizations of the Goldstone mode resulting from the spontaneous breaking of translation invariance by the domain wall in the symmetry reduced one-dimensional theory. These modes are present on any membrane-like object in a translation invariant theory. The amplification only relied on a strong deformation of the effective potential binding these fluctuations to the wall. Such deformations will be extremely common in domain walls formed in scalar field theories, and should also arise in collisions between other membrane-like objects such as D-branes. Therefore, we expect qualitatively similar results to be obtained in a wide variety of collisions involving membrane-like objects.
The linearized approach taken here cannot tell us what the ultimate fate of the exponentially growing modes will be. Since the modes have k ⊥ = 0, they do not respect the planar symmetry of the background. This suggests that this symmetry will be badly broken once the fluctuations become large. The treatment of the fully nonlinear field evolution is the subject of two companion papers [13] and [14]. We will explore nearly planar symmetric domain wall and nearly SO(2,1) symmetric bubble collisions using high resolution massively parallel scalar field lattice simulations. These investigations will demonstrate that the nonlinear evolution leads to a complete breakdown of the original symmetry of the background, including a dissolution of the walls and production of a population of oscillons in the collision region. This entire process is unique to more than one spatial dimension and is a completely new effect that has not previously been considered in either domain wall or bubble collisions.
of the system. We assume the field profile takes the form of an interacting kink-antikink pair This ansatz ignores production of radiation, excitations of the shape mode, and distortion of the kinks due to their mutual interaction. However, as emphasized in the main text, we wish to separate out the effects of the actual collision from the subsequent evolution and approximation (A.1) does this. Our goal is to obtain an effective Lagrangian for the locations, r, of the kink and antikink. This should be of the form for a pair of relativistic point particles interacting through a potential, along with corrections induced by the finite thickness of the kinks. For simplicity, we assumeγ = 0. Terms involvingγ arise only from the kinetic term for the fields, and we are dropping a finite thickness correction by ignoring them. Substituting (A.1) into the Lagrangian for the field, we find The required integrals are obtained by considering C f (z)dz for f (z) = zsech n (z + r)sech n (z − r) and the contour C given by α coth(2α)(12 − 20 coth 2 (2α)) − 8 3 + 10 coth 2 (2α) .

(A.4)
Note that the interactions between the kink and antikink depend on their relative speeds as well as on their separations. This is because these interactions are generated by integrals of the overlap between the kink and antikink. As their speeds increase, the kinks Lorentz contract which changes the amount of overlap. When the kink and antikink are far apart this overlap is exponentially small. Therefore, we make an additional approximation: when performing the integrals we keep overall γ multipliers only on the terms with no r dependence. The effective Lagrangian becomes with an the effective potential and The multiplier on the relativistic kinetic term can be identified as 2M where M = 2 √ 2 3 is the energy of a single kink at rest. This potential vanishes exponentially fast for ωr 1. We only consider bound motions in this paper, so we must have E − 4 √ 2 3 < 0. Therefore, at large r the walls move nonrelativistically and we can set γ ≈ 1. This is incorrect for ωr 1, but for small separations the kink and antikink are close to each other and interacting strongly. In the small separation regime, the kink profiles will be deformed, making our ansatz invalid. Setting γ = 1, the resulting effective potential, given by  Figure 17. The effective potential V ef f (r) and noncanonical contribution to the kinetic term K(r) for our effective single-particle Lagrangian describing the separation of the kink and antikink pair. Also included are the individual contributions from the gradient energy (V grad ) and potential energy (V pot ) in the original scalar field Lagrangian. For comparison, we have also included the asymptotic potential V (r) Energy conservation implies a minimum value for r, but this equation allows r → −∞. We cure this by cutting off the logarithmic divergence: which enforces the condition r(T /2) = r min where V ef f (r min ) = V ef f (r max ) and r min < 0.
In figure 18 we compare the accuracy of this analytic approximation to a full solution of the equations for the Lagrangian (A.5). We have approximated γ = 1 in K(r) and V ef f (r) but have otherwise included both the noncanonical kinetic correction and relativistic corrections. Our approximation (A.11) is very accurate for most of the evolution. The above procedure can be generalized to the asymmetric well. However, the main effect of setting δ = 0 is to break the Z 2 symmetry so that the two potential minima are no longer degenerate. This introduces a contribution to the effective potential of form [V (φ f ) − V (φ t )]r for ωr 1, which induces a constant force driving the wall and antiwall towards each other.

B Numerical Approach and Convergence Tests
In this appendix we summarize the numerical techniques used in this paper. For timeevolution in our codes, different schemes were used for the nonlinear background (3.2) and the linear fluctuations (3.3). For the one-dimensional nonlinear wave equation (3.2), we used a 10th order accurate Gauss-Legendre quadrature based method. This is a specific choice of an implicit Runge-Kutta process. Given an initial condition y t to dy/dt = H(y) at time t, the solution at time t + dt is obtained by solving where a ij and b i are numerical constants defining the process. For Gauss-Legendre methods, a ij and b i 's are solutions to The c i 's are the roots of P ν (2c − 1) where P ν (x) is the Legendre polynomial of degree ν. These relations arise by approximating the time-evolution integrals using Gauss-Legendre quadrature. Because of the excellent convergence properties of quadrature approximations, the result is an order 2ν integrator. 9 Explicit formulae for ν up to 5 are given in Table 2 of Butcher [55], but it is easier in practice to simply solve (B.3) numerically. For the linear fluctuation equation (3.3), we employ Yoshida's [56] operator-splitting technique that was introduced into the preheating community by Frolov and Huang. For further details see, for example, [56][57][58]. For this set of integrators, the solution to df /dt = H(f ) is first written as f (t+dt) = e Hdt f (t), where H should now be interpreted as an operator acting on f . We decompose H = i H i so the action of each individual H i on f is simple to compute accurately. Finally, we re-express the time evolution operator U ≡ e Hdt as a product of exponentials for the individual H i operators, e Hdt = U (w M )U (w M −1 ) . . . U (w 0 )U (w 1 ) . . . U (w M )+ O(dt n+1 ), where U (w i ) ≡ e w i H 1 dt/2 e w i H 2 dt e w i H 1 dt/2 is a second-order accurate time-evolution operator for time-step w i dt and we have specialized to the case of an operator split into only two parts for simplicity. Via clever choices of the number and value of the numerical coefficients w i , integrators of various orders n may be constructed. For this paper, we use an O(dt 6 ) method with coefficients w i given by The Gauss-Legendre and Yoshida approaches are both symplectic integrators for Hamiltonian systems. Therefore we find it convenient to use Hamilton's form for the evolution equations. With the exception of the collective coordinate location for the bouncing walls in the double well, all of the Hamiltonians used in this paper can be split into two exactly solvable pieces so that H = H 1 + H 2 . For reference, we provide these splits below, even though they are not required for the Gauss-Legendre method used to solve the nonlinear background wave equations. For the planar walls the split terms are (up to an overall normalization) The linearized fluctuations evolve in the Hamiltonian The Hamiltonian for the SO(2, 1) invariant bubbles is In all three cases, π f,i represents the canonical momentum for field f at lattice site i. The operator G[φ i ] is a discrete approximation to (∂ x φ(x i )) 2 . We now describe the spatial discretization of the system. For all production runs we used a Fourier pseudospectral approximation for the field derivatives. The only derivative appearing in the various equations of motion is the one-dimensional Laplacian along the collision axis ∂ xx . Therefore, in practice the system was evolved in real space, with the Laplacian evaluated in Fourier space through the use of the FFT. Although the resulting FFT and inverse FFT are numerically more expensive than a finite-difference approximation, the continuum limit is approached much more rapidly as seen in figure 20. This is especially important when computing Floquet exponents, as our approach requires solving 2N individual PDE's in order to form the fundamental matrix where N is the number of grid points. As well, in order to maintain a fixed accuracy in the time-integration, the ratio dx/dt should be kept constant meaning that the total work required scales as N 3 for a finite-difference approximation and N 3 log(N ) for a pseudospectral approach. 10 Since the pseudospectral calculations converge using much fewer lattice sites than the finite-difference stencils, the pseudospectral approach ends up requiring less CPU time yet is (orders of magnitude) more accurate.
To provide independent verification of our results, we also performed several runs using finite-difference discretizations of G[φ i ]. The Hamiltonian was discretized directly, thus ensuring a consistent discretization of ∇ 2 φ and (∇φ) 2 . We tested with both second-order accurate and fourth-order accurate (B.9) finite-difference stencils, where dx is the lattice grid spacing. The corresponding Laplacian stencils L[φ i ]/dx 2 satisfying i G[φ i ] + φ i L[φ i ] = 0 (on periodic grids) are then the familiar second-order accurate and fourth-order accurate centered differences.

B.1 Convergence Tests
The combination of high-order time-integrations and spectrally accurate derivative approximations leads to a rapid convergence of both the nonlinear field evolution used to study the background dynamics and the Floquet exponents obtained by solving the perturbation equations. Several measures of convergence are shown in figure 19 for the nonlinear wave equation with initial conditions as in the bottom right panel of figure 3. In the top row we show the pointwise convergence of the solution as we vary the number of grid points N (or equivalently the grid spacing) and the time-step dt, thus independently testing our Fourier spatial discretization and our Gauss-Legendre integrator. We consider two closely related measures, where φ (p) denotes the numerical solution for the pth approximation (here either the number of grid points or the time step), and we take N (p+1) /N (p) = 2 = dt (p) /dt (p+1) . To compare solutions with different spatial resolutions, we took all sums over the grid from the N = 2048 simulation. The top left panel, shows that the solution rapidly converges (to the level of machine precision) as we increase the spatial resolution, exactly as we expect for a properly resolved pseudospectral code. The top right panel shows that the growing error at late times is due to errors in the time-stepping rather than in the spatial discretization. This is not really a time-stepping issue, but a demonstration that making a pointwise comparison of the fields is not necessarily the best measure of convergence. In particular, the errors that accumulate at late time occur because we have an oscillating localized blob of field. Small errors in the oscillation frequency then lead to errors in the instantaneous value of the field. These errors manifest themselves as accumulating changes in the pointwise differences at late times. As a further test of our time evolution we check the conservation of the energy density ρ = T 00 = φ2 2 + (∇φ) 2 2 + V (φ) and field momentum P x = T 0x = −φ∂ x φ for a range of time steps dt. T µν is the energy-momentum tensor of φ and · denotes a spatial average. For all choices of time step, the field momentum is conserved to machine precision, while for dt ≤ dx/5 the energy is also conserved to machine precision.
To understand the accuracy of our Floquet exponents, and to demonstrate the great gains in accuracy obtained via a pseudospectral approach relative to a finite-difference scheme, we show convergence plots for the maximal Floquet exponent in figure 20. Here we can directly compare the individual Floquet exponents, so we plot ∆µ (p) ≡ |µ (p+1) max − µ (p) max | .  Figure 19. Convergence of our one-dimensional lattice code for the double-well potential with δ = 1/30, mL = 1024 and various choices of grid spacing dx and time step dt. We plot the two norms defined in (B.12). The accumulating errors at late times are due to small errors in the oscillation frequency and initial phase of the oscillon that has formed at the origin. Decreasing the time step below dx/5 does not lead to a decrease in this error, suggesting that it arises due to machine roundoff.
In the bottom row we demonstrate the convergence of both energy (bottom left) and momentum (bottom right) of the system for the same choices of dt as in the top right plot. Since both (conserved) quantities are constant to very high accuracy, this demonstrates that our time-stepping has converged.
left panel for various numerical schemes. The top right panel shows the convergence rate as the time-step is varied. As expected for a sixth-order accurate integrator, the error decreases rapidly, although not quite as quickly as for the Gauss-Legendre integrator. Also of note is that the error decreases uniformly for all values of k ⊥ (except for those values that are already at machine roundoff levels) indicating that important features such as the locations of stability bands where µ max = 0 are not shifting around as the time-step is varied. In the bottom row we show similar convergence plots as the number of grid points are varied while holding the simulation box size fixed. We show results for a pseudospectral approximation, a second-order accurate finite-difference stencil, and a fourth-order accurate finite-difference stencil. As promised, the psuedospectral method converges much more rapidly than the finite-differencing methods. Also of note is the uniform convergence for all k ⊥ with the pseudospectral approximation, whereas the convergence is non-uniform for the finite-difference methods. Taking the far right 4th order chart as an example, there are several extreme  figure 4. In the top right panel, we show the difference in the numerically determined values of µ max T breather holding the number of grid points fixed (at N = 64) while varying the discrete time step dt. A pseudospectral derivative approximation was used for each run and a sixth-order Yoshida integrator for the time evolution. The bottom row shows the convergence properties as the grid spacing is decreased, using a pseudospectral (bottom left), secondorder finite-difference (bottom centre) and fourth-order finite-difference (bottom right) approximation for the Laplacian. In all three panels, we took dx/dt = 20 and used a sixth-order accurate Yoshida scheme. Because of the rapid convergence, the instability charts (with the exception of the N = 32 case) are visually indistinguishable from the top left panel.
spikes in the region k 2 ⊥ (1 + v −2 ) for which the difference between the N = 128 and N = 256 approximation is of order machine precision, but then rises to the 10 −3 level when comparing to the N = 512 solution. The ultimate source of these appearing and disappearing spikes is a slight shifting of the edges of the stability bands as the resolution is varied.