Semiclassical description of soliton-antisoliton pair production in particle collisions

We develop a consistent semiclassical method to calculate the probability of topological soliton-antisoliton pair production in collisions of elementary particles. In our method one adds an auxiliary external field pulling the soliton and antisoliton in the opposite directions. This transforms the original scattering process into a Schwinger pair creation of the solitons induced by the particle collision. One describes the Schwinger process semiclassically and recovers the original scattering probability in the limit of vanishing external field. We illustrate the method in (1+1)-dimensional scalar field model where the suppression exponents of soliton-antisoliton production in the multiparticle and two-particle collisions are computed numerically.


Introduction
An intriguing possibility of observing nonperturbative phenomena in particle collisions and in the early Universe triggered development of powerful semiclassical methods for description of false vacuum decay [1,2,3], instanton-like transitions [4,5] and multiparticle production [6,  7]. Although these processes occur with exponentially small rates [8,9,10,11,12,13,14] in weakly interacting theories, they may be important for baryogenesis [15,5], stability of the Higgs vacuum [16], or violation of the baryon and lepton numbers at future colliders [17]. We consider nonperturbative transitions of a new type: classically forbidden creation of topological soliton-antisoliton pairs in N -particle collisions, see the sketch of the N = 2 process in Figs. 1a,b. At weak coupling the probabilities of these transitions are expected to be exponentially suppressed. Indeed, the solitons can be crudely regarded [18] as bound states of N S ∝ 1/g 2 particles, where g 1 is the coupling constant. The cross section of their pair production in the collision of few particles involves the multiparticle factor [6] P few ∼ g 4N S (2N S )! ∼ exp(−const/g 2 ) due to roughly (2N S )! tree-level diagrams exemplified in Fig. 1c. This argument, however hand-waving 3 , suggests a general expression for the probability P N (E) ∼ e −F N (E)/g 2 (1.1) of the inclusive process N particles → solitons + particles; here F N (E) is the multiparticle suppression exponent. The suppression should disappear at sufficiently large number of colliding particles N 2N S when creation of the solitons proceeds classically. Presently the form (1.1) is supported by unitarity arguments [9] and recent calculation [19] in a quantum mechanical model describing the soliton moduli space. However, no reliable field theoretical method for evaluating the exponent F N (E) is known. Below we propose a general semiclassical method of this kind which is applicable, in particular, to the processes with N = 2 initial particles. Recently several dynamical mechanisms for enhancing the rate of soliton-antisoliton production in particle collisions were proposed [20]. They are supported by classical arguments which cannot be directly extended to quantum level. Our method for computing the exponent in Eq. (1.1) will be valuable for tests of these mechanisms and their phenomenological applications.
So far the processes in Figs. 1a,b eluded semiclassical treatment. The reason can be traced back to the attraction between the topological soliton and antisoliton [21]: taken at rest, they accelerate towards each other and annihilate 4 into many particles. Thus, the potential barrier between the soliton-antisoliton pair and multiparticle states is absent. As a consequence, classically forbidden transitions from particles to solitons cannot be directly described by the powerful semiclassical methods developed for tunneling processes.
In this paper which extends Ref. [22] we solve the above difficulty by relating production of the solitons from particles to their Schwinger pair creation in external field. We illustrate our method with explicit numerical calculations and present details of the numerical techniques.
The above relation is established in the following way. By definition the topological soliton and antisoliton can be equipped with the topological charges ±q; we denote by J µ the respective topological current, ∂ µ J µ = 0. We introduce an external U (1) field A ext µ coupled to this current 5 , where d is the number of spacetime dimensions. We consider the field A ext µ = (−E i x i , 0) with the constant "electric" component E which pulls the soliton and antisoliton in the opposite directions, see Fig. 2. Then the solitons are Schwinger pair produced like the opposite charges in the ordinary electric field. In our method one starts from the well-known semiclassical solution describing spontaneous creation of the soliton pair at E = 0. Then one adds N particles with energy E to the initial state and obtains solutions for the collision-induced 4 We do not consider integrable models where annihilation of the soliton and antisoliton may be forbidden by conservation laws. 5 The interaction in Eq. (1.2) does not need to be consistent at the quantum level. It is introduced as an auxiliary semiclassical tool and will be switched off in the end of the calculations.
Schwinger processes, cf. [23]. Finally, one raises the particle energy above the kinematic threshold 2M S of soliton-antisoliton production, where M S is the soliton mass, and takes the limit of vanishing external field E → 0. In this way one arrives at the semiclassical solutions describing creation of the solitons from particles. The semiclassical suppression exponent F N (E) is computed using these solutions.
To be concrete, imagine that our method is applied to pair production of 't Hooft-Polyakov monopoles [24] in particle collisions. The modification term (1.2) in this case couples gauge-invariant magnetic current [25] J µ to a small external field E. One can show [26] that the latter is equivalent to external magnetic field. Then the monopole-antimonopole pairs are created in the Schwinger process in accordance with the electric-magnetic duality. One therefore obtains the semiclassical solutions relevant for the monopole-antimonopole production in particle collisions from the solutions [26,27] describing their creation in the magnetic field.
The main technical difficulty of our method is related to the semiclassical description of the collision-induced Schwinger processes. Below we use the approach of [28,11,12] which involves a family of complex semiclassical solutions satisfying a certain boundary value problem. The suppression exponent F N (E) is calculated as a functional on these solutions. In the end of the calculation we send E → 0.
An obstacle to the semiclassical method appears in the phenomenologically interesting case of N = 2 colliding particles: the two-particle initial state cannot be described semiclassically. We overcome this obstacle by appealing to the Rubakov-Son-Tinyakov conjecture 6 [28] which states that the multiparticle exponent F N (E) does not depend on the initial particle number N at semiclassically small values of the latter i.e. at N 1/g 2 , see [30,31,32] for confirmations. This means that the two-particle exponent F 2 (E) coincides with F N (E) in the limit g 2 N → 0. Calculating semiclassically 7 the latter exponent and taking the limit, we obtain F 2 (E).
We illustrate the above semiclassical method in the (1+1)-dimensional scalar field model, where the coupling constant 8 g is small, x µ = (t, x), and the scalar potential V (φ) has two degenerate vacua φ = φ ± , see Fig. 3a, solid line. The model (1.3) possesses a pair of kink-like 6 An alternative strategy involving singular solutions is suggested in [29,10,7]. 7 Recall that the semiclassical description is valid at N 1, which can hold even in the limiting region φ - static solutions interpolating between the vacua: topological soliton (S) and antisoliton (A) shown in Fig. 3b. Below we consider their pair production in the N -particle collisions. Following the general strategy outlined above, we introduce an external field A ext µ = (−Ex, 0) coupled to the topological current 9 J µ = µν ∂ ν φ, where E is negative in accordance with Fig. 2 and µν is the antisymmetric symbol in two dimensions. Integrating by parts, we rewrite the modification term (1.2) as This interaction can be absorbed in the redefinition of the scalar potential V → V + g 2 Eφ.
After that φ − and φ + become false and true vacua, V (φ − ) > V (φ + ). One concludes that the Schwinger pair production of the kink-like solitons in the external field E is equivalent [2] to the spontaneous false vacuum decay via the nucleation of true vacuum bubbles. The latter is a celebrated tunneling process with well-known semiclassical description [3].
In what follows we modify the scalar potential V (φ) giving negative energy density V (φ + ) = −δρ to the true vacuum, δρ ∝ |E|. We calculate the rate of false vacuum decay accompanied by the N -particle collision [33,10,11,14] and remove the modification δρ → 0 in the end of the calculation.
Let us visualize the potential barrier between the soliton-antisoliton pair and multiparticle sector of the modified system. Kink-like soliton and antisoliton at distance l attract each other with Yukawa force F att ∝ e −m + l , where m 2 + = V (φ + ), cf. [21]. Besides, they are pushed apart by the pressure δρ. This leads to the potential energy of the static soliton-antisoliton g 2 N 1. 8 Field rescaling φ → gφ brings this constant in front of the interaction terms in the action. 9 The related topological charge is pair shown in Fig. 3c. The barrier U (l) separates the multiparticle sector of the theory from the static soliton-antisoliton pairs with l > l cr , where l cr ∼ |log(const · δρ)|/m + is the critical distance corresponding to the barrier top. The pairs with l < l cr annihilate and therefore should be attributed to the multiparticle sector. Unstable static solution φ cb (x) describing the soliton and antisoliton at a distance l = l cr is the famous critical bubble [1,2,3]; its energy E cb determines the height of the potential barrier between the solitons and multiparticle states. As δρ → 0, the critical bubble turns into infinitely separated soliton and antisoliton at rest: l cr → +∞ and E cb → 2M S . In this limit the potential barrier is hidden below the kinematic threshold E = 2M S of soliton pair production. False vacuum decay at E = N = 0 (no initial particles) is described by the bounce solution [3]. In the main body of the paper we numerically 10 find similar semiclassical solutions φ s (t, x) for the decay accompanied by the collision of N particles at energy E. We observe several regimes of transitions summarized in Fig. 4a. The shaded region E < N m − , where m − is the particle mass in the false vacuum, is kinematically forbidden. In the opposite case of high energies and N 4/g 2 (region III in Fig. 4a) transitions proceed classically and we obtain F N (E) = 0. The latter classical regime is investigated in [34], see also [35].
At smaller multiplicities the false vacuum decay is classically forbidden, F N (E) > 0. We find that the properties of the respective semiclassical solutions φ s (t, x) are qualitatively different at energies below and above the height E cb of the potential barrier (regions I and II in Fig. 4a). Solutions with E < E cb resemble the bounce: they describe formation of large true vacuum bubbles and imply infinite suppression 11 F N (E) ∝ δρ −1 in the limit δρ → 0. Highenergy solutions from the region II in Fig. 4a have smaller spatial extent. They approach the critical bubble φ cb (x) plus outgoing waves as t → +∞. Thus, at infinitesimally small δρ they describe creation of widely separated soliton-antisoliton pairs at rest. We explicitly demonstrate that the respective value of the suppression exponent F N (E) approaches a finite value as δρ → 0.
Our numerical results for the suppression exponent of producing the soliton pair from particles at δρ = 0 are shown in Fig. 4b. The two-particle exponent F 2 (E) is obtained by evaluating the additional limit g 2 N → 0. One observes that at high energies F N (E) is finite and (almost) energy-independent which is the expected behavior [9,10,14].
The paper is organized as follows. We introduce the semiclassical method in Sec. 2 and numerical technique in Sec. 3. Numerical solutions and results for the suppression exponent are described in Sec. 4. We conclude in Sec. 5. Technical details are presented in Appendices.

The semiclassical boundary value problem
Consider false vacuum decay induced by a collision of N particles at energy E. The inclusive probability of this process equals whereÛ is the quantum evolution operator and the infinite time limits t i, f → ∓∞ are assumed. The sum in Eq. (2.1) runs over all N -particle initial states with energy E in the false vacuum and final states |f containing a bubble of true vacuum. In the approximate equality we introduced the suppression exponent F N (E), c.f. Eq (1.1). At small g 2 the action (1.3) is parametrically large and any path integral with exp(iS) in the integrand can be evaluated in the saddle-point approximation. In Appendix A we review the semiclassical method for calculating the probability (2.1): introduce a path integral for P N (E) and derive equations for its (generically complex) saddle-point configuration φ s (t, x), see Refs. [28,37,5] for details. The leading semiclassical exponent exp(−F N /g 2 ) is given by the value of the integrand at φ = φ s . Let us describe the saddle-point equations for φ s (t, x) which will be solved in the next Sections. This complex configuration should extremize the action (1.3) i.e. satisfy the classical field equation where V ≡ dV /dφ. We will consider Eq. (2.2a) along the complex-time contour in Fig. 5a corresponding to evolution of particles prior to collision (part AB of the contour), tunneling via formation of the true vacuum bubble (part BC) and expansion of the bubble to infinite size (part CD). The expected structure of the semiclassical solution along the contour is visualized in Fig. 5b. The duration T of Euclidean evolution is a parameter of the solution. The field equation is supplemented with the boundary conditions related to the initial and final states of the process. Namely, in the asymptotic future the solution should be real, due to the inclusive final state in Eq. (2.1). In the asymptotic past i.e. at t = iT + t and t → −∞, the configuration φ s should evolve linearly in the false vacuum, where ω 2 k = k 2 + m 2 − and we expect that the relevant semiclassical solutions are x → −x symmetric like the soliton pair in Fig. 3b. The waves in Eq. (2.2c) represent free onshell particles in the initial state of the process. The initial saddle-point condition relates negative-and positive-frequency components of these waves, where θ is the second parameter of the solution. In the special case θ = 0 Eq. (2.2d) implies real-valued φ s at t → −∞. The initial state in this case is inclusive i.e. optimized with respect to the multiplicity N at fixed energy E, cf. Eq. (2.2b). The semiclassical solutions at θ = 0 are called periodic instantons [38], we will consider them in the Sec. 4. At θ > 0 the solutions are complex at t → −∞. Moreover, at θ → +∞ Eq. (2.2d) reduces to the Feynman boundary condition f k = 0. The initial state in this limit is indistinguishable from the vacuum and therefore contains semiclassically small number of particles, N 1/g 2 . Equations (2.2) form the T /θ boundary value problem [28] which will be used to find the semiclassical solutions φ s . The parameters T and θ of the solutions are Lagrange multipliers conjugate to the energy E and initial multiplicity N , respectively. The latter quantities are given by the standard expressions In what follows we find φ s for all possible values of T and θ and compute (E, N ) by Eqs. (2.3). The suppression exponent in Eq. (2.1) is evaluated as a functional on the semiclassical solution φ s (t, x), where the last two terms represent contributions from the nontrivial initial state of the process. Note that the semiclassical parameter g does not enter the boundary value problem (2.2). Thus, the exponent (2.4) depends on the combinations g 2 E and g 2 N rather than individually on g, E and N , cf. Eqs. (2.3). This feature is in agreement with the Rubakov-Son-Tinyakov conjecture [28] discussed in the Introduction: the semiclassical exponent does not depend on N at N 1/g 2 if the limit exists. Then F 2 (E) is the suppression exponent of the two-particle processes. We remark that the method of Lagrange multipliers implies the following Legendre transforms for T and θ, see Appendix A for derivation. Below we will keep in mind that T and θ are proportional to the partial derivatives of F N (E). Let us comment on the delusively simple final state of the process. First, recall that the configurations φ s (t, x) should describe false vacuum decay i.e. contain the bubble of true vacuum -expanding or critical -at t → +∞ (pink/gray region in Fig. 5b). This property is nontrivial and will be used as a selection rule for the relevant semiclassical solutions. Second, one might think that the semiclassical solutions are real at the real time axis: starting from real φ s and ∂ t φ s in the asymptotic future and solving the classical field equation backwards in time, one arrives at real φ s (t, x) at t ∈ R. This simple logic is, however, not applicable if the solution approaches the critical bubble φ cb (x) in the asymptotic future [39]. Indeed, the latter is unstable and contains exponentially growing and decaying modes δφ ± (t, x) ∝ e ±|ω − |t in the spectrum of linear perturbations. At large positive times one obtains φ s ≈ φ cb + A − δφ − + real waves, and the coefficient A − is complex in general. Then the overall solution, although complex-valued at the real time axis, satisfies Eq. (2.2b) asymptotically. We will demonstrate that the semiclassical solutions with E > E cb approach the critical bubble in the infinite future. In this case the asymptotic reality (2.2b) cannot be imposed at finite real t; we overcome this difficulty using the methods of [39,40,32].
To summarize, we arrived at the practical method for evaluating the suppression exponent of soliton pair production in particle collisions. One starts by describing the induced false vacuum decay at δρ > 0: finds a family of the solutions φ s (t, x) to Eqs. (2.2), relates their parameters (T, θ) to (E, N ) by Eqs. (2.3) and computes the suppression exponent F N (E) using Eq. (2.4). The exponent of the soliton-antisoliton production is recovered in the limit δρ → 0 above the kinematic threshold E > 2M S . Besides, in the limit g 2 N → 0 one obtains the exponent F 2 (E) of the two-particle processes, see Eq. (2.5).

Choosing the potential
Nonlinear semiclassical equations (2.2) do not admit analytic treatment and we solve them numerically. To this end we specify the scalar potential in dimensionless units, and a = 0.4. This function has a double-well form with false and true vacua at φ = φ − ≡ −1 and φ = φ + > 0, respectively. We have V (φ − ) = 0 and fix the energy density V (φ + ) ≡ −δρ of the true vacuum by tuning the parameter v. In particular, at v ≈ 0.75 the vacua are degenerate; the respective scalar masses are m − ≈ 1, m + ≈ 7.6. Below we start from the solutions at δρ > 0 (larger v), then send δρ → 0. Function (3.1) is shown in Fig. 3a at δρ = 0 and δρ = 1 (solid and dashed lines, respectively). The soliton and antisoliton configurations in the case of degenerate vacua are demonstrated in Fig. 3b. Let us explain the choice (3.1) of the scalar potential. First, we do not consider the standard φ 4 theory because evolution of the kink-antikink pairs there is known to exhibit chaotic behavior [41]. That chaos is related to the fact that the spectrum of linear perturbations around the φ 4 kink contains two localized modes representing its spatial translations and periodic vibrations of its form. The modes accumulate energy during the kink-antikink evolution which is therefore described by two 12 collective coordinates. This evolution is chaotic like the majority of two-dimensional mechanical motions. Irregular dynamics is difficult [42,43,44,45] for the semiclassical methods, we avoid it by changing the form of the scalar potential. We checked that the model (3.1) is not chaotic 13 , i.e. there is a single localized mode in the spectrum of linear perturbations around the soliton.
Second, the potential (3.1) is almost quadratic at φ < 0 and essentially nonlinear at positive φ. This property is convenient for numerical methods. Indeed, the initial conditions (2.2c), (2.2d) should be imposed in the asymptotic past where φ s (t, x) < 0 describes wave packets linearly evolving in the false vacuum (diagonal lines in the left part of Fig. 5b). In the model (3.1) the wave packets remain free almost up to their collision point, and the initial condition can be imposed at relatively small negative times (realistically, at Re t i ∼ −(6 ÷ 7)/m − ). Long nonlinear evolution in other models can be costly for numerical 12 The other two coordinates representing center-of-mass motion and P -odd vibrations of the kinkantikink pair decouple due to momentum conservation and x → −x symmetry. 13 Another type of irregular behavior [46] appearing at m + < m − is also not met in our model. x w a v e s w a v e s Figure 6: Rectangular lattice superimposed with the schematic form of the semiclassical solution φ j,i . The time sites t j cover the contour ABCD in Fig. 5a. Solution at x < 0 is reconstructed using the x → −x symmetry.
methods. Note that this problem with slow linearization is specific to (1+1)-dimensional systems and should be absent in higher dimensions.

Discretization
We discretize the semiclassical equations (2.2) using the rectangular (N t + 3) × N x lattice in Fig. 6. The sites t j and x i of the lattice cover the contour ABCD in Fig. 5 and the interval 0 ≤ x ≤ L x , respectively. The solution at x < 0 can be restored from the reflection symmetry φ s (t, x) = φ s (t, −x). Our lattice is as close to uniform as possible. Namely, , are constant at the real and imaginary parts of the contour. Typically, we choose |∆t j | ∆x. This property will be useful for formulating the boundary conditions. Our unknowns φ j, i ≡ φ s (t j , x i ) are the (complex) values of the semiclassical solution at the lattice sites.
We replace the derivatives in the action (1.3) with the second-order finite-difference expressions taking values at the links; below we mark all link quantities with half-integer indices. We also trade the time and space integrals for the second-order trapezoidal sums, e.g.
where the first and second expressions are used for the link and site quantities, respectively. The last sum involves the "averaged" spacings ∆t j = (∆t j + ∆t j−1 )/2 and ∆t −1 = ∆t −1 /2, For completeness we present this function in Appendix B. The lattice field equations are then obtained by extremizing S(φ) with respect to the field values at the interior sites of the contour ABCD: Note that the equations at the last spatial site i = N x − 1 correspond to the Neumann boundary condition ∂ x φ = 0 imposed at x = L x . Note also that Eqs. (3.4) can be derived by writing the lattice path integral for the probability and evaluating it in the saddle-point approximation, like in the continuous approach of Appendix A. The final boundary condition (2.2b) implies that the field is real at the last two time sites, cf. Eq. (3.2). In Sec. 2 we remarked that this condition cannot be used for solutions approaching the critical bubble at t → +∞. That case will be considered separately. Discretization of the initial conditions (2.2c), (2.2d) is far less trivial. Let us first simplify the discussion and consider the system in discrete space {x i } and continuous time t. Recall that the semiclassical evolution is linear in the beginning of the process t ∼ t −1 . We therefore introduce a deviation of the field from the false vacuum √ ∆x i and leave only quadratic terms in the action. We arrive at the system of N x coupled oscillators, Here we introduced the integration constants f n , g n and shifted the time by iT to keep contact with Eq. (2.2c). By direct inspection one observes that the the index n, vectors ξ (n) i and constants f n , g * n are the lattice analogs of k, cos(kx) and f k , g * k in Eq. (2.2c), respectively. Thus, the discrete version of the initial condition (2.2d) is f n = e −θ g n . The quantities entering this condition are extracted from ψ i (t) as = δ n, n and omitting irrelevant phase factor. To obtain the lattice form of the initial condition, we discretize the time derivatives in Eqs. (3.8). The simplest way is to perform the substitutions (3.2), (3.3) in the quadratic action (3.6), and then define ∂ t ψ i at the very first time site as cf. the derivation of the initial conditions in Appendix A.
We see that f n and g n in Eqs. (3.8), (3.9) are the linear combinations of ψ i (t −1 ), ψ i (t 0 ) and their complex conjugates. This means that the initial condition f n = e −θ g n can be explicitly rewritten as a set of 2N x real equations on Discretization leaves us with (N t + 1) × N x complex field equations (3.4) and 4N x real boundary conditions (3.5), (3.10). This is precisely the required number of equations to determine (N t + 3) × N x complex unknowns φ j, i . The semiclassical equations will be solved in the subsequent Sections.
Finally, let us introduce lattice expressions for the energy, initial particle number and suppression exponent. The full nonlinear energy is directly discretized using the substitutions (3.2), (3.3). We monitor its conservation to control the discretization errors. The initial energy and particle number (2.3) in the discrete case have the form, see Appendix B. The suppression exponent F N (E) is computed using the discrete action S(φ j, i ) and quantum numbers (3.12) in Eq. (2.4); the last integral in this expression is discretized in the standard way.

Fixing the time translations
The continuous semiclassical equations (2.2) preserve time 15 shifts. Namely, φ s (t − t 0 , x) with real t 0 satisfies the equations if φ s (t, x) does. This feature, if left unnoticed, leads to a numerical artifact. Indeed, since the time translations are explicitly broken in the lattice system, the position t 0 of the numerical solution is out of direct control: it is fixed by the discretization and finite-volume effects. Depending on the lattice parameters, it may become arbitrary large and cause divergence of the numerical procedure. To cure the above artifact, we notice [12] that due to the continuous time-shift symmetry one of the lattice equations starts to depend on the others at ∆t → 0, t −1 → −∞. Indeed, in this limit the total energy of the solution is conserved. It is real due to the final boundary conditions (3.5) and, on the other hand, proportional to the sum n ω n f n g * n in Eq. (3.12). Now, consider the initial conditions f n = e −θ g n which imply, in particular, N x relations arg f n = arg g n . One of the the latter is redundant because the sum n ω n f n g * n is automatically real. We can exclude the redundant equation without affecting the solution.
Following this line of reasoning, we drop the phase of the initial condition f n 0 = e −θ g n 0 at a given n = n 0 relating only the absolute values To keep the number of lattice equations equal to the number of unknowns, we add one condition fixing the time translation invariance of the solution. Provided the set of the lattice equations is solved, the arguments of f n 0 and g n 0 will be automatically equal up to discretization errors. In realistic numerical calculations we use Eq. (3.13) for a highly populated mode with 16 n 0 = 2 or 3; the equality of phases is then satisfied with accuracy better than |arg f n 0 − arg g n 0 | 10 −3 for all our solutions. Let us now discuss the additional equation fixing the time translation invariance. It is convenient to use different conditions at energies below and above the barrier height E cb . At low energies we impose the relation [11] Lx at the point C of the time contour in Fig. 5a. Equation (3.14) places the turning point ∂ t φ = 0 of the solution, if there is one, at a given time t = t C . Even if the solution does not have turning points, Eq. (3.14) keeps the singularities of the solution (crossed circles in Fig. 5a) away from the time contour. At high energies the structure of the semiclassical solutions changes, and we use different constraint for the time translations. Namely, we fix the center-of-mass coordinate R of the wave packet at the start of the process, see Fig. 6, where ρ E is the energy density entering Eq. (3.11): In realistic calculations we take R ≈ 0.7L x to guarantee that the initial wave packet is inside the lattice range at t = t −1 .
Constraints (3.14), (3.15) are discretized in the standard way using the substitutions (3.2) and (3.3). We repeat that exchanging one complex equation f n 0 = e −θ g n 0 for the two real -Eq. (3.13) and one of Eqs. (3.14), (3.15) -one keeps the number of equations equal to the number of unknowns.

Solving the equations
Let us select the appropriate values of the lattice parameters. To this end we regard the schematic form of the semiclassical solution φ j, i in Fig. 6.
First of all, the site numbers N t , N x should be large enough to ensure small discretization errors O(∆t 2 ), O(∆x 2 ). Note, however, that the characteristic frequencies of the solutions are estimated as ω k ∼ E/N , see Eqs. (3.12). They become higher at larger energies and smaller multiplicities. In these regions the solutions are sharp and require fine lattice resolution. On the other hand, large N t and N x are costly for numerical algorithms. In practical calculations we use lattices N t × N x = 3200 × 150 and 11000 × 500 below and above E cb , respectively. We keep track of the discretization errors and stop 17 obtaining solutions whenever the relative accuracies become worse than 1%.
Next, we want to keep the nonlinear part of the solution at x, Re t ∼ 0. Then the initial time Re t i ≡ t −1 should be large negative to ensure linear evolution at the start of the process. We find that 18 Re t −1 = −(6 ÷ 8) is large enough: in this case the relative contributions of nonlinear interactions at t ∼ t −1 are smaller than 10 −3 . Besides, the spatial volume |x| ≤ L x should encompass the entire solution. Since the initial wave packets propagate inside the lightcone (diagonal lines in Fig. 6), we require L x |Re t −1 |. In practice it is convenient to fix 19 L x = 7, then choose t −1 to keep the wave packets at |x| < L x . We checked that our results are insensitive to the value of L x at the relative level 10 −4 . Finally, we explained in Sec. 2 that the solutions with E < E cb are real at the real time axis. In this case we reduce the part CD of the time contour to two sites t Nt , t Nt+1 where the condition (3.5) is imposed. At E ≥ E cb the solution becomes real-valued only asymptotically at t → +∞, and we are obliged to extend the contour to large positive times. We choose t Nt+1 ∼ |Re t −1 |. Then the waves emitted in the interaction region remain within the lattice range, see Fig. 6.
We performed several tests of the numerical procedure, see Appendix C for details. The overall conclusion is that the linearization and finite-volume effects cause relative errors smaller than 10 −3 , while the relative discretization accuracies remain below 1%.
We proceed by describing the numerical procedure to solve the set of algebraic lattice equations. Our choice of the numerical method is limited because the field equation (2.2a) is of hyperbolic and elliptic types at the Minkowski and Euclidean parts of the time contour, respectively. In addition, the semiclassical solution φ s (t, x) is complex and satisfies the boundary rather than initial conditions. The most convenient numerical technique in this case [11,12] is based on multidimensional Newton-Raphson method. To simplify notations, let us denote the complete set of lattice field equations (3.4) and boundary conditions (3.5), (3.10), (3.13), (3.15) by F l, k (φ) = 0. In the Newton-Raphson method one starts from the approximation φ (0) j, i to the solution and repeatedly refines it by finding the correction 17 See [14] for the study of the high-energy region. 18 Recall that in our units m − ≈ 1. 19 At low energies we use L x = 8 ÷ 15 due to larger sizes of the respective true vacuum bubbles.
δφ ≡ φ s − φ (0) . The latter is obtained from the set of linear equations in the background of φ (0) . After obtaining the correction, one redefines the approximation φ (0) → φ (0) + δφ and repeats the procedure. The iterations stop once δφ becomes smaller than the predetermined numerical error. The Newton-Raphson method converges quadratically [47]. Typically, the acceptable precision δφ 10 −10 is achieved in 5 − 6 interactions. However, this method is extremely sensitive to the very first approximation φ (0) : even slightly incorrect choice of the latter may cause divergence of the iterative procedure.
We therefore use a careful strategy for finding the numerical solutions. First, we obtain the "simpler" solution at some particular values of the parameters (T, θ) = (T 0 , θ 0 ). The obvious candidate is the bounce [3]. It satisfies the semiclassical boundary value problem (2.2) at T 0 = +∞, θ 0 = 0 and, on the other hand, can be computed using half-analytic methods. In the next Section we will discuss an extended one-parametric family of "simpler" solutions. Second, we apply the Newton-Raphson method to find the solution at (T, θ) = (T 0 + ∆T, θ 0 + ∆θ) using the one at (T 0 , θ 0 ) as the first approximation. At sufficiently small ∆T , ∆θ the method has to converge. Third, starting from the new solution, we change (T, θ) by a small step, again, and obtain yet another solution. We repeat this procedure until a complete two-parametric family of the numerical solutions is found.
We finish this Section with a remark that the Newton-Raphson method requires numerical solution of the sparse linear system (3.16). This is the most time-consuming part of the numerical procedure. We perform it using two alternative algorithms described in Appendix D. Our "fast" algorithm arrives at the solution in CPU time t CP U ∝ N t · N 2 x operations but has poor stability properties. We find that it works for the high-energy solutions but accumulates round-off errors at E E cb , see explanation in Appendix D. In the latter region we use stable (yet slower) algorithm involving t CP U ∝ N t · N 3 x operations. Both our algorithms are highly parallelizable and implemented at the multiprocessor cluster.

Periodic instantons
We start by considering the periodic instantons [38] -semiclassical solutions at θ = 0 and arbitrary T . The boundary conditions (2.2b), (2.2c) and (2.2d) in this case imply reality of φ s (t, x) in the beginning and at the end of the process i.e. along the Minkowski parts AB and CD of the time contour in Fig. 5a. Further insight is gained if we assume that the solutions have turning points at the corners B and C of the contour, Then the semiclassical evolution along the Euclidean part BC also proceeds with real-valued φ s (t, x). The solutions satisfying Eq. (4.1) are 2T -periodic in Euclidean time.
Physical arguments [38] suggest that all relevant periodic instanton solutions satisfy the Ansatz (4.1). Indeed, consider a somewhat different process: transition between the sectors of the false and true vacua at temperature β −1 . Its rate Γ(β) is obtained by integrating the multiparticle probability (2.1) with the Boltzmann exponent, where the prefactors are ignored. The integral (4.2) is saturated near the saddle point where we used the Legendre transforms (2.6) in the first equalities. One sees that the saddlepoint parameters θ = 0 and T = β/2 in Eq. contains an exponentially growing mode δφ − (x) with eigenfrequency ω 2 − < 0. This means that the configuration solves the field equation at small A − . The approximate solution (4.5) is periodic in Euclidean time with turning points at t = 0 and t = πi/|ω − |. It satisfies the boundary condition (4.5) and therefore represents the periodic instanton with T = π/|ω − |.
To compute numerically the approximate configuration (4.5), we solve the static field equation and find the critical bubble φ cb (x). The spectrum of its linear perturbations is given by the matrix-diagonalization routine of Sec. 3.2, see the inset in Fig. 7. Picking up the negative mode δφ − with ω 2 − < 0 and choosing A − = 0.3, we construct the approximate solution (4.5). Now, we are ready to find the entire family of the periodic instantons. We numerically solve the field equation along the Euclidean part BC of the time contour with the boundary conditions (4.1). We use the Newton-Raphson method described in the previous Section. The very first solution (Fig. 8c) is obtained at T = π/|ω − | + ∆T , where ∆T is small. In this case the configuration (4.5) serves as a zeroth-order approximation for the numerical procedure. Next, we increase the parameter T in small steps finding one solution at a time. (squares in the lower panel of Fig. 8).
One sees a clear distinction between the low-and high-energy solutions in Figs. 8a and 8c. While the latter resembles the critical bubble, the former is nearly rotationally-invariant and has large Euclidean period T . At E → 0 the solutions approach the bounce thus describing spontaneous decay of the false vacuum. The periodic instantons are absent in the opposite region E > E cb . In Appendix E we remind that the sizes of the true vacuum bubbles become infinite at small δρ justifying the celebrated thin-wall approximation [3,33,10,36,48]. The respective solutions φ s (t, x) can be obtained analytically. Their action equals and stays constant at larger periods, see Appendix E and [36]. Note that the typical values of T and 2Im S in Eq. (4.6) are proportional to 1/δρ implying infinite suppression F N (E) → +∞ at δρ → 0. This comes with no surprise, since transitions between the degenerate vacua are energetically forbidden below the threshold E = 2M S of soliton pair creation. On the other hand, at finite T and δρ → 0 the rate (4.4), (4.6) reduces to the Boltzmann probability of finding the soliton pair in the thermal ensemble at temperature β −1 = (2T ) −1 , cf. [49]. In Fig. 9 we compare the actions of the periodic instantons at different δρ (dashed lines) with Eq. (4.6) (solid line). One observes nice agreement which becomes almost perfect if one extrapolates the numerical results to δρ = 0 (empty points in Fig. 9).
Re t − Im t

Solutions below E cb
Solving the field equation backwards and forwards in time, we continue the periodic instantons to the parts AB and CD of the time contour. We thus obtain the complete solutions, see the one in Fig. 10a. At large negative times they describe wave packets in the false vacuum which collide at t ∼ iT producing expanding bubbles of the true vacuum. We compute the in-state quantum numbers (E, N ) of the solutions by Eqs. (2.3). The line of the periodic instantons is shown with filled squares in the (E, N ) plane of Fig. 11.
Starting from the known solutions at θ = 0, we find the ones with positive θ. To this end we increase the value of θ in small steps keeping T = const. At each step we numerically solve 21 the boundary value problem (2.2). An example of the solution Re φ s (t, x) with θ > 0 is shown in Fig 10b. It is complex and contains sharp wave packets at large negative time, see the color representing arg φ s .
Evaluating the quantum numbers (E, N ) of the solutions by Eqs. (2.3), we mark them with points in the initial data plane of Fig. 11. One sees that the numerical data cover the region E < E cb and N > 1.4/g 2 : at small N the solutions become sharp and require better lattice resolution. On the other hand, the high-energy region will be explored in the next Section. The suppression exponent F N (E) is computed using Eq. (2.4).
The thin-wall arguments of the previous Section suggest that the exponent is inversely proportional to δρ at small values of the latter. Thus, its Laurent series expansion starts with the singular term, In Appendix E we deduce considering dynamics of the thin-wall bubbles, see also [36]. The function (4.8) is shown by solid line in Fig. 12. It does not depend on N and decreases with energy reaching zero at in Appendix E we argue that the thin-wall approximation breaks down in that region. Presumably, this means that the limit δρ → 0 of the suppression exponent exists at E 2M S . Let us compare the numerical results for F N (E) with the thin-wall expression (4.8). Since the periodic instantons with θ = 0 have been already studied, we consider the opposite case θ → +∞ or g 2 N → 0. In this limit F N (E) coincides with the exponent F 2 (E) of transitions from the two-particle initial states, see the conjecture (2.5). We extrapolate the numerical data for F N (E) to g 2 N = 0 and obtain the dashed lines in Fig. 12; the details of this extrapolation will be discussed in Sec. 4.6. One observes that the numerical graphs approach Eq. (4.8) at smaller δρ and coincide with it after the additional extrapolation to δρ = 0 (empty points).

Going to E > E cb
Since the thin-wall bubbles do not describe transitions at E > E cb , one expects to find new properties of the semiclassical solutions in that region. In realistic calculations we observe somewhat different effect: at energies above 22 E cb our numerical method either diverges or produces unphysical "reflected" solutions exemplified in Fig. 13. The latter satisfy the semiclassical equations but approach the false vacuum at t → +∞. They apparently cannot describe transitions between the vacua. One immediately identifies [39] the root of the problems: a condition forcing the solutions to interpolate between the two vacua is not present in the semiclassical boundary value problem (2.2). As a consequence, the numerical procedure can produce unphysical solutions even if the ones with correct properties exist. The problem of fixing the qualitative behavior of the saddle-point configurations is rather general. We solve it using the -regularization technique of Refs. [39,40,32]. To this end we recall that in our numerical method the solutions with shrinking bubbles are continuously obtained from the physical ones by decreasing the parameter T below some value T * (θ). The "boundary" solutions at T = T * (θ) are very special: their true vacuum bubbles neither shrink nor expand but survive to t → +∞. The main idea of the -regularization is to exclude all "boundary" configurations from the set of accessible semiclassical solutions. Once this is achieved, one cannot obtain solutions with shrinking bubbles from the physical ones by continuous deformations. Then at E > E cb we will find solutions with correct properties (if they exist).
With this logic in mind, we add [39] a small imaginary term to the classical action, where the modification parameter and functional T int are positive-definite. Importantly, we require special properties of T int [φ]: it should take finite values on any configuration interpolating between the vacua and diverge if φ(t, x) contains a static finite-size bubble in the final state. Then the latter "boundary" configurations with static bubbles represent singularities of the modified action S → +i∞ and cannot coincide with its extrema. To the contrary, the semiclassical solutions in the modified system extremize S and do not belong to the class of "boundary" configurations. We conclude that the modified solutions cannot change qualitative properties in the course of continuous deformations; finding their continuous family at E > E cb and sending → +0, one arrives at correct high-energy solutions.
Let us construct the functional T int [φ] with the desired properties. We choose where f (x) = e −x 2 /2σ 2 f and W int (u) = u 4 e −u 2 ; a = 0.4. The function f (x) restricts the spatial integral to the central region |x| σ f , where σ f = 0.4 in our numerical calculations. At the same time, W int is vanishingly small at φ ≈ φ ± and takes positive values between the vacua, φ − < φ < φ + . Accordingly, the time integral in Eq. (4.10) diverges if the configuration φ(t, x) contains a static finite-size bubble at t → +∞. However, if the configuration is physical and the bubble expands, the value of the field at |x| σ f tends to φ + at large times and the t-integral converges. We conclude that the functional (4.10) discriminates between the "boundary" and physical configurations. Note that the modification (4.9), (4.10) simply deforms the scalar potential, introducing minimal changes of the numerical code. We remark that the regularization (4.9) is pretty general: it was successfully applied in quantum mechanics [39,40,32,45,50], field theory [12] and gravity [51]. Besides, it can be justified by adding a constraint to the path integral with the Faddeev-Popov trick [40,32], cf. [52].
Following the above recipe, we pick up some physical solution at E < E cb , introduce regularization (4.11) with = 5 · 10 −2 and find the modified solution with the same T and θ. After that we decrease T obtaining the set of modified solutions at E > E cb , see Fig. 14a and cf. Fig. 13. As expected, all these solutions contain expanding bubbles at t → +∞ and therefore describe transitions between the vacua. We finally decrease the modification parameter to infinitesimally small values 10 −3 , Fig. 14b, and continue to explore the plane of initial data by changing the parameters T and θ in small steps. Eventually, we obtain all possible solutions at E > E cb , see   Fig. 15. The initial wave packets in these solutions are composed of high-frequency modes, their true vacuum bubbles are small, cf. Fig. 10. Besides, the respective periods T of Euclidean evolutions are short. This allows us to use the faster numerical algorithm of Appendix D. Recall that the solutions become singular in the limits of small N and high E; for our best lattice this bounds the accessible region to N 1/g 2 and E 14/g 2 .
We stress that the transition to the region E > E cb is smooth at > 0. In particular, the parameters N and T of the modified solutions change smoothly along the lines θ = const, see Fig. 17. However, this crossover becomes sharper at smaller suggesting continuous rather than differentiable changes in the = 0 solutions across E cb . Now, consider the limit → +0 of the semiclassical solutions. Figure 14 displays two solutions with E > E cb at different values of the regularization parameter . The true vacuum bubble in the solution with smaller remains static for a notably longer period before it starts to expand. This suggests that the original solution with = 0 contains a static finite-size bubble in the final state (t → +∞). The latter property becomes apparent once we plot the t = const sections of the solution Re φ s (t, x) with infinitesimally small , see Importantly, we find that the asymptotics (4.12) holds for all original ( = 0) solutions above E cb , not just the ones at the boundary of this region. Thus, all high-energy solutions are unstable and cannot be obtained by direct numerical method. The unusual behavior of the solutions at E > E cb suggests that the respective processes of false vacuum decay proceed in two stages. First, the critical bubble is created with exponentially small probability; the energy excess E − E cb is dropped in the form of the outgoing waves. Second, the critical bubble, being classically unstable, decays producing a growing bubble of true vacuum with probability of order one. Similar tunneling mechanisms appear in multidimensional quantum mechanics [44,39,40,32,50] and other models of field theory [12,13].
Finally, consider the limit → +0 of the quantum numbers (E, N ) and exponent F N Re t − Im t at fixed T and θ. We find that these quantities are regular functions of : the data points in Fig. 18b are well fitted with the quadratic polynomials (dashed lines). We quadratically extrapolate them to = 0 and obtain the exponent F N (E) of false vacuum decay at high energies. The numerical errors related to this extrapolation procedure are negligibly small. The suppression exponent is plotted at g 2 N = 2 and δρ = 0.4 in Fig. 19.

Classical over-barrier transitions
An important test of our semiclassical method involves classically allowed decay of the false vacuum filled with many particles. These processes proceed with unsuppressed probabilities, F N (E) = 0, their initial quantum numbers (E > E cb , N ) occupy the upper right corner in Fig. 15. Importantly, the minimal particle number N = N min (E) required for such transitions can be obtained in two ways: by studying the real classical solutions and by finding the region F N (E) = 0 from the classically forbidden side. Comparing the two results, one checks the semiclassical method at high energies. We studied the classically allowed decay of the multiparticle states in the false vacuum in [34] using the stochastic sampling technique of [35]. Namely, we constructed many sets of random Cauchy data {φ(t i , x), ∂ t φ(t i , x)} in the false vacuum and obtained a classical solution for each set. We selected the solutions arriving to the true vacuum, i.e. those containing expanding bubbles at t → +∞. Finally, we calculated the initial quantum numbers (E, N ), Eqs. On the other hand, the function N min (E) can be obtained using the semiclassical results of the previous Section. One notes that any real solution to the classical field equations satisfies the boundary value problem (2.2) at T = θ = 0 and gives F N (E) = 0. Thus, at T, θ → +0 the semiclassical solutions become real and their quantum numbers (E, N ) approach the boundary N = N min (E) of classical over-barrier transitions. Besides, since F N (E) vanishes at the latter boundary, where the Legendre transforms (2.6) were used in the last equality. Relation (4.13) implies that the limit T, θ → 0 should be performed at ϑ ≡ θ/T = const, where ϑ parametrizes the curve N min (E). An exemplary semiclassical solution with small T and θ is plotted in Fig. 16a, its quantum numbers are marked by the circle in Fig. 15. In Fig. 20 we plot the boundary N min (E) extracted from the semiclassical results (solid line). It almost coincides with that obtained from the real classical solutions. This supports our semiclassical method in the high-energy region E > E cb . 4.5 Soliton-antisoliton production: δρ → 0 So far we have considered the collision-induced decay of the false vacuum. In this Section we send δρ → +0 and extract the exponent F N (E) of the soliton-antisoliton pair production in particle collisions. Keeping E > E cb and > 0, we decrease the energy density of the true vacuum to δρ ∼ 0.02 ÷ 0.06. We arrive at the semiclassical solutions exemplified in Fig. 21a. Their true vacuum bubbles expand at small constant velocities without any apparent acceleration. Sending, in addition, → +0, we obtain the solutions arriving at the static critical bubble φ cb (x). Recall that the spacial size of φ cb is logarithmically large at small δρ.
As expected, we find that the semiclassical exponent F N (E) is not very sensitive to δρ in the high-energy region E > E cb . Moreover, at small δρ its dependence is linear, see Fig. 21b. This means that the singular term in the "thin-wall" expansion (4.7) is absent and the limit δρ → 0 of the exponent exists. Linearly extrapolating F N (E) to δρ = 0 at fixed E and N (dashed line in Fig. 21b), we obtain the exponent of soliton-antisoliton production in particle collisions.
In Fig. 22 we plot the semiclassical exponent at g 2 N = 2: dashed lines represent the numerical data at different values of δρ, solid line is the result of extrapolation to δρ = 0. The results for F N (E) at δρ = 0 are also shown in Fig. 4b at different values of g 2 N . One observes that the exponent decreases with energy approaching constant at E 2M S . In [14] we argued that this is the expected behavior for the collision-induced tunneling at E → +∞. Note that numerical extrapolation to δρ = 0 is harder to perform near the threshold E cb = 2M S because the asymptotic expansion of F N in δρ has different form at smaller energies. However, the thin-wall arguments of Appendix E suggest near-threshold behavior where c 1 and c 2 are unknown functions of N . We find that the numerical results in Figs. 22 and 4b are consistent with the asymptotics (4.14) (dotted lines in both figures).

Two-particle processes
Now, consider creation of the soliton-antisoliton pairs in the two-particle collisions (N = 2). These processes are natural from the viewpoint of collider physics but cannot be described by direct semiclassical methods. We compute their leading exponent F 2 (E) using the Rubakov-Son-Tinyakov conjecture (2.5), i.e. evaluating the limit g 2 N → 0 of the multiparticle exponent F N (E). Recall that the semiclassical solutions develop a singularity in this limit. Thus, we cannot address them directly. In what follows we extrapolate the numerical results to g 2 N → 0 using an educated guess on the behavior of the multiparticle exponent. One finds that vanishingly small initial particle number is achieved at θ → +∞: in this case the initial condition (2.2d) reduces to the Feynman one, and the initial state of the process becomes semiclassically indistinguishable from the vacuum. Besides, since the combination e −θ enters the semiclassical equations, one expects regular expansion g 2 N = n 1 · e −θ + n 2 · e −2θ + . . . ⇒ θ = − log(g 2 N ) + θ 0 + θ 1 (g 2 N ) + . . . , (4.15) where n i and θ i are the energy-dependent Taylor coefficients. Now, we integrate the Legendre transform (2.6) at fixed energy, 16) where F N → F 2 and g 2 θN → 0 at g 2 N → 0; besides, we integrated by parts in the second equality. Substituting the expected behavior (4.15) into Eq. (4.16), we obtain the asymptotic form of the exponent at small g 2 N . This behavior, if confirmed, automatically implies that the limit g 2 N → 0 of the multiparticle exponent exists.
In Fig. 23a we compare the numerical data for the quantity F N + g 2 N θ (points) with the expectation (4.17) (dashed line). One observes a consistent fit involving two parameters, F 2 and θ 1 . Changing the number of data points in the fit, we learn that the result for F 2 (E) is stable with relative precision of order 1%.
The final numerical graphs of the two-particle exponent F 2 (E) are shown in Fig. 23b: dashed and solid lines are obtained at δρ > 0 and δρ → 0, respectively. The last graph corresponding to soliton-antisoliton production is repeated in Fig. 4b.

Concluding remarks
In this paper we developed a new semiclassical method to calculate the probability of the topological soliton-antisoliton pair production in particle collisions. Our main idea was previously reported in [22], here we presented the details, confirmations and generalizations of the method.
In spite of all technicalities, the essence of our approach is simple: it is based on expectation that the semiclassical solutions describing classically forbidden transitions continuously depend on parameters of the transitions. 23 With this idea in mind, we introduced a small external field E coupled to the topological charges of the soliton and antisoliton. We started from the well-known solutions describing spontaneous Schwinger creation of the soliton pairs in the field E. Then we added colliding particles to the initial state of the process and obtained the solutions for the collision-induced Schwinger processes. Finally, we gradually increased the collision energy E above the threshold 2M S of soliton pair production and switched off the field, E → 0. In this way we continuously arrived at the solutions describing soliton-antisoliton production in particle collisions.
We illustrated the above semiclassical method with the explicit numerical calculations in the (1 + 1)-dimensional model of a scalar field. In particular, we demonstrated that the above-mentioned semiclassical solutions continuously depend on the parameters E and δρ ∝ |E| of the induced Schwinger process. We showed that these solutions reproduce the thin-wall results at low energies and correctly describe classical over-barrier transitions at E > 2M S and high initial multiplicities. Finally, we computed the semiclassical suppression exponents of the soliton-antisoliton pair production in the N -particle and two-particle collisions, see Fig. 4b.
We believe the semiclassical approach of this paper will be useful in other models / setups of particle and condensed matter physics.
where |a , |b are the eigenstates of the false-vacuum annihilation operatorsâ k with eigenvalues a k , b k and normalization a|b = exp dk a * k b k . Integrals over T and θ run along the imaginary axes. Expression (A.1) can be proven [37] by acting on the Fock stateŝ a † k 1 . . .â † kn |0 in the coherent-state representation. We transform Eq. (A.1) into the configuration representation using where φ i (k) stands for the spatial Fourier transform of φ i (x). Here and below the prefactors are omitted. We obtain the matrix elements the quadratic functional of φ i and φ i . Next, we write the probability (2.1) in the path integral form, where the dots denote the time derivatives coming from the variations of the action, e.g.
δS/δφ f =φ f /g 2 . Finally, differentiation with respect to T and θ gives equations where B is defined in Eq.
Rewriting the last term with help of Eq. (A.9) and the boundary condition (2.2d), we arrive at the standard expression (2.4). We finish derivation with two remarks. First, the solution φ s and all equations will be considered along the complex-time contour in Fig. 5 corresponding to tunneling. In this case the initial time t i is taken complex because the exponent (2.4) does not depend on it anyway. Second, the functional (A.10) depends on E, N explicitly and via the saddle-point configuration {φ s , T, θ}. However, since we have already extremized this functional with respect to φ, T and θ, its E and N -derivatives are given by Eqs. (2.6).

B Lattice formulation
In this Appendix we list the discretized action and finite-difference semiclassical equations (2.2), see [11,12] for similar approaches.
Our lattice {t j , x i } is defined in Sec. 3.2. Performing the substitutions (3.2) and (3.3) in the action (1.3), we obtain, Recall that the lattice field equations are the derivatives of this expression at the internal points 0 ≤ j ≤ N t of the time contour. Assuming for simplicity 1 ≤ i ≤ N x − 2, we obtain, where V is the derivative of the scalar potential. Equations at the spatial boundaries, i.e. at i = 0 and i = N x − 1, are derived in similar manner. We assume that the evolution is linear in the beginning of the process. In this case the potential reduces to V ≈ m 2 − (φ − φ − ) 2 /2. Substituting it into the above action, we obtain Eq. (3.6) with discrete time and We directly compute the real eigenvectors ξ (n) i and eigenvalues ω 2 n of this symmetric matrix by QR decomposition. The coefficients f n and g n of the linear solution (3.7) are then given by Eqs. (3.8), where the time derivative (3.9) simplifies, Using explicit expressions for f n and g n , we rewrite the initial condition f n = e −θ g n in the matrix form, where the discrete operator ∂ t is defined in Eq. (B.3). One reads off the matrices M R and M I in Eq. (3.10) from the above equations. We directly discretize the nonlinear energy (3.11), It conserves, i.e. does not depend on j up to O(∆t 2 ) numerical errors. Somewhat different discretization procedure is natural for the energy of the linear system (3.6) at t = t −1 . We start from the expression in continuous time and discrete space. Substituting the solution (3.7), we obtain the first of Eqs. (3.12). The respective formula for the initial multiplicity is then easily deduced from Eqs. (2.3).
We finish this Appendix with the lattice expression for the suppression exponent 24 F N (E) = −g 2 (2ET + N θ) + 2g 2 Im S + Im where the action is given by Eq. (B.1).

C Tests of the numerical procedure
We subdued the lattice solutions to a number of consistency tests which support the numerical methods of Sec. 3 and allow us to estimate the numerical errors. For a start, we checked sensitivity of the results to the spatial cutoff L x . To this end we increased the cutoff from L x = 7 to 14 at a fixed lattice spacing ∆x = L x /(N x − 1). The integral quantities E, N and F N stayed independent of L x up to relative errors of order 10 −4 . 24 The last integral in this expression is O(∆t) and O(∆x 2 ) accurate. In practice the first-order correction in ∆t is negligible because |∆t| ∆x. All other lattice expressions in our study are second-order both in space and time.
Another source of numerical errors is related to the finite extent t −1 ≤ t ≤ t Nt+1 of the temporal lattice. Recall that the semiclassical solutions should describe free waves in the false vacuum at t ∼ t −1 . This property is conceptually important, it was used in the derivation of the initial condition (2.2d). We estimate the effect of nonlinear interactions in the beginning of the process by evaluating the energy of the linearized system given by Eqs. (3.12), (3.8) at different time sites t j and comparing it with the full conserved energy (3.11), see Fig. 24a. The graphs stay close at t ∼ t −1 ≈ −8.5, separating in the nonlinear region t −3. The relative error due to nonlinear interactions in the beginning of the process is estimated as |E full − E linear |/E full 10 −3 . One concludes that the semiclassical evolution at t ∼ t −1 is, indeed, free thanks to the clever choice of the scalar potential in Sec. 3.1.
Let us turn to discretization effects which should be O(∆t 2 ) and O(∆x 2 ) small in the second-order finite-difference scheme. First of all, one observes that conservation of the full energy in Fig. 24a is violated at the level of ∆E full /E full 10 −3 . This is the effect of time discretization because energy conservation is restored at ∆t → 0 and finite ∆x. Second, we directly estimate the finite-difference errors comparing numerical results at different lattices, see Figs. 24b and 25. Recall that our reference-point lattices are different at energies below and above E cb : N t × N x = 3200 × 150 and 11000 × 500, respectively. In the former case the relative errors in E, N and F N stay below 1% reaching maximum at the smallest N . Fine lattice resolution at high energies gives errors well below 1% at E ∼ E cb ≈ 6.1 and high N . The errors grow, however, to 1% at E 14/g 2 and/or small multiplicities. In what follows we exclude the results with g 2 E > 14 and g 2 N < 1 due to improperly high discretization errors. To conclude, we keep the finite-difference effects below the relative level of 1%. Figure 25 shows that most of the integral quantities linearly depend on ∆x 2 ∝ N −2 x and ∆t 2 ∝ N −2 t , like they should in the second-order discretization scheme. The only exception is the energy E (lower panel in Fig. 25b) which receives small nonpolynomial correction from adiabatic high-frequency waves in the solution. This effect is negligible in the region E 14/g 2 which we consider (see, however, the study of the high-energy solutions in [14]).
In Sec. 3.3 we traded the initial condition arg f n 0 = arg g n 0 at a given n 0 for the artificial constraint fixing the time-translation invariance. We argued that the omitted condition will be automatically satisfied once the other lattice equations are solved correctly. In Fig. 26 we demonstrate that this is, indeed, the case: at n 0 = 2 the related numerical error is smaller than 10 −3 .
A good piece of our qualitative and quantitative results is based on the Legendre transforms (2.6) which should hold with acceptable numerical precision. In Fig. 27 we plot the partial derivatives (∂ E F N , ∂ N F N )/g 2 of the suppression exponent and compare them to −2T and −θ (dashed horizontal lines). One observes that the Legendre transforms hold with absolute precision ∆T, ∆θ 10 −3 .
To summarize, the relative finite-volume and discretization effects in our solutions are smaller than 10 −3 and 10 −2 , respectively. The other numerical errors are vanishingly small.

D Solving the linear system
Each iteration of the Newton-Raphson method involves solution to the system of (N t + 3) × N x 10 6 complex linear equations (3.16). This is the major CPU time-consuming part of the numerical procedure. We decrease its computational cost using the sparse structure of the lattice field equations.
In what follows we suppress the spatial indices working with the N x -dimensional vector field φ j = (φ j, 0 , . . . , φ j, Nx−1 ). Recall that the linear system (3.16) is obtained by substituting φ j = φ (0) j + δφ j into the lattice equations and expanding them to the linear order in δφ j . In particular, the discrete field equations (B.2) take the form, where F j is their left-hand side at φ = φ (0) , while C j , L j and R j are the N x × N x coefficient matrices which can be deduced from Eqs. (B.2), e.g. (C j ) i, i = ∂F j, i /∂φ j, i φ (0) . Note that Eq. (D.1) is sparse: it relates δφ j−1 , δφ j and δφ j+1 . Besides, the matrices L j and R j are diagonal, while C j is three-diagonal. We will use both these facts while solving the linear system. Our "stable" algorithm [31] eliminates equations from the set (D.1). Namely, suppose we express δφ j from the j-th equation, and substitute it into the neighbouring (j ∓ 1)-th equations. The latter will keep the form (D.1) albeit with new coefficient matrices. In particular, the (j − 1)-th equation will relate δφ j−2 , δφ j−1 and δφ j+1 with Repeatedly using Eqs. (D.3), we eliminate all field equations except for the very first and last ones at j = 0, N t . To make this procedure more stable, we first apply (D.3) to the equations with odd j, then eliminate odd equations from the remaining set, etc. Once we are left with the equations at j = 0 and N t relating δφ −1 , δφ 0 , δφ Nt , δφ Nt+1 , we add the linearized boundary conditions and solve the resulting linear system by the direct method of LU decomposition. The complete solution {δφ j } is restored from the boundary values of δφ and Eqs. (D.2).
The above elimination process is apparently ineffective. Indeed, the substitutions (D.3) do not preserve the sparse form of the coefficient matrices thus requiring general matrix multiplications and inversions at the second stage of the elimination process. One concludes that cN t N 3 x operations with c ∼ 1 are needed for obtaining the solution. Note, however, that the disadvantages of this "slow" algorithm are compensated by its exceptional stability properties. We exploit it at E < E cb where the faster procedure accumulates round-off errors and causes divergence of the Newton-Raphson iterations.
Our "fast" method benefits from the sparse form of the coefficient matrices in Eq. (D.1). It is based on Cauchy problem for the lattice field equation. Consider the N x × (2N x + 1) matrix Ξ j satisfying the analog of Eq. (D.1), where the last matrix term in the right-hand side contains the vector F j in the last column; the bold 0 and 1 denote zero and unit N x × N x matrices. We solve the Cauchy problem for the "propagator" Ξ with the initial conditions This procedure involves multiplications of Ξ j by the diagonal and three-diagonal matrices L j , R −1 j and C j , i.e. cN t N 2 x operations, c ∼ 1. Given the propagator Ξ, we introduce a convenient representation of the solution, where the (2N x + 1)-vector in the right-hand side is composed from the boundary fields δφ −1 and δφ 0 . One can check that the vector (D.6) satisfies Eq. (D.1). Equation (D.6) relates, in particular, δφ Nt , δφ Nt+1 to δφ −1 and δφ 0 . Adding these two relations to the set of boundary conditions, we obtain a closed linear system for δφ −1 , δφ 0 , δφ Nt and δφ Nt+1 . We solve the latter using the LU decomposition method. Given δφ −1 and δφ 0 , we Cauchy-evolve Eq. (D.1) determining the solution {δφ j }.
In the "fast" method we arrive at the solution in cN t N 2 x operations, a factor of N x faster than in the "slow" algorithm. The N 3 x and N 2 x scalings of the CPU times in our "slow" and "fast" codes are illustrated in Figs. 28a,b.
Let us point at the reason behind the poor stability properties of the "fast" algorithm at low energies. We argued in Sec. 4 that the distinctive property of the semiclassical solutions at E < E cb is long periods T of their Euclidean evolutions. Linear perturbations δφ j grow exponentially in Euclidean spacetime magnifying the initial round-off errors. They cannot be correctly evolved within the "fast" Cauchy approach as opposed to the homogeneous "slow" algorithm. As a consequence, we exploit the "fast" procedure only for almost-Minkowskian solutions at E > E cb . Both our algorithms can be performed in parallel. To this end one divides equations (D.1) into the subsets with index ranges j = 0 . . . N k , N k . . . 2N k , etc., and performs computations in every subset at a separate processor. At the second stage of the algorithm one combines the data from all processors. For example, in the parallel version of the "fast" algorithm the propagators Ξ (k) j in all subsets are computed, and the original propagator Ξ is restored as their product. In Fig. 29 we demonstrate that the elapsed real time for obtaining Ξ scales as c/N proc with the number of processors.

E Thin-wall approximation
The semiclassical solutions describing false vacuum decay can be found analytically if the sizes of their true vacuum bubbles are parametrically large compared to the widths of the bubble walls. We will see that this thin-wall regime [3,33,10,36] occurs at δρ → 0 and E < E cb . Here A is the spacetime area occupied by the true vacuum, L is the length of its boundaries; in the second equality we expressed A and L in terms x 1 (τ ), x 2 (τ ) denoting the τ -derivatives with dots. Anticipating the limit δρ → 0, we identified the wall tension with the soliton mass M S . In what follows we treat the action (E.1) as a functional of x 1 (τ ) and x 2 (τ ). Now, let us construct the periodic instantons, i.e. solutions to the Euclidean field equations with the boundary conditions (4.1). One notes that these are the extrema of the Euclidean action in the interval τ ∈ [−T, 0]. Indeed, the Neumann conditions (4.1) mean that the action is stationary with respect to the boundary values of the fields at τ = −T and 0. Working in the thin-wall approximation, we extremize the action (E.1) by varying x 1 (τ ) and x 2 (τ ). One can draw some intuition here by noting that the latter action is similar to the static energy of a soap bubble in two space dimensions. Its extremum is achieved if the bubble walls x 1 and x 2 form circular arcs with fixed curvature radius Extremization with respect to the boundary values of x 1 and x 2 gives conditionsẋ 1 =ẋ 2 = 0 at τ = −T and 0. The only exception appears if the arcs x 1 and x 2 meet at one point as in Fig. 30b; then the boundary condition changes toẋ 1 = −ẋ 2 . The stationary thin-wall configurations at T < R b and T > R b are shown in Figs. 30b and 30c, respectively. Substituting them into the action (E.1), one obtains, where α is the half-angle between the arcs in Fig. 30b. One obtains Eq. (4.6) for the Minkowski action S ≡ iS E in terms of the half-period T . Let us find out when the thin-wall approximation is trustworthy. One naively expects this regime to be solid at small δρ because all bubble sizes are proportional to R b ∝ 1/δρ. However, the distance between the bubble walls in Fig. 30b vanishes at α → 0 as α 2 R b . Comparing it with the typical wall width m −1 + , one obtains the correct condition for the thin-wall approximation. This inequality breaks down when the solution approaches the critical bubble. Indeed, consider the thin-wall energy E = 2M S cos α , (E.4) which is easily deduced from Eq. (E.1) and Fig. 30b. It tends to 2M S as α, T → 0. Thus, at small δρ the thin-wall approximation is trustworthy below the barrier height E cb but breaks at E E cb . Next, we consider transitions from the semi-inclusive initial states with fixed energy E and multiplicity N . Note that the periodic instantons describe particular processes of this kind with N = N P I (E), see Fig. 11. Their in-states are obtained by continuing the solutions from τ ≡ it = −T to the initial part AB of the time contour in Fig. 5a. One notices [36], however, that the τ = −T sections of the thin-wall configurations in Figs. 30b,c do not depend on δρ. This suggests that the initial states of the periodic instantons and their matrix elements with the states of smaller multiplicities are also independent of δρ. Thus, decreasing the initial multiplicity from N = N P I to N = 2 one at best gets O(δρ 0 ) correction to the exponent F N (E). We conclude 25 that the leading 1/δρ part of F N (E) does not depend on N . Expressing the action (E.2) in terms of energy (E.4) and substituting it into Eq. (2.4), one obtains the thin-wall result (4.7), (4.8) for the suppression exponent at E < E cb .
Finally, let us guess the behavior of F N (E) near the point E ≈ 2M S where the asymptotic expansion in δρ breaks down. Inequality (E.3) shows that the thin-wall approximation is valid at where we used Eq. (E.4) and m + ∼ g 2 M S . One expects that the terms of the thin-wall expansion (4.7) become comparable at the boundary of this interval. In particular, where the explicit form (4.8) was used. We express δρ from Eq. (E.5), and obtain F N,0 ∼ const · |E − 2M S | 1/2 . Assuming some regular contribution besides this singular term, we arrive at Eq. (4.14).