Hill crossing during preheating after hilltop inflation

In 'hilltop inflation', inflation takes place when the inflaton field slowly rolls from close to a maximum of its potential (i.e. the 'hilltop') towards its minimum. When the inflaton potential is associated with a phase transition, possible topological defects produced during this phase transition, such as domain walls, are efficiently diluted during inflation. It is typically assumed that they also do not reform after inflation, i.e. that the inflaton field stays on its side of the 'hill', finally performing damped oscillations around the minimum of the potential. In this paper we study the linear and the non-linear phases of preheating after hilltop inflation. We find that the fluctuations of the inflaton field during the tachyonic oscillation phase grow strong enough to allow the inflaton field to form regions in position space where it crosses 'over the top of the hill' towards the 'wrong vacuum'. We investigate the formation and behaviour of these overshooting regions using lattice simulations: Rather than durable domain walls, these regions form oscillon-like structures (i.e. localized bubbles that oscillate between the two vacua) which should be included in a careful study of preheating in hilltop inflation.

In "hilltop inflation", inflation takes place when the inflaton field slowly rolls from close to a maximum of its potential (i.e. the "hilltop") towards its minimum. When the inflaton potential is associated with a phase transition, possible topological defects produced during this phase transition, such as domain walls, are efficiently diluted during inflation. It is typically assumed that they also do not reform after inflation, i.e. that the inflaton field stays on its side of the "hill", finally performing damped oscillations around the minimum of the potential. In this paper we study the linear and the non-linear phases of preheating after hilltop inflation. We find that the fluctuations of the inflaton field during the tachyonic oscillation phase grow strong enough to allow the inflaton field to form regions in position space where it crosses "over the top of the hill" towards the "wrong vacuum". We investigate the formation and behaviour of these overshooting regions using lattice simulations: Rather than durable domain walls, these regions form oscillon-like structures (i.e. localized bubbles that oscillate between the two vacua) which should be included in a careful study of preheating in hilltop inflation.

Introduction
Inflation has proven to be a very successful paradigm in early universe cosmology. It successfully solves the horizon and flatness problems of Big Bang cosmology and at the same time provides an explanation for the adiabatic, nearly scale invariant and Gaussian primordial curvature perturbations implied by observations of the cosmic microwave background and large scale structure [1][2][3][4][5].
At the end of inflation, the universe is dominated by the inflaton field's potential energy which is converted into radiation in a process called reheating. Reheating typically happens in at least two stages: an initial "preheating" phase of non-perturbative particle production, and a later stage of perturbative inflaton decay and thermalization. Together, these processes describe the transition of the inflationary, vacuum dominated phase to the radiation dominated universe of standard ΛCDM cosmology. A detailed understanding of the reheating phase is required not only for precise predictions of the primordial spectrum (as the expansion history after inflation affects the matching of inflationary and CMB perturbations), but also for calculating the non-thermal production of e.g. a baryon asymmetry, dark matter or other relics from inflaton decays.
Despite its success in cosmology, it is not yet clear how exactly inflation is realized in terms of particle physics. One of the many interesting possibilities is hilltop inflation, also called "new inflation" (see e.g. [6][7][8][9][10][11]): the inflaton φ starts at a local maximum of the scalar potential at φ = 0 and slowly rolls towards a minimum of the scalar potential at φ = 0, spontaneously breaking some particle physics symmetry like a family symmetry [12]. The Figure 1. Inflaton potential with the points of maximum tachyonic mass (φ m ), the end of tachyonic preheating (φ i ) and the potential minimum (v). Inflation happens at the hilltop for |φ| v 2 /( √ 24m pl ) as the inflaton field slowly rolls down towards the minimum at φ = v.
initial conditions for hilltop inflation (φ starts near a local maximum) can be easily generated by a coupling of the inflaton to some matter field χ, which allows for a period of preinflation [8,9,13,14] and also provides an efficient decay channel for reheating [14]. Other proposals for generating the initial conditions include [15][16][17][18].
Preheating in hilltop inflation has been studied numerically in [19] and analytically in [20]. Both papers find that perturbations are amplified by both tachyonic preheating and parametric resonance, and that perturbations grow very large at scales corresponding to the inflaton mass. Detailed analytical calculations for the linearised perturbations are presented in [20]; however, when the perturbations eventually grow too large, the linear approximation breaks down and any further calculation must include the non-linear interactions. [19] studies preheating numerically using lattice simulations, using a slightly different (logarithmic) hilltop potential. They study the evolution of particle number densities n k and find that the growth of perturbations is terminated by the non-linear interactions when δφ 2 ∼ φ 2 .
In this paper, we study both the linear and the early non-linear stages of preheating in hilltop inflation. Beyond confirming the earlier results of [19] and [20], we study the evolution of the inhomogeneous inflaton field in position space. We are particularly interested in the possibility that the large inflaton fluctuations might push the inflaton field over the local maximum of the scalar potential towards the "wrong" vacuum, which could in principle result in the formation of domain walls. We find that such "hill crossing" indeed occurs for a wide range of parameters. However, it turns out that these regions of "wrong" vacuum do not form domain wall networks, but that they instead generate localized oscillating bubbles, which might be interpreted as oscillons [21][22][23][24][25][26][27][28][29].
The paper is structured as follows. In section 2, we introduce the hilltop inflation model studied in this paper. The linear phase of preheating is discussed in section 3 based on [20], with special emphasis on the exponential growth of perturbations which suggests that perturbations might grow large enough to push the inflaton field over the hilltop in some regions of space. However, the linear approximation breaks down for such large perturbations, so to check whether hill crossing can actually happen, we study the non-linear preheating dynamics using lattice simulations in section 4. We then summarize our results in section 5.

Hilltop inflation
In this paper, we consider the hilltop inflation potential with v m pl . This potential is plotted in fig. 1. Hilltop inflation happens near φ = 0, with φ slowly rolling down the hilltop, while the universe expands by at least N * ∼ 60 e-folds to solve the horizon and flatness problems. We choose φ > 0, which can always be arranged by a field redefinition φ → −φ. Slow-roll inflation ends shortly after η = m 2 pl V /V = −1, which happens at the inflaton field value |φ η | v 2 / √ 24, and the inflaton quickly rolls towards the minimum at φ = v. 4

Primordial perturbations
The scalar perturbation amplitude matches the observed value A s = 2.2 × 10 −9 for 2) and the predictions for the spectral index n s and tensor-to-scalar ratio r are in some tension with the most recent Planck bounds n s = 0.965 ± 0.005 (at 68% CL) [1]. The predicted value of n s can be larger if there is an additional tachyonic mass term for φ, which occurs naturally in supergravity implementations of hilltop inflation due to corrections from the Kähler potential, or if the φ 4 /v 4 in the inflaton potential is replaced by some higher power φ p /v p with p > 4.

Initial conditions from preinflation
The initial conditions for hilltop inflation (φ 0) can be naturally generated e.g. by a phase of preinflation [8,9,13,14]. During preinflation, another scalar field χ slowly rolls down its potential, while φ is driven to zero by a large mass term generated from couplings like V int = λφ 2 χ 4 [14]. Eventually, when χ reaches very small values, the potential for φ takes the form of eq. (2.1), and φ starts to slow-roll away from zero either towards φ = v or φ = −v, starting with an initial displacement from φ = 0 due to quantum fluctuations. This second slow-roll phase, which can be described as single-field hilltop inflation in φ with the potential (2.1), lasts more than N * ∼ 60 e-folds. As a consequence, any topological defects generated during the initial symmetry breaking are inflated away, and at the end of inflation the entire observable universe has either φ > 0 or φ < 0.

Preheating
Preheating in hilltop inflation has been shown to happen in two steps [20]: tachyonic preheating and tachyonic oscillations.
2. For v 10 −5 m pl , φ does not quickly relax at the minimum. Instead, it oscillates between φ 0 and φ 2 1/4 v. These "tachyonic oscillations" cause a strong growth in small-wavelength perturbations around a scale k ∼ k peak 7(m pl /v) 3/4 H 0 , which grow until the perturbations become non-linear (for v/m pl 10 −2 ) or until Hubble damping reduces the oscillation amplitude so much that the inflaton oscillations remain in the approximately quadratic part around the minimum.

Linear preheating
In this section, we discuss the phase of linear preheating after inflation. During this phase, the inflaton field is still almost homogeneous with small fluctuations δφ which can be treated as small perturbations (keeping only terms to leading order in δφ). The initial phase of preheating can be studied by solving the time evolution of these fluctuations in Fourier space.
We start with a brief review of the theory of linear preheating and its application to the hilltop inflation potentials in eq. (2.1). For these potentials, preheating starts with tachyonic preheating, followed by a period of tachyonic oscillations. We summarize the analytic results from [20] for these two phases and show numerical results for the spectra produced by linear preheating. We then discuss how these power spectra suggest that for 10 −5 v/m pl 10 −2 , the perturbations could push the inflaton field above the hilltop towards φ(x) < 0, a process which could potentially produce domain walls.

Equations of motion
The equation of motion of the inflaton field φ(t, x) in a FRW background (neglecting backreaction from the metric perturbations 5 ) is given bÿ At the end of inflation, the inflaton field is almost homogeneous, so we can expand it in terms of a homogeneous mean field φ(t) and a perturbation field δφ(t, x): As long as the perturbations remain small, we can linearise eq. (3.1) by dropping all terms of O(δφ 2 ):φ (3.3b) The partial differential equation (3.3b) for δφ(t, x) can be transformed into ordinary differential equations by transforming to Fourier space: where φ k is the Fourier mode with comoving momentum k.

Initial conditions from Bunch-Davies vacuum
The initial conditions are given by matching to the Bunch-Davies vacuum at early times: For a numerical implementation, where matching at t → −∞ is not practical, "early times" means that the Hubble constant at the time of matching must still be approximately constant, and that the inflaton mass V (φ) should be negligible compared to k 2 /a 2 .

Inflaton fluctuation amplitude and power spectrum
Given eqs. (3.4) and (3.5), one can calculate the evolution of the scalar field perturbations for each mode k. The result can be translated into position space perturbations by an inverse Fourier transformation. For the inflaton field's variance, this gives where the limits of the integral are the IR and UV cutoffs, 6 and for d = 3 spatial dimensions the spectrum P φ can be calculated as The results from linear preheating are reliable only until higher orders in δφ become important. At that point, the fully non-linear eq. (3.1) must be solved. This is usually done by solving the equation numerically on a discrete spacetime lattice, as we will discuss in section 4. For the hilltop potential (2.1), we can estimate that the linear approximation is valid approximately until δφ 2 φ 2 , or equivalently max

Tachyonic preheating
When V (φ) < 0, all modes with k/a < −V (φ) experience exponential growth. This period of tachyonic growth lasts the longest for modes with k H 0 , which are already tachyonic at the end of inflation. Higher k modes become tachyonic later, and therefore We see that for small v/m pl 10 −1 , the oscillations continue to bring φ deeply into the tachyonic region for several oscillations, which explains why the perturbations grow non-linear within a few oscillations. For large v/m pl = 10 −1 , however, Hubble damping quickly dampens the oscillation amplitude so that after a few oscillations, the inflaton field never re-enters the tachyonic region. Preheating is much less efficient in this case: the perturbations remain small and are eventually redshifted away as the inflaton field oscillates around the quadratic part of its potential.
achieve a smaller final amplitude because they grow for a shorter time.
After tachyonic preheating 7 , the power spectrum around the Hubble scale is independent of v. For aH 0 < k k m , the power spectrum falls off as P φ ∝ (k/a) −2 , and for k k m , the spectrum is just the vacuum spectrum P φ k 2 /(4π 2 a 2 ). The variance of the inflaton field is dominated by the large contribution from the Hubble scale: For very small v/m pl 10 −6 , this implies δφ 2 (x) v 2 , so the perturbations become non-linear during tachyonic preheating -the inflaton field becomes very inhomogeneous and the description in terms of a background field plus perturbations breaks down. For larger v/m pl 10 −5 , eq. (3.9) implies δφ 2 v 2 and the linear approximation should remain viable.

Tachyonic oscillations
During the subsequent phase of tachyonic oscillations, the background field oscillates around the minimum. For v/m pl 10 −1 , these oscillations take the inflaton back into the tachyonic region with φ < φ i . In this case, one can estimate the inflaton field value after j complete oscillations as These oscillations trigger strong growth of perturbations at a characteristic peak scale For small v/m pl 10 −1 , the perturbations grow by many orders of magnitude and become non-linear within a few oscillations around the minimum. For large v/m pl 10 −1 , however, the oscillations are quickly damped due to Hubble damping (see fig. 2). After a few oscillations, the amplitude decreases so much that the inflaton field never re-enters the tachyonic region, which makes the preheating much less efficient. In this case, the perturbations remain so small that we expect the linear approximation to remain valid throughout the entire preheating phase. 8 Note that the perturbations with k k peak , which are amplified most during tachyonic preheating, do not continue to grow during the oscillation phase. Even though they are amplified every time that the inflaton field rolls down the tachyonic part of the potential, this amplification is exactly cancelled by an exponential damping sustained every time the field rolls back up. The infrared part of the spectrum therefore oscillates: it is very large near the minimum at φ ∼ v, and it is negligible at φ j (i.e. near the hilltop).

Numerical results for the power spectra
Integrating the mode equations (3.4) numerically, we can calculate the power spectrum P φ (k) of the inflaton perturbations. Fig. 3 shows P φ /φ 2 for various values of v evaluated at the time when the inflaton crosses the minimum (φ = v andφ > 0), while fig. 4 shows the power spectra closest to the hilltop at the end of each full oscillation, whereφ = 0.
All of the plots show the strong growth of power for modes with k ∼ k peak with each successive oscillation. For v/m pl ≤ 10 −2 , it only takes a few oscillations for the perturbations to grow non-linear (i.e., P φ (k peak ) φ 2 ). Only for large v/m pl = 10 −1 , the perturbations remain small due to Hubble damping, as explained above.
For the infrared modes, fig. 3 clearly shows the tachyonic preheating spectrum with P φ (aH 0 ) ∼ 10 −13 m 2 pl and P φ (k) ∝ k −2 . The infrared modes do not continue to grow during the tachyonic oscillation phase: for v/m pl ≤ 10 −2 , the spectra after the j-th oscillation lie on top of each other. For v/m pl = 10 −1 , we even see that the infrared modes slowly decay due to Hubble damping and redshifting. Near the hilltop, the tachyonic preheating spectrum almost vanishes, as one can see in fig. 4. The infrared part of the spectrum actually oscillates during the tachyonic oscillation phase, growing to maximum values at φ v and shrinking to minimal values at φ φ j , whereas P φ (k peak )/φ 2 continuously grows throughout each oscillation 9 until the evolution becomes non-linear.  Figure 3. Linear power spectra evaluated at φ = v andφ > 0 for the first few oscillations (first to last oscillation: blue to red). Tachyonic preheating enhances modes with k < k max 0.435m min , with the strongest amplification in the infrared (k H 0 ), whereas the oscillations enhance modes with much larger momenta around the peak scale k peak 7(m pl /v) 3/4 aH 0 , which is denoted by a vertical dashed line. For v/m pl = 10 −1 , we do not show every oscillation, because linear preheating lasts for very many oscillations; the dark red line in this case corresponds to over 1000 oscillations.

End of linear preheating and hill crossing
Linear preheating predicts that for 10 −5 v/m pl 10 −2 , there is a rapid growth of perturbations at the scale k peak , producing a large fluctuation amplitude This could indicate either that the fluctuations become large enough to push the inflaton field over the hilltop towards φ( x) < 0, or that linear perturbation theory breaks down and non-linear interactions stop the growth of perturbations before this happens. To find out which of these options is realized, we need to go beyond the linear approximation. As we will show below, it turns out that the fluctuations can really push φ( x) to negative values at some points in space, and the potential subsequently drives φ( x) → −v in these regions. The inflaton field is then no longer well described by a homogeneous field φ plus small perturbations δφ( x). To study the formation and evolution of the regions with φ( x) < 0, we need to solve the field equations of motion non-perturbatively using numerical lattice simulations.  Figure 4. Linear power spectra normalized to the mean inflaton value φ j , evaluated at the end of each oscillation when the inflaton is closest to the hilltop (first to last oscillation: blue to red). We see that the infrared peak from tachyonic preheating is negligible near the hilltop, whereas the perturbations around the peak scale k peak (vertical dashed line) can grow amplitudes larger than φ j , indicating that on these scales, the perturbations can push the inflaton towards φ(x) < 0. For v/m pl = 10 −1 , we do not show every oscillation, because linear preheating lasts for very many oscillations; the dark red line in this case corresponds to over 1000 oscillations.

Non-linear dynamics
In this section, we present and discuss the results of numerical lattice simulations of the non-linear stage of preheating for v/m pl = 10 −2 and v/m pl = 10 −5 .
In the previous section, we have seen that for 10 −5 v/m pl 10 −1 , the inflaton fluctuations get amplified during the initial descent towards φ = v at scales k tac ∼ aH 0 . Then the system oscillates between φ v and φ > v during which time the spectrum gets amplified at scales k peak ∼ (m pl /v) 3/4 aH 0 , either until non-linearities become dominant (for v/m pl 10 −2 ) or until the oscillations become damped and redshifted by the Hubble expansion (for v/m pl 10 −1 ).
In order to study the non-linear phase of preheating, we numerically solve the non-linear equations of motionφ where∇ are gradients with respect to the comoving coordinatesx. We use a modified version of the program LATTICEEASY [32] to solve these equations on a discrete spacetime lattice. The hierarchy between k tac and k peak implies that the lattice must have much more than (m pl /v) 3/4 points per dimension to include both scales in the simulation. For small v/m pl ∼ 10 −5 , this becomes impractical due to our limited computing power.
However, from the linear analysis we expect the fluctuations at k tac to remain relatively small, as only the fluctuations at the peak scale k peak grow during the oscillations in the linear phase. At the end of the linear phase, most of the perturbations will be around k ∼ k peak . We therefore ignore the large wavelength modes at k tac and only include scales around k peak in order to resolve the dominant small distance fluctuations.
For v/m pl = 10 −2 , we could even keep most of the infrared modes, because the hierarchy between k tac and k peak is smaller. However, the linear power spectrum at k tac is very small for v/m pl = 10 −2 , so even in this case we prefer to cut off some infrared power in favour of a larger UV cutoff, which is more important to resolve the small-distance fluctuations produced by the hill crossing.
To resolve a reasonable range around k peak , we focus on 2 spatial dimensions, although we also did some lower-resolution runs in 3 dimensions (see section 4.3) and saw the same qualitative behaviour.

Initialization of the lattice
The initial time of the lattice simulation must be chosen such that the system is still well described by the linear equations (3.4). We choose to start after η = −1 but shortly before ε = 1 2 m 2 pl (V /V ) 2 = 1, with initial conditions and parameters shown in table 1. The initial field fluctuations φ 2 k lin for v = 10 −2 m pl and v = 10 −5 m pl are extracted from the linear equations (3.4). As shown in fig. 5, the power spectra around k peak at that time are mainly dominated by the vacuum contribution. The tachyonic contribution becomes important for scales k 0.1k peak , corresponding to the infrared tail of our initial spectra.
As usual, the norm of the fluctuations is initialized as a stochastic variable obeying the probability distribution [33,34] The initial Fourier components of φ and their derivatives are given by 10 −2 0.058 8.9 × 10 −10 100k peak 4096 0.83 10 −5 0.0098 2.7 × 10 −11 100k peak 4096 4.7 × 10 −3 where α + , α − , θ + and θ − are random real numbers uniformly distributed between 0 and 1. 10 In setting the derivativeφ k , we assumed that φ k ∝ e ±ikt . Although this is true only for the vacuum fluctuations and not for the infrared part k 0.1k peak of the the spectra in fig. 5, our results are not affected by this assumption. Indeed, the lattice simulations' results during the linear phase match those of the linear analysis.
The fluctuations (4.4) are then Fourier transformed to position space and a discretized version of the equations of motion is solved on a lattice with spacing δx = 2π/k uv , with k uv = 100k peak as the ultra-violet cutoff. The size of the lattice is given by L = 2πN/k uv where N is the number of points per dimension. This defines the infrared cutoff k ir = k uv /N . In section 4.2 we present results of simulations in D = 2 spatial dimensions with N = 4096 and in section 4.3 we briefly discuss simulations in D = 3 spatial dimensions with N = 256.
In order to run the simulations, we must choose suitable rescalings. We define our program units to be This corresponds to choosing LATTICEEASY rescalings A = 1, B = k peak , r = 3/2 and s = 0.

Results of lattice simulations in 2 dimensions
Based on the linear analysis of section 3.4, we expect that for v = 10 −1 m pl the dynamics remain linear throughout the oscillations phase. Our lattice simulations confirm this result, with the variance growing to a maximal value of δφ 2 10 −5 v 2 . For v 10 −6 m pl , on the other hand, preheating is expected to become non-linear already during the tachyonic preheating stage [20], which should prevent a tachyonic oscillations phase, so that no spectral peak at k peak is expected.
For these reasons, we focus on the lattice simulation results for v/m pl = 10 −2 and v/m pl = 10 −5 . In 2 spatial dimensions we can include N = 4096 points per dimension, which is sufficient to resolve the relevant scales around k peak for both of our choices of v.  [35].

Position space slices
We clearly see that overshooting to φ < 0 indeed occurs, creating small negative regions with φ < 0 which are separated by a distance of roughly λ peak = 2π/k peak . The size d − of these negative regions is significantly smaller than the typical distance λ peak , especially for v = 10 −2 m pl . This is what one would expect from overshooting due to fluctuations dominated by a specific wavelength λ peak : the regions which overshoot are those close to the maxima of this wave, which explains both their small size d − < λ peak and the typical distance λ peak .
During the early stage around the initial overshooting, the dynamics is very different for v = 10 −2 m pl and v = 10 −5 m pl . For v = 10 −2 m pl , the initial overshooting occurs while φ ∼ v. The initial negative regions are small, and they very quickly start oscillating towards φ ∼ −v. For v = 10 −5 m pl , the initial overshooting occurs close to the hilltop with φ v, with larger initial size of the negative regions and a more complicated subsequent evolution.
In both cases, we eventually end up with small overshooting bubbles (blue), roughly separated by the distance λ peak , which are connected by filaments in which the field is close to the hilltop (yellow). These filaments then vanish, and only the blue overshooting bubbles remain.
For v = 10 −2 m pl , we can clearly see that the overshooting bubbles oscillate between φ < 0 and φ > 0. These oscillations are localized and mostly in phase: we can clearly see that the blue regions vanish in the 5th, 7th and 9th lattice slice, and most of them reappear in the 6th and 8th slice at approximately their original positions. Over time, the oscillations get damped and some of the blue regions evaporate. The oscillations of the different bubbles also get out of phase over time, which is why the 9th slice has more blue bubbles than the 5th: in the 5th slice, all bubbles have just oscillated towards φ > 0 at the same time, whereas in the 9th slice, some bubbles are lagging behind (or ahead) and are still (or already) at φ < 0.
For v = 10 −5 m pl , the blue bubbles in fig. 7 do not clearly exhibit any coherent oscillations. The largest bubbles indeed tend to not oscillate but quickly fragment into smaller structures. However, some of the smaller blue bubbles oscillate between φ < 0 and φ > 0, but these oscillations are strongly out of phase so that fig. 7 mostly shows the overall decay of the oscillations. 12 We might interpret these oscillating bubbles as "oscillons". It has been shown in the literature that such oscillons can form from collapsing bubbles in double-well potentials, and that they can be very long-lived, depending on their properties (like radius and initial field displacement) [21]. While we observe that some of the oscillon-like bubbles with φ < 0 vanish within a few oscillations, our simulation does not extend to sufficiently late times to check whether the others might remain stable over a longer timescale, or whether the bubbles which disappeared continue as long-lived oscillons which are not "hill crossing". We note that recently, for a double-well potential similar to ours, oscillons (not "hill crossing" ones) have been observed which remain long-lived on cosmological timescales [28].
In summary, we find that hill crossing occurs and produces regions with φ < 0. However, these regions do not form domains settling at φ = −v separated from the other vacuum by domain walls, but instead they behave like localized oscillations between φ ∼ −v and φ ∼ v. While some of these oscillon-like bubbles evaporate over time, our simulation does not clearly indicate whether or not some of them may be longer-lived.

Mean, variance and fraction of negative field regions
The left plots in fig. 8 show the mean φ and perturbation amplitude δφ 2 of the inflaton field. We see the first few coherent oscillations of the background inflaton field during which the perturbations remain relatively small. When the perturbation amplitude becomes comparable to the mean, the coherent background oscillations break down and hill crossing occurs. This happens during the 7th oscillation for v = 10 −2 m pl and during the 3rd oscillation for v = 10 −5 m pl , in agreement with our expectations from the linear analysis. For v = 10 −2 m pl , some coherent oscillations in φ and δφ 2 remain after the hill crossing, as the negative bubbles oscillate in phase between φ < 0 and φ > 0. For v = 10 −5 m pl , the averages φ and δφ 2 do not oscillate because the individual bubbles' oscillations are out of phase and therefore roughly average out.
The hill crossing is most obvious in the right plots of fig. 8, which show the ratio of points with φ < 0 as a percentage of total points on the lattice. For v = 10 −2 m pl , this plot also clearly shows both the in-phase oscillations of these points between φ < 0 and φ > 0, and the overall decay of the negative regions over time. One can also see that the oscillations gradually get out of phase: this is why the minima in this figure become larger during later oscillations.
For v = 10 −5 m pl , the right plot of fig. 8 shows no clear oscillations because the individual bubbles oscillate out of phase. However, one can still see both the initial overshooting and the subsequent decay of the overshooting bubbles.
We also see that just after hill crossing, the ratio of points with φ < 0 is larger for smaller v. We have observed this general tendency in several lattice runs, though the exact ratio varies with the random phases chosen when initializing the lattice, see section 4.1. Fig. 9 depicts the power spectra for different points in time, starting from around the time of hill crossing. We have checked that until shortly before hill crossing, the power spectra in 12 The oscillations are clearly visible in the animations containing more time slices, which are available online [35]. the lattice simulation match the results from the linear analysis.

Power spectrum
Around the time of hill crossing, we can see that the non-linearities wash out the oscillatory features in the spectrum for k > k peak , and that they transfer power from k peak to larger k. For example, hill crossing transfers some power by producing small overshooting regions of size d − < λ peak from perturbations of a longer wavelength λ peak . As explained above, this is expected because overshooting should happen near the local maxima of the dominant perturbation modes, and such regions around the local maxima are necessarily smaller than the wavelength.
For v = 10 −2 m pl , the spectrum remains peaked at k peak , whereas for v = 10 −5 m pl , the spectrum develops much more power at small scales k > k peak . This fits well to the lattice slices in figs. 6-7, where for v = 10 −2 m pl the bubbles roughly retain the size they had during the initial overshooting, whereas for v = 10 −5 m pl the bubbles tend to fragment and smaller features develop.
Eventually, the spectrum flattens over the k-range covered by the lattice (especially for v = 10 −5 m pl ). This limits the maximum time for which our lattice calculations remain valid. Extending the lattice simulation beyond that time would require more lattice points N , and therefore much more computing power.

Results of lattice simulations in 3 dimensions
In order to confirm the validity of our results in 3D, we repeated the simulations for v = 10 −2 m pl and v = 10 −5 m pl with N = 256 3 points. Because of the lower resolution, we decreased the ultraviolet cutoff to k uv = 50k peak . The infrared cutoff is therefore k ir 0.2k peak . The initial conditions for the field mean and fluctuations are the same as in the 2D simulations, see table 1 and fig. 5.
The 3D simulations show the same qualitative behaviour as the 2D simulations. Hill crossing happens shortly after the end of the linear phase, and we reproduce fig. 8 throughout the coherent oscillations and the initial overshooting phases. For v = 10 −2 m pl , the 2D results are reproduced until the end of the simulation at t = 60k −1 peak . On the other hand, for v = 10 −5 m pl , the simulation is only reliable until t ∼ 35k −1 peak , because afterwards the spectrum is dominated by the smallest scale fluctuations in the lattice.
Time slices of the 3D simulations are shown in fig. 10. 13 For both v/m pl = 10 −2 and 10 −5 , the initial hill crossing phase is characterized by filament-shaped regions where φ < 0. Eventually the filaments disappear and are replaced by localized oscillations between φ ∼ −v and φ ∼ v, with typical separation λ peak = 2π/k peak , as we saw in the 2D simulations.
Furthermore, for v = 10 −2 m pl , we can see a tendency of most overshooting regions to develop towards a spherical shape. For v = 10 −5 m pl , the overshooting regions are too small for their shape to be accurately resolved in the lattice in 3D.
Recapitulating, the 3D simulations confirm the results of the 2D simulations. The linear and initial hill crossing phases are accurately reproduced for both v/m pl = 10 −2 and 10 −5 .
For v = 10 −2 m pl the 3D simulation remains well-resolved until the final time of the 2D simulation and we see the overshooting regions form and tend towards a spherical shape. For v = 10 −5 m pl , the resolution becomes insufficient shortly after the initial hill crossing.

Summary and conclusions
In this paper, we have studied preheating after "hilltop inflation", where inflation takes place when the inflaton field rolls slowly from close to a maximum of its potential (i.e. the "hilltop") towards its minimum. When the inflaton potential is associated with a phase transition, possible topological defects produced during this phase transition, such as domain walls, are efficiently diluted during inflation. It is typically assumed that they also do not reform after inflation, i.e. that the inflaton field stays on its side of the "hill" (where without loss of generality φ > 0), finally performing damped oscillations around the minimum of the potential.
To investigate the preheating stage, we first integrated the mode equations during the linear phase of preheating, during which inhomogeneities are still small enough to be described as small perturbations on a homogeneous background. When the inhomogeneities get larger, i.e. during the subsequent non-linear stage, we solved the inflaton field's classical equations of motion with lattice simulations in 2D and 3D. We were particularly interested in the possibility that large inflaton fluctuations during preheating might push the inflaton field over the local maximum of the scalar potential towards the "wrong" vacuum, which could in principle result in the formation of domain walls.
We found that for 10 −5 v/m pl 10 −2 the fluctuations of the inflaton field during the tachyonic oscillation phase indeed grow strong enough to allow the inflaton field to form regions in position space where it crosses "over the top of the hill" towards the "wrong vacuum" where φ < 0. The negative regions with φ < 0 are separated by a distance of roughly λ peak = 2π/k peak , as one would expect from "hill crossing" due to fluctuations dominated by a specific wavelength λ peak . Rather than forming durable domain walls, these regions develop into localized bubbles that oscillate between φ ∼ −v and φ ∼ v, which might be interpreted as oscillons [21][22][23][24][25][26][27][28][29]. These oscillon-like bubbles should be included in a careful study of preheating in hilltop inflation.
As next steps beyond the work presented here, it will be interesting to study the evolution of these oscillon-like structures in more detail, to study their fate (and lifetime) after the hill crossing phase. Furthermore, it will be interesting to study the cosmological consequences of the hill crossing phase, especially when a second "matter field" is included which couples to the inflaton and receives its mass dynamically from the vacuum expectation value of the inflaton field (as e.g. in [14]). This may lead to explosive matter particle production at the hill crossing regions, but may also trigger a faster decay of the oszillating bubbles. Finally, in order to derive precise predictions of explicit hilltop inflation models it will be crucial to reliably calculate the expansion history of the universe, which depends on the evolution and lifetime of the hill crossing regions. v = 10 −2 m pl , N = 1024, D = 2, k uv = 100k peak . The simulation clearly shows the initial overshooting, formation and dissipation of a filament structure, and eventually localized bubbles oscillating in phase between φ ∼ −v and φ ∼ v, as discussed in section 4.2.1. A higher time resolution movie of the evolution can be downloaded at [35].  fig. 8, where the slices correspond to the blue dots in the lower right plot. The legend is shown in units of v. The simulation clearly shows the initial overshooting, formation and dissipation of a filament structure, and eventually localized bubbles with φ < 0 which decay over time, as discussed in section 4.2.1. Many bubbles also oscillate between φ ∼ −v and φ ∼ v (especially the smaller ones), though their oscillation is out of phase and thus only visible with a higher time resolution. A higher time resolution movie of the evolution can be downloaded at [35].  δφ 2 and the mean φ in units of v as a function of k peak t for v = 10 −2 m pl (above) and v = 10 −5 m pl (below). We can clearly see how the mean performs six or two oscillations, respectively, before the variance grows to values δφ ∼ φ and non-linearities become important. This matches the expectations from the linear analysis, see figs. 3 and 4. Right: The fraction of lattice points where φ(x) < 0 as a function of k peak t for v = 10 −2 m pl (above) and v = 10 −5 m pl (below). The blue dots correspond to the slices in figs. 6 and 7. We clearly see that around the onset of the non-linear stage, hill crossing happens for both values of v. For v = 10 −2 m pl , the crossing regions oscillate (initially in phase) between φ ∼ −v and φ ∼ v, which makes the ratio of points with φ < 0 oscillate. Over time, the amplitude of this oscillation is damped, as some of these oscillating bubbles decay. They also get out of phase, which explains why the minima of the ratio of negative points get larger during the later oscillations. For v = 10 −5 m pl , the oscillations are less pronounced: some of the smaller bubbles oscillate, larger bubbles tend to fragment into smaller substructures. In this case, the bubbles' oscillations are also out of phase already initially, which is why the oscillations are not clearly visible in the plots (which average over all of the many bubbles in our lattice).   For v = 10 −2 m pl (above), the spectra are shown at times 40, 44, 47 and 60 (blue to red) in units of k −1 peak . The spectra are progressively shifted towards the ultraviolet, although the tail is still dominated by the vacuum fluctuations at the end of the simulation. We conclude that the relevant scales are well resolved throughout the entire simulation.
For v = 10 −5 m pl , the spectra are shown at times 25, 26, 28, 30, 32 and 40 (blue to red) in units of k −1 peak . The ultraviolet tail progressively grows and by the end of the simulation the peak of the spectrum is only one order of magnitude larger than the ultraviolet tail, which limits our ability to study the late-time behaviour of the overshooting bubbles. However, during most of the evolution the neglected UV power is much smaller than the power at the well-resolved scales around k peak . v = 10 −2 m pl , N = 256, D = 3, k uv = 50k peak . v = 10 −5 m pl , N = 256, D = 3, k uv = 50k peak . Figure 10. φ(x) at fixed times from the 3D lattice simulations for v = 10 −2 m pl (above) and v = 10 −5 m pl (below). The contours correspond to points where φ = 0 (yellow) and φ = −0.8v (blue). For v = 10 −2 m pl , the left box corresponds to t = 45.3k −1 peak and 5.8% of points with φ < 0, and the right box to t = 53.7k −1 peak and 3.2% of points with φ < 0. The initial filament-shaped negative regions evolve into bubbles with a tendency towards a spherical shape. As in the 2D simulations, these bubbles perform localized oscillations between φ ∼ −v and φ ∼ v. The oscillations are visible in the animation containing more time slices, available online [35].
For v = 10 −5 m pl , the left box corresponds to t = 26.7k −1 peak and 10.7% of points with φ < 0, and the right box to t = 30k −1 peak and 4.4% of points with φ < 0. The initial filament-shaped negative regions are larger than for v = 10 −2 m pl . They fragment to smaller structures in a manner similar to our 2D simulations; the right box matches the corresponding 2D result, which is roughly the middle slice in fig. 7. However, due to the limited resolution in 3D we cannot resolve the subsequent evolution well enough to ascertain that they form oscillating bubbles as they did in 2D.