{'fraction_non_alphanumeric': 0.0465159121263759, 'fraction_numerical': 0.014138995072076636, 'mean_word_length': 4.253810791192838, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Soliton pair creation in classical wave scattering

We study classical production of soliton-antisoliton pairs from colliding wave packets in (1+1)-dimensional scalar field model. Wave packets represent multiparticle states in quantum theory; we characterize them by energy E and particle number N. Sampling stochastically over the forms of wave packets, we find the entire region in (E,N) plane which corresponds to classical creation of soliton pairs. Particle number is parametrically large within this region meaning that the probability of soliton-antisoliton pair production in few-particle collisions is exponentially suppressed.

Recently [16,17] a question of kink-antikink pair production in classical wavepacket scattering was addressed 1 , cf. Refs. [19][20][21][22]. The interest to this question stems from the fact that, within the semiclassical approach, wave packets describe multiparticle states in quantum theory. Studying kink-antikink creation from wave packets one learns a lot about the quantum counterpart process: production of nonperturbative kink states in multiparticle collisions. The prospect of Refs. [16,17] was to describe a class of multiparticle states leading to classical formation of kinks.
Due to essential nonlinearity of classical field equations the process of kinkantikink creation cannot be described analytically 2 and one has to rely on numerical methods. A difficulty, however, is related to the space of initial Cauchy data which is infinite-dimensional in field theory. Because of this difficulty the analysis of Refs. [16,17] was limited to a few-parametric families of initial data.
In this paper we explore the entire space of classical solutions describing solitonantisoliton pair creation from wave packets. To this end we sample stochastically over the sets of Cauchy data and obtain large ensemble of solutions, cf. Refs. [23,24]. We select solutions evolving between free wave packets and soliton-antisoliton pair and compute the energies E and particle numbers N of the respective initial states. In this way we obtain the region in (E, N) plane corresponding to classical creation of solitons. We are particularly interested in solutions from this region with the smallest N.
The model we consider is somewhat different from the standard φ 4 theory used in Refs. [16,17]. We do study evolution of a scalar field in (1 + 1) dimensions but choose nonstandard potential V (φ) shown in Fig. 1a, solid line. The reason for the unusual choice is chaos in kink-antikink scattering in φ 4 theory: it would be a venture to try applying new method in a potentially nontrivial chaotic model. We will comment on generalizations of our technique in the Discussion section.
With the above set of classical solutions we test the method of Ref. [25] where classically forbidden production of kinklike solitons in the same model was studied. Namely, we compare the boundary of the "classically allowed" region in (E, N) plane with the same boundary obtained in Ref. [25] from the classically forbidden side. Coincidence of the two results justifies both calculations.
The paper is organized as follows. We introduce the model in Sec. 2 and explain the stochastic sampling technique in Sec. 3. In Sec. 4 we present numerical results which confirm, in particular, results of Ref. [25]. We conclude and generalize in Sec. 5.

The model
The action of the model is where φ(t, x) is the scalar field; semiclassical parameter g does not enter the classical field equation We assume that the potential V (φ) has a pair of degenerate minima φ − and φ + . Then there exists a static solution of Eq.
We consider classical evolutions of φ(t, x) between free wave packets in the vacuum φ − and configurations containing soliton-antisoliton pair. Initial and final states of the process are shown schematically in Fig. 2.
We restrict attention to Pin:  In what follows we solve Eq. (2.2) numerically. To this end we introduce a uniform spatial lattice {x i }, i = −N x , . . . N x of extent −L x ≤ x i ≤ L x . At lattice edges x = ±L x we impose energy-conserving Neumann boundary conditions ∂ x φ = 0. We also introduce a uniform time step ∆t. Typically, L x = 15, N x = 400, ∆t = 0.03. Discretization of Eq. (2.2) is standard second-order 3 . We take advantage of the reflection symmetry x → −x and use only one half of the spatial lattice.
We consider the potential where dimensionless units are introduced; W (x) = e −x 2 (x + x 3 + x 5 ), a = 0.4. The value of v is chosen to equate the energy densities of the vacua, v ≈ 0.75. The potential (2.3) is depicted in Fig. 1a, solid line. We denote the masses of linear excitations in the vacua φ − and φ + by m − and m + , respectively. Our choice of the potential is motivated in two ways. First, we have already mentioned that kink dynamics in the standard φ 4 theory is chaotic [10][11][12]. The source of chaos hides in the spectrum of linear perturbations around the φ 4 kink. The latter contains two localized modes: zero mode due to spatial translations and first excited mode representing kink periodic pulsations. Localized modes accumulate energy during kink evolution which is thus described by two collective coordinates. Mechanical model for these coordinates is chaotic [12], just like the majority of twodimensional mechanical models.
We get rid of the chaos by choosing the potential (2.3) where the spectrum of linear perturbations around the soliton contains only one localized mode. Due to this property soliton motion is described by one-dimensional mechanical system which cannot be chaotic.
Let us compute the spectrum of the soliton in the model (2.3). Consider small perturbations φ − φ S (x) = δφ(x) · e ±iωt in the background of the soliton. Equation (2.2) implies, where U(x) = V ′′ (φ S (x)) and nonlinear terms in δφ are neglected. Discretization turns the differential operator in Eq. (2.4) into a symmetric (2N x + 1) × (2N x + 1) matrix; we compute the eigenvalues {ω Another, unrelated to the soliton spectrum, mechanism of chaos was proposed recently in Ref. [26]. This mechanism works under condition m + < m − which is not met in our model.
The second reason for the choice (2.3) is linearization of classical solutions at large negative times. Interaction terms should be negligible in the initial part of the classical evolution; otherwise initial wave packets cannot be associated with the perturbative Fock states. However, (1+1)-dimensional solutions linearize slowly due to wave dispersion. Brute force linearization would require large lattice which is a challenge for the numerical method. Our model is specifically designed to overcome this difficulty. At a ≪ 1 the potential (2.3) is quadratic everywhere except for the small region φ ≈ 1. Wave packets move freely in this potential if their tops are away from φ = 1, see Fig. 2. After collision the wave packets add up coherently and hit the interaction region φ ≈ 1. Below we find that all classical solutions of interest behave in the described way. We use a = 0.4 which is small enough to provide, for the chosen lattice size, linearization at the level of 1%.
It is worth noting that the problem with slow linearization is absent in multidimensional theories because amplitudes of spherical waves in D > 2 decay at power laws with distance. In general (1+1)-dimensional model linearization can be achieved artificially by switching off the interaction terms of the potential at |x| > L int . This corresponds to a physical setup where interaction takes place in a sample of length 2L int .
We compute the energy E and particle number N of the initial wave packets in the following way. Since the wave packets move freely in the vacuum φ − , where we took into account the reflection symmetry and introduced the amplitudes a k ; ω 2 k = k 2 + m 2 − . Given the representation (2.5), one calculates E and N by the standard formulas, Expression for N can be thought of as a sum of mode occupation numbers n k = |a k | 2 , where the latter are defined as ratios of mode energies ω k |a k | 2 and energy quanta 4 ω k . Note that the energy E is conserved; it can be calculated at arbitrary moment of classical evolution as In the case of free evolution this expression coincides with the first of Eqs. (2.6). Needless to say that Eqs. , ω k } of linear excitations above the vacuum φ − . In this way one obtains lattice analogs of the standing waves cos(kx) and frequencies ω k = k 2 + m 2 − . Arbitrary linear evolution in the vacuum φ − has the form where we used the eigenmode basis with normalization One extracts the amplitudes a k from the classical solution φ(t, x) at large negative t by decomposing φ(t, x), ∂ t φ(t, x) in the basis of δφ k (x) and comparing the coefficients of decomposition with Eq. (2.8). Summing up the energies and occupation numbers of different modes, one obtains where Eq. (2.9) is taken into account.

The method 3.1 Modification of the potential
It is difficult to select solutions containing soliton-antisoliton pairs in the infinite future. On the one hand, numerical methods do not allow us to extend φ(t, x) all the way to t → +∞. On the other hand, soliton and antisoliton attract; taken at rest, they accelerate towards each other and annihilate classically into a collection of waves. Thus, we never can be sure that φ(t, x) contains solitons at t → +∞, even if lumps similar to soliton-antisoliton pairs are present at finite times. We solve this difficulty by changing the value of v in Eq. (2.3) and thus adding small negative energy density (−δρ) to the vacuum φ + , see Fig. 1a, dashed line. This turns solitonantisoliton pair into a bubble of true vacuum φ + inside the false vacuum φ − [27][28][29].
Large bubbles expand at δρ > 0 since attraction between the solitons in this case is surmounted by the constant pressure δρ inside the bubble. Thus, at δρ > 0 we simply look whether solution φ(t, x) contains large bubbles at finite t. In the end of calculation, however, we have to consider the limit δρ → 0. We remark that solutions containing soliton-antisoliton pairs at t → +∞ can be identified by other methods. Our way, besides being particularly simple, has the following advantage: at δρ > 0 there exists a critical bubble [27] -unstable static solution φ cb (x) lying on top of the potential barrier between the true and false vacua. Given the critical bubble, one easily constructs classical evolutions between the vacua. Indeed, in the critical bubble attraction between the soliton and antisoliton is equal to repulsion due to δρ. Being perturbed, it either starts expanding or collapses forming a collection of waves in the vacuum φ − . Thus, adding small perturbation to the critical bubble and solving classical equations of motion forward and backward in time, one obtains the classical solutions of interest. Critical bubble at δρ = 0.4 is depicted in Fig. 3a. Let us obtain a particular solution describing creation of expanding bubble from wave packets at δρ > 0. We solve numerically Eq. (2.4) with U(x) = V ′′ (φ cb (x)) and find the spectrum of linear perturbations {δφ k } around the critical bubble. This spectrum is shown in the inset in Fig. 3a; it contains precisely one negative mode δφ neg (x), ω 2 neg < 0 due to changes in the bubble size. The latter mode describes decay of the critical bubble, where we fix t 0 = 0, B neg = 10 −2 in what follows. Using the configuration (3.1) and its time derivative as Cauchy data at t = 0, one solves numerically Eq. (2.2) forward and backward in time and obtains φ(t, x), see 6 Fig. 3b. The latter solution interpolates between free wave packets above the vacuum φ − and expanding bubble. We compute the values of (g 2 E, g 2 N) ≈ (6.  In Fig. 3b we check the linearity of evolution at t → −∞ by comparing the linear and exact energies of φ(t, x), Eqs. (2.10) and (2.7). As expected, the linear energy coincides with the exact one at large negative times and departs from it when wave packets collide. In what follows we estimate the precision of linearization as fractions of percent.

Stochastic sampling technique
In this Section we sample over classical solutions with bubbles of true vacuum in the final state, see Refs. [23,24]. It is hard to pick up the initial Cauchy data for such solutions: most of the initial wave packets scatter trivially and do not produce expanding bubbles. Instead, we consider the data {φ(0, x), ∂ t φ(0, x)} at t = 0. We decompose these data in the basis of perturbations around the critical bubble, where the negative mode is treated separately. Note that any functions φ(0, x), ∂ t φ(0, x) can be written in the form (3.2). For each set of Cauchy data {A neg , A k , B neg , B k } we solve numerically Eq. (2.2) and obtain classical solution φ(t, x). Due to instability of the critical bubble there is a good chance to obtain transition between the vacua φ − and φ + .
Next, we study the region in (E, N) plane corresponding to classical formation of expanding bubbles from colliding wave packets. We are particularly interested in the lower boundary N = N min (E) of this region. Let us organize the artificial ensemble of solutions describing transitions between the vacua. The probability of finding each solution in our ensemble is proportional to where E and N are the energy and initial particle number of the solution; τ and ϑ are fixed numbers. At large positive ϑ solutions with the smallest N dominate in the ensemble and we obtain the boundary N min (E) with good precision. Value of τ controls the region of energies to be covered. We use Metropolis Monte Carlo algorithm to construct the ensemble (3.3). In our approach solutions are characterized by the coefficients in Eq. (3.2); condition A neg = 0 is used to fix the time translation invariance of Eq. (2.2). The algorithm starts from the solution (3.1) describing decay of the critical bubble; it has B neg = 10 −2 , A k = B k = 0. Denote the energy and particle number of this solution by (E 0 , N 0 ). We pick up a random coefficient from the set {A k , B k , B neg } and change it by a small step. The latter step is a Gauss-distributed random number with zero average and small dispersion 7 σ. Substituting the modified set of coefficients into Eqs. (3.2), we find φ(0, x) and ∂ t φ(0, x). Then, solving numerically the classical field equation, we obtain the entire solution φ(t, x). We compute the values of (E, N) by Eqs. (2.10). Solution is rejected if it does not interpolate between the vacua φ − and φ + ; otherwise we accept it with the probability p accept = min 1, e −τ ∆E−ϑ∆N , (3.4) where (∆E, ∆N) are differences between the new values of (E, N) and the values (E 0 , N 0 ) for the solution we started with. If the new solution is accepted, we write it down and use its parameters {A k , B k , B neg }, (E, N) → (E 0 , N 0 ) for the next cycle of iterations. After many cycles we obtain the ensemble (3.3) of accepted solutions. A typical run of the Monte Carlo algorithm is shown in Fig. 4a where the accepted solutions are marked by dots in (E, N) plane. The algorithm starts in the vicinity of the critical bubble, then moves to smaller N and finally arrives to the boundary N min (E) where the majority of solutions is found.

Numerical results
We perform Monte Carlo runs at different values of τ and ϑ until the entire curve N = N min (E) is covered with solutions. In total we obtained 2 · 10 7 solutions, where the value of τ was ranging between 0 and 10 4 ; ϑ = 10 3 , 10 4 , 5 · 10 4 . The boundary N = N min (E) is constructed by breaking the energy range into small intervals ∆E = 0.01 and choosing solution with minimal N inside each interval. This is the result we are looking for: N min (E) gives the minimum number of particles needed for classically allowed production of bubbles. It is plotted in Fig. 4b, solid line. As expected, N min (E) starts from (E, N) = (E cb , N cb ) and decreases monotonously with energy. At high energies N min (E) is approximately constant. Note that the particle number is parametrically large in the "classically allowed" region, N min ∼ 1/g 2 . This means, in particular, that the probability of producing the bubble from few-particle initial states is exponentially suppressed.
Given the boundary N min (E), we check results of Ref. [25] where classically forbidden transitions between N-particle states and states containing the bubble were considered. The probability of these processes is exponentially suppressed in the semiclassical parameter, where F N (E) is suppression exponent. One expects that this exponent vanishes in the "classical" region N > N min (E). We extract the boundary of the set F N (E) = 0 from the results of Ref. [25] and plot this boundary in Fig. 4b (dashed line). It coincides with N min (E) within 0.5% accuracy; the agreement justifies both calculations. Let us look at solutions with almost-minimal initial particle number, N ≈ N min (E). Two such solutions are plotted 8  moving in the vacuum φ − . After collision the wave packets emit waves and form the bubble.
The most surprising part of the evolutions in Fig. 5 is emission of waves during the bubble formation. One assumes that the role of these waves is simply to carry away the energy excess which is not required for the creation of bubble. Indeed, solutions at different E look alike, cf. Figs. 5a and 5b; besides, N min (E) is independent of energy at high values of the latter. In Refs. [30][31][32] it was assumed that there exists certain limiting energy E l which is best for bubble creation. Then, classical solutions with minimal particle number at E > E l are sums of two parts: non-trivial soft part describing bubble production at E = E l and trivial hard part -waves propagating adiabatically in the soft background. Hard waves carry away the energy excess E − E l without changing the initial particle number; this is achieved at small wave amplitudes and high frequencies, see Eqs. (2.10).
Numerical results do not permit us to judge whether the limiting energy exists. We can, however, confirm the conjectured structure of high-energy solutions. Consider the energies ǫ k = ω k |a k | 2 /g 2 of the modes at t → −∞, where the amplitudes a k and frequencies ω k are defined in Eq. (2.8). In Fig. 6 we plot these energies for the two solutions depicted in  tail of excited high-frequency modes is seen, however, in the graph representing the high-energy solution 5b. The tail carries substantial energy while the corresponding modes propagate adiabatically and do not participate in nonlinear dynamics. This is precisely the behavior conjectured in Refs. [30][31][32].
It is reasonable to assume that solutions at arbitrarily high energies have the same "hard+soft" structure. Then, N min (E) stays constant as E → +∞ and classical formation of bubbles is not possible at any energies unless the initial particle number is larger than N min (E = +∞). Finally, we consider the limit δρ → 0. A typical solution at small δρ is depicted in Fig. 7a. It describes creation of soliton and antisoliton which move away from  Fig. 7b (dashed lines). They are almost indistinguishable; thus, the limit δρ → 0 exists. Extrapolating N min (E) to δρ = 0 with linear functions, we obtain the region in (E, N) plane for classically allowed production of soliton-antisoliton pairs (above the solid line in Fig. 7b). All initial wave packets leading to classical creation of soliton pairs have the values of (E, N) within this region. The region in Fig. 7b is qualitatively similar to the regions at δρ > 0; in particular, N min (E) is constant at high energies.

Discussion
In this paper we studied multiparticle states leading to classically allowed production of soliton-antisoliton pairs in (1+1)-dimensional scalar field model. We characterized these states with two parameters -energy E and particle number N; we have found the corresponding "classically allowed" region in (E, N) plane. There were two main ingredients in our technique. First, we added constant pressure δρ pulling soliton and antisoliton apart. This modification led to appearance of the critical M M F att −g m H g m H H bubble -unstable static solution lying on the boundary between the perturbative states and soliton-antisoliton pair. Second, we applied stochastic sampling over Cauchy data in the background of the critical bubble. We thus obtained large ensemble of classical solutions describing formation of soliton-antisoliton pairs. Calculating the values of (E, N) for each solution we found the required "classically allowed" region.
Our method is naturally generalized to higher-dimensional models. For example, consider formation of t'Hooft-Polyakov monopoleantimonopole pairs in four-dimensional gauge theories [33,34]. Constant force dragging monopole and antimonopole apart is provided [35] by external magnetic field H, see the figure. At H = 0 there exists a direct analog of the critical bubble [36]: unstable static solution where the attractive force F att between the monopole and antimonopole is compensated by the external forces ±g m H (g m is a magnetic charge of the monopole). One performs Monte Carlo simulation in the background of this static solution and obtains many classical evolutions between free wave packets and monopoleantimonopole pairs.
A particularly interesting application of our technique might be found in the study of kink-antikink production in (1+1)-dimensional φ 4 theory. In this model the boundary N min (E) is lowered [16,17] due to chaos. We do not expect any difficulties related to nontrivial dynamics of solitons. However, classification of initial data for kink-antikink formation may require modification of our method.