Numerical study of multiparticle scattering in λϕ4 theory

We study numerically classical collisions of waves in λϕ4 theory. These processes correspond to multiparticle scattering in the semiclassical regime. Parameterizing initial and final wavepackets by energy E and particle numbers Ni, Nf we find classically allowed region in the parameter space. We describe properties of the scattering solutions at the boundary of the classically allowed region. We comment on the implications of our results for multiparticle production in the quantum regime.


Introduction
Multiparticle production in weakly coupled bosonic field theories has long been of great interest (see e.g. [1,2] for reviews), which was initially related to instanton-induced processes [3,4]. Tree level calculations performed in (3 + 1)-dimensional scalar λφ 4 theory show [5][6][7][8][9] that multiparticle amplitudes of 1 → N f processes near the threshold grow factorially with the number of produced particles N f . This behaviour persists at one-loop level [10] and violates unitarity at sufficiently large multiplicities N f ∼ 1/λ. It was observed [11], that corresponding multiparticle cross section of 1 → N f processes can be conveniently written in the following functional form in the limit λ → 0 with λN f and E fixed, where E is the average kinetic energy of the final particles. Inspired by this exponential behaviour several semiclassical techniques for calculation of the multiparticle cross sections were developed in [12][13][14][15][16]. They involved singular solutions of classical equations of motion and allowed to extend the observations obtained with perturbative calculations to higher energies. Still in spite of those vast efforts the most reliable results were obtained in the regime of small λN f where corresponding cross sections turn out to be exponentially suppressed. Unitarity bounds [2,17,18] indicate that such suppression of the probability (1.1) takes place at all particle numbers. Interesting insights come from studies of an analogue of the scattering amplitudes in (0+1)-dimensional theory, i.e. for the anharmonic oscillator with a quartic potential [19][20][21], where it was found that the perturbative factorial growth is replaced by the exponential suppression for λN f > ∼ 1. Let us note that recently an interest to the problem of multiparticle production has been renewed [22][23][24][25][26][27][28][29][30][31][32]. In particular, calculations of refs. [28,32] performed in λφ 4 theory JHEP11(2018)068 with spontaneously broken Z 2 symmetry, exhibited that, contrary to the common belief, the probability of 2 → N f scattering process near the threshold is not suppressed at λN f 1 but exhibits an exponential growth at least near the threshold. This conclusion was used to invite the Higgsplosion scenario for solution of the hierarchy and fine-tuning problems of the Standard Model. Very recently authors of refs. [33,34] argued that the Higgsplosion scenario and in particular the exponential growth of the multiparticle production probability is not consistent with basic principles of quantum field theory and at least requires some modification. Therefore behaviour of the probability (1.1) in λφ 4 field theory (in broken and unbroken phases) at large λN f is still an open problem even in the weak coupling regime which will be assumed in the paper.
It was long ago understood that the processes f ew → N f for large N f can be studied starting from the processes N i → N f where both initial N i and final N f particle numbers are large, i.e. of order 1/λ. In this case, one can consider classical counterpart of the quantum scattering processes, i.e. collisions of classical wave packets. If the ingoing and outgoing waves are in the linear regime at t → ±∞ corresponding initial and final field configurations can be associated with coherent states having average energy E, initial N i and final N f particle numbers. The probability of corresponding quantum scattering in the semiclassical regime is not exponentially suppressed and the whole family of such solutions span a classically allowed region in the space of parameters E, N i and N f . At a given energy E * and final particle number N f one can try to minimize the initial particle number N i with respect to initial conditions. If the minimum goes to zero for some energy then the probability of the processes 2 → N f is not exponentially suppressed at E > E * [18]. Otherwise the existence of a nontrivial minimum at N i ∼ 1/λ would indicate on an exponential suppression of the 2 → N f scattering probability for E < ∼ E * . This idea of exploring the classical region was used previously for studies of several different processes induced by collisions of particles. Among of them are the false vacuum decay [35], baryon number violating processes in the Standard Model [36,37] and soliton-antisoliton pair production in the (1 + 1)-dimensional scalar field theory [38].
In this paper we study classical solutions describing scattering of the wave packets in the unbroken λφ 4 theory. In particular, we try to approach the boundary of the corresponding classically allowed region. The motivation for this study is two-fold. First, by identifying classically forbidden and allowed regions we can unambiguously tell which multiparticle scattering processes are suppressed in the semiclassical limit and which are not. In particular, our results give clear indication for the suppression of the probability of f ew → N f processes even in the regime λN f > ∼ 1, λ 1 which is in agreement with unitarity arguments. Note that it would be very interesting to apply the method of classical solutions used in our study to explore the classically allowed region in the (E, N i , N f ) parameter space in the spontaneously broken λφ 4 theory in view of direct relevance to the Higgs physics. Another part of our motivation is related to the boundary of the classically allowed parameter space and corresponding classical solutions. Using suitable semiclassical methods one can start obtaining solutions of the classical equations of motion which describe multiparticle scattering processes N i → N f in the classically forbidden region in the semiclassical approximation and calculate the semiclassical exponent similar to those JHEP11(2018)068 in eq. (1.1). As one approaches the boundary of the classically allowed region in the (E, N i , N f ) parameter space the semiclassical exponent is expected to nullify and comparison with the boundary and corresponding classical solutions obtained in the present study would be a valuable check of the semiclassical procedure which is known to plague from multiple branches of solutions.
Classical scattering of waves in relation to multiparticle production was previously studied to some extent in the (1+3)-dimensional φ 4 model [39], (1+1)-dimensional abelian Higgs model [40] and non-abelian gauge theories [41][42][43]. In particular, authors of ref. [39] constructed an initial wave packet consisting of a few high frequency free field modes and examine its evolution. They found no significant energy transfer to low frequency modes which would indicate to the production of many "quanta". In our study we use stochastic sampling technique to scan numerically over classical solutions describing the multiparticle scattering. We limit ourselves to spherically symmetric solutions which reduces the system to a (1 + 1)-dimensional model and makes the problem numerically feasible. We find it is technically more convenient to fix not final N f but initial particle number N i and energy E and find the minimal and maximal values of the particle number N f of outgoing wave packets which is obtained by solving classical equations of motion. Both approaches are equivalent due to the time reversal symmetry. With our numerical methods we find a nontrivial classically allowed region for the multiparticle scattering processes. Namely, we are able to determine the upper and lower boundaries N min f and N max f as functions of energy E for a set of fixed values of N i . We observe that the change in particle number in the scattering reaches values up to 22% for the studied energy range. We examine properties of the classical solutions corresponding to the boundaries.
The rest of the paper is organized as follows. Section 2 is devoted to the description of the model and introduction of notations and most relevant quantities. In section 3 we describe numerical method which we utilize to study solutions to the classical equations of motion describing scattering of wave packets with a particular interest to those "providing" maximal change in their particle numbers. In section 4 we present our numerical results. Section 5 is reserved for a discussion and conclusions.

The model
We consider (3 + 1)-dimensional model of a scalar field with the action The parameter m 2 is taken to be positive which corresponds to the unbroken phase of this theory. By making rescaling x µ → m −1 x µ and φ → m 2 λ φ this action can be cast into In this form λ is a semiclassical parameter which does not enter equations of motion. We are interested in a particular type of classical solutions, which describe collisions of wave JHEP11(2018)068 packets. In what follows we limit ourselves to spherically symmetric field configurations. In this case it is useful to make the following redefinition and we obtain the action of (1 + 1)-dimensional theory on a half-line with the boundary condition χ(t, 0) = 0 and the position dependent interaction. The corresponding equation of motion reads and the energy related to the field configuration χ(t, r) is In what follows we solve the classical e.o.m. (2.5) numerically. To do this we restrict our solutions to a finite space interval [0, R] where R is sufficiently large. At r = R we impose Neumann boundary condition ∂ r χ = 0. This restriction allows us to expand the field configuration as follows c n (t) 2 R sin k n r, where k n = π 2R (2n + 1), n = 0, 1, . . .

(2.7)
Inserting this expansion in the action (2.4) we obtain the following equations of motion c n + ω 2 n c n + I n = 0 , n = 0, 1, . . . (2.8) where ω 2 n = k 2 n + 1 and the interaction term reads dr χ 3 (t, r) r 2 sin k n r , n = 0, 1, . . . Below we consider solutions which linearize at initial and final times. In this case they can be interpreted as describing a multiparticle scattering. In the linear regime the timedependent Fourier components c n (t) can be conveniently written via positive and negative frequency components as follows 2ω n a n e −iωnt + a * n e iωnt as t → −∞ (2.11)

JHEP11(2018)068
and c n (t) → 1 √ 2ω n b n e −iωnt + b * n e iωnt as t → +∞. (2.12) Using these representations one can compute the energy of colliding wave packets as Here |a n | 2 and |b n | 2 can be thought as occupation numbers for initial and final states and their sums N i = 4π λ n |a n | 2 , N f = 4π λ n |b n | 2 (2.14) are initial and final particle numbers, respectively. For convenience we introduce the following shorthand notations For the numerical implementation we truncate the expansion (2.7) at n = N r and solve the evolution equations (2.8) using Bulirsch-Stoer method, see e.g. [44]. Some of the results have been verified with the 4-th order Runge-Kutta method with a very small time step. For convenience, we introduce the uniform spacial lattice r i , i = 0, . . . , N r with r 0 = 0 and r Nr = R. In what follows we use several values of R = 20, 30 and 50 and N r = 400, 600 and 1000 to study the dependence of our results on the lattice. The time interval is taken to be somewhat smaller than ∼ 2R. We use FFTW implementation of discrete Fourier transformation [45] to compute the field configuration χ(t, r) as (2.7) and the interaction term (2.9).

The method
In this section we describe a method used to explore the classical transitionsÑ i →Ñ f at fixed energyẼ. In particular, we are interested in solutions of the classical e.o.m. which maximize |Ñ f −Ñ i | at fixed values ofÑ i andẼ.

Initial conditions
We take initial conditions which correspond to a wave packet which is localized well away from the interaction region and which is propagating towards r = 0. Technically, we choose an interval [r 1 , r 2 ], where r 1 , r 2 are chosen points belonging to the spacial lattice, i.e. r 1 = r i 1 and r 2 = r i 2 for some i i < i 2 , and construct following function Herek n = πn r 2 −r 1 ,ω n = k2 n + 1, n = 1, . . . , i 2 − i 1 and f n are initial Fourier amplitudes. By construction the configuration (3.1) nullifies at the boundaries r = r i 1 and r = r i 2 .

JHEP11(2018)068
Next, to make the initial configuration smooth enough we introduce an artificial smearing multiplying (3.1) by the function where d and a are chosen parameters. We checked that our results have a very mild dependence on the precise way of the smearing. In what follows we fix d = 0.2 and a = 0.1 for concreteness. Next, we obtain Fourier amplitudesf n of the smeared configuration by inverting (3.1) and take the initial condition for the e.o.m. (2.8) (i.e. the field configuration and its first time derivative at initial time t = t i ) from This corresponds to the initial wave packet propagating towards r = 0. In what follows we set t i = 0. We consider only solutions whose evolution is linear at initial (and final) times.

Going to the classical boundary
Using the initial field configuration we calculate its energyẼ and initial particle number N i by using eqs. (2.10), (2.11), (2.13) and (2.14). Our aim is to find the classical solutions describing scattering of wave packets in which the particle number changes as much as possible at a given energyẼ. Since we solve initial value problem it is more convenient to fix the initial particle numberÑ i and find the absolute minimum (or maximum) ofÑ f with respect to chosen set of initial conditions, i.e. f n entering (3.1). We fix the initial particle number by the respectful normalization of the Fourier amplitudesf n . The final particle numberÑ f is a highly non-linear function of initial dataf n and it may have several local extrema. To find its absolute minimum (maximum) we use the stochastic sampling technique in combination with the simulated annealing method, see [44,46]. Namely, we generate an ensemble of the classical solutions with a fixed initial particle numberÑ i weighted by the probability where If large positive β and ξ are taken the ensemble will be dominated by solutions having small F thus driving their distribution towards the lower boundaryÑ min f (Ẽ * ) of the classically allowed region. To reach the upper boundaryÑ max f (Ẽ * ) signs of β and ξ should be negative. To generate the ensemble (3.4) we use the Metropolis Monte Carlo algorithm. We start with a solution specifying by a randomly chosen set of f n , see (3.1), which have a fixed initial particle numberÑ i . This solution has an energyẼ and final particle numberÑ f which is found after solving the classical equations of motion. Next, we choose a few randomly chosen amplitude numbers and change corresponding amplitudes f n by small amounts,

JHEP11(2018)068
f n → f n = f n + ∆f n . We used to take three amplitudes at a time which shows a good performance of our algorithm. The quantities ∆f n are chosen to be normally distributed with the standard deviation a/ √ω n , where a is a small number. In our calculations its value is varied from 10 −4 to 0.1. The modified set of amplitudes f n is rescaled in a proper way and used to construct an initial wave packet having the same initial particle number N i and some value of energyẼ , see section 3.1. Then we evolve the system forward in time far enough until it reaches the linear regime and we can find the final particle number N f of the new solution. We compute ∆F ≡ F − F using (3.5). The new set of amplitudes f n is accepted with probability p accept = min 1, e −∆F (3.6) and used for the next iteration. Typically, we fixÑ i andẼ * and perform several runs starting with different β and ξ. Each run spans of order 10 3 − 10 4 iterations. The value of ξ remains a constant during a single run, while β which is analogue of inverse temperature is gradually increased from its initial value β 0 according [47] to where i is the iteration number. Values of β 0 and ξ are temporarily increased from 10 to 10 6 and from 10 −4 to about 1, respectively. We stop the procedure when relative change in the value of F during a run becomes smaller than 10 −3 for the same ξ and it does not increases with the variation of a and β 0 . Let us note, that the stochastic sampling technique was previously used to study classically allowed sphaleron transitions [36,37] and classical soliton-antisoliton pair production in particle collisions in (1+1)-dimensional scalar field theory [38], where its efficiency in determination of the boundary of the classically allowed region was confirmed. the lineÑ f = 1 on the (Ẽ,Ñ f ) plane. Initially, we take it relatively large value of a = 0.1 and start our stochastic sampling with small β 0 about 10.0. During single run its value is changed according to (3.7). We repeat the runs with a larger starting value of β. The value of a which determines size of the Fourier amplitude changes ∆f n is also decreased after several runs. At the same time the value of ξ is increased to fix the energy of the solution to be nearẼ * . The obtained domain of the classical solutions in (Ẽ,Ñ f ) plane, see figure 1, has smooth envelopes,Ñ min f (E) andÑ max f (E), which represent the boundary of the classically allowed region of the lattice version of the continuous model. We see that the maximal possible change in the particle number |Ñ f −Ñ i | during the classical evolution increases with the energy of solutions but does not exceed 12% in the chosen energy range. Let us stress that due to the time reversal symmetry one can always interchange the initial N i and finalÑ f particle numbers in the discussion of the classically allowed region in the figure 1. For all boundary solutions which we have found withÑ i = 1 the initial wave packet has a relatively sharp part near its end. This part becomes more pronounced with increase of its energy as seen in figure 3. The initial wave packet propagates towards the interaction region and its spiky parts increases considerably near the point r = 0. The reflected wave packet have the form very similar to the incoming one. The spiky part of the incoming wave packet is the last to arrive into the interaction region and the first to leave it.

Numerical results
According to the used extremization procedure which involves the stochastic sampling, the found classical solutions describe the scattering of waves withÑ i = 1 and a maximum value of |Ñ f −Ñ i | at a given energyẼ. It is interesting to study how the particle number changes during the time evolution. In section 2 we introduced the particle numbers for the initial and final wave packets. However, we can calculate an instantaneous particle number N (t) by using eq. (2.11) (or (2.12)) as a definition of time-dependent positive and negative frequency components at arbitrary time t and applying formulas similar to (2.14). The quantityÑ (t) coincides withÑ i andÑ f at initial and final times where the evolution of the field is linear. On the right panel of the figure 4 we plot the instantaneous particle number for the boundary solution shown in figure 3. Comparing it with the time evolution of the field presented on figure 3 we see that the actual change of the particle number occurs precisely when the sharpest part of the initial wave packet reaches the interaction region. For instance, the deviations ofÑ at time stamps t = 12.4 and t = 21.3 (see figure 3) from its asymptotic valuesÑ i andÑ f are less than 0.0005, which are much smaller thenÑ f −Ñ i ≈ 0.072. Such behaviour ofÑ (t) give us support that found extrema of F are independent of choice of the interval [r 1 , r 2 ] which is used to generate initial field configurations as far as it is sufficiently large and include the part of the wave packet responsible for the most of the change inÑ (t). Also we verified that numerical error in |Ñ f −Ñ i | related to the smearing is less than 10 −3 . In a way similar toÑ (t) we define a linearized energyẼ lin (t) which is used to check the linearity of the solution by comparing it with the exact energy (2.10). Such comparison is presented in figure 4 (left panel) for the solution from figure 3. The total energy is conserved on the entire solution with the accuracy less than 10 −3 , while it coincides with the linearized energy better than 10 −4 at initial and final times. A nontrivial change in particle number implies a redistribution of the energy between low and high frequency modes. In figure 5   where ∆k = π R . We observe expected softening of the energy distributions for the final wave packet as compared to the initial one. However, the effect is quite small even for boundary solutions which describe processes with maximal change of the particle number. Also we see that with increase of the collision energyẼ the modes with larger k n (i.e. higher frequencies) become filled in and atẼ > ∼ 11.0 the tail of energy distributions reaches the largest value of k n for the chosen lattice size and spacing. This is a clear indication that a finer lattice is required to obtain reliable results at these and higher energies. At the same time solutions with energies not far from the threshold energyẼ th =Ñ i consist mostly of nonrelativistic modes and a larger space interval R is needed to reach linear regime at initial and final times. For this reason we do not consider solutions with energies lower than 1.5.
Next, we study the dependence of the obtained results for the boundary of the classically allowed region in the (Ẽ,Ñ f ) plane forÑ i = 1 and for the boundary solutions on parameters of the lattice, namely on the spacial cutoff R and on the number of modes N r . For comparison to the case with R = 20 and N r = 400 we repeat the same numerical procedure to find the boundary of the classically allowed region for the case R = 20, N r = 600 with a smaller lattice spacing and for the case R = 30, N r = 600 with a larger space interval but the same lattice spacing. The results are presented in figure 6, where we show lower and upper boundariesÑ min f (E) andÑ max f (E) for the different cases. The space intervals for selection of the initial wave packets are taken to be [6.7, 19.2] in the former case and [8.5, 29.8] in the latter. We observe that different boundaries coincide with accuracy better than 4 · 10 −3 . A deviation is seen for the case N r = 600, R = 30 for large energiesẼ >  difference lies in the enlarged soft oscillating part of the solution. We found that the time evolution of the particle numberÑ (t) for these solution is similar to that of presented in figure 4. Namely, the actual change inÑ starts with arrival of the sharp tail of the incoming wave train into the interaction region. Energy distributions over Fourier modes for the solutions with the same energy obtained on different lattices are almost indistinguishable as seen from figure 8.  very similar properties. Let us note that the time reversal symmetry translates a lower boundary solution with initialÑ i and finalÑ f particle numbers into an upper boundary solution with the initial and final particle numbers interchanged.
Now we turn to discussion of numerical results for another values of the initial particle numberÑ i . In figure 11 we show the upperÑ f (E) and lowerÑ f (E) boundaries of the classically allowed region forÑ i = 0.1 and the relevant energy interval 1.5Ñ i , 15.0Ñ i with the values ofẼ * defined by (4.1). We see that the maximum relative difference of the initial and final particle numbers, i.e.
, is more than two order of magnitude smaller than for the caseÑ i = 1 for energiesẼ < ∼ 15Ñ i . Numerically, the particle number changes no more than 0.33% for the chosen energy range. In figure 12 we show examples of the initial wave packets and instantaneous particle number evolution for the upper and lower boundary solutions forẼ ≈ 0.6. The space interval for the initial wavepacket [r 1 , r 2 ] is shown by the dashed lines. Qualitative properties of these solutions are very similar to those obtained forÑ i = 1. Namely, the solutions have a spiky part whose transition through the interaction region results in actual change in particle number and a soft oscillating part. As in the previous case we obtain several branches of the boundary solutions which differ by corresponding shifts in time.
Let us now turn to the caseÑ i = 10 in which we obtain numerical results concerning boundary solutions which are qualitatively different from those obtained at smaller values of N i = 0.1, 1.0. In figure 13 boundary of the classically allowed region is shown with the thick solid (red) and thick dashed (blue) lines for the energy interval 1.5Ñ i , 15.0Ñ i with the same set (4.1) ofẼ * . We find that although the form of the boundary is similar to that of for   Figure 14. From top to bottom: a) initial wave packet; b) evolution of the instantaneous particle numberÑ (t); c) distribution of energy over the modes k per wave number for the initial field configuration; d) the same as in c) but for the final field configuration. Left and right panels correspond to the lower boundary solutions withẼ ≈ 65 and 80, respectively.

JHEP11(2018)068
on left panels. The oscillating pattern in k indicates that the incoming (and outgoing) field configuration looks as a sum of two separated in space wavetrains each having a similar smooth Fourier image. In particular, the distance δr between the spikes in the initial configuration and the oscillation period δk in the initial energy distribution are related approximately by δr · δk ≈ 2π for all solutions of this branch. For example, the solution withẼ ≈ 65 presented on left panels of figure 14 has δr ≈ 2.8 and δk ≈ 2.3. Comparing the energy distributions k for the field configurations at initial and final times one observes a small energy transfer from low to high frequency modes. At energies larger than about 70 our numerical results show that the absolute minimum (and maximum) ofÑ f is delivered by a different branch of the solutions whose initial (and final) space field configuration contains already three spikes. On the right panels in figure 14 we present the initial field configuration φ(t = 0, r), time evolution ofÑ (t), and energy distributions over the Fourier modes k for the initial and final field configurations for a boundary solution withẼ ≈ 80 having three spikes. On this branch of solutions the initial and final field configurations looks as a sum of three separated in space wavetrains with similar Fourier image. As in the cases of smallerÑ i , the instantaneous particle numberÑ (t) undergoes the most dramatic changes when spiky parts of the field configurations reach the interaction region. We find that the two-spike solutions atẼ > ∼ 70 and the three-spike solutions atẼ < ∼ 70 still represent local minima ofÑ f . To figure that out we take a two-spike lower boundary solution atẼ ≈ 65 as a seed for our Monte-Carlo procedure with E * > 70. To ensure that the field configuration does not jump to another (three-spike) branch we take a relatively small value of a which governs the size of changes in the amplitudes f n . In this way we find that two-spike branch of solutions atẼ > ∼ 70 which is shown in figure 13 by thin red (solid) line. In similar way we find a continuation of the three-spike branch of the solutions at E < ∼ 70 which is shown by this blue (dashed) line in figure 13. Upper boundary solutions have very similar properties. 1 Let us note that in the chosen energy range the maximal changes in particle number reach values around 20%.
There have been also found another branches of solutions which deliver only local extrema toÑ f . In the linear regime these solutions look as a sum of two or three singlespike wavetrains which differ from those already described by somewhat different distances between the wavetrains. In figure 15 we show two examples of such solutions. Left panels show a solution of this type near the upper boundary with two spikes andÑ f >Ñ i . The distance between the wavetrains in the initial field configuration is larger than that of for the boundary two-spike solution. Similarly, the right panels show an example of the solution near the lower boundary delivering a local minima toÑ f with three spike. Again we observe somewhat larger distance between the wavetrains in the initial field configuration, cf. figure 14 (upper right panel). These two branches of solutions are shown in figure 13 by dots. Let us note, that forÑ i = 10 we do not find any field configuration which would deliver (even local) extrema to the functionalÑ f (E) and whose initial (final) wave packet would have a single spike only. We tried to find it on purpose by taking as a JHEP11(2018)068  Figure 15. The initial wave packets (left panels) and evolution of the instantaneous particle number (right panels) for the solutions withÑ i = 10 andẼ ≈ 80 corresponding to local minima of N f . Upper and lower panels correspond toÑ f >Ñ i andÑ f <Ñ i . seed some initial field configuration which was obtained forÑ i = 1 as a boundary solution which is of single-spike type and increasingÑ i by small steps toÑ i = 10. We find that at intermediate values ofÑ i around 5-6 the initial field configurations of the boundary configurations prefer to be splitted into more than one wavetrains to produce larger change in particle number Ñ f −Ñ i .
The picture becomes even more complicated at largerÑ i . In figure 16 we show our numerical results for the boundary of the classically allowed region forÑ i = 30 taking R = 30, N r = 600 (solid blue line). We observe that the difference |Ñ f −Ñ i | does not exceed 22% in this case in the chosen energy interval. We find that the corresponding boundary solutions contain already 4-7 spikes in their initial (and final) wavetrains. Examples of the initial field configurations of the boundary solutions are shown in figure 17. At the same time we find that the number of local extrema ofÑ f increases dramatically at largeÑ i . Corresponding solutions have similar number of wavetrains but with different distances between them. This greatly complicates the task of finding the boundary of the classically allowed region because, on the one hand, many of our numerical runs get stuck in such local minima and, on the other, the values ofÑ f at the sameẼ for different branches of near-boundary solutions are at some cases seem to be close to our numerical accuracy. For these reasons at present we are unable to classify in details the boundary solutions as it has been done for the cases with smallerÑ i . In the figure 16 we present only the envelopes f ,Ẽ (2) , exists then there should also exist a solution which has energy and particle numbers equal to the following sums Such a solution can be constructed explicitly by taking it as a sum of the individual solutions sufficiently separated in space-time. This observation, in particular, means that the width of the classically allowed region, i.e. |Ñ max f −Ñ i | or |Ñ min f −Ñ i | at fixed ratioẼ/Ñ i should grow faster than a linear function with increase of initial particle numberÑ i . On the figure 18 we show 2 the dependence of the quantity |Ñ min f −Ñ i | as a function ofÑ i for two values ofẼ/Ñ i . For comparison a linear function is shown in thin (blue) line. One can see that the energy dependence of |Ñ min f −Ñ i | tends to be linear at largeÑ i . We found that |Ñ max f −Ñ i | has similar behaviour. This can be considered as an indication on the existence of some limiting boundary of the classically allowed region in the plane (Ẽ/Ñ i ,Ñ f /Ñ i ) at largeÑ i or in the plane (Ẽ/Ñ f ,Ñ i /Ñ f ) at largeÑ f due to the symmetryÑ i ↔Ñ f .

Discussions and conclusions
In this paper we considered the classical scattering of wavepackets in the unbroken λφ 4 theory. Initial and final field configurations were characterized by energy E, initial N i and final N f particle numbers. We found that although the maximal change of particle number in the scattering grows with energy at fixed value of N i (or N f ), the value of does not exceed 22% for the considered energy range E < ∼ 10N f m. In the quantum counterpart of the problem the initial and final wavepackets correspond to initial and final coherent states with given average values of energy and particle number. Our classical JHEP11(2018)068 approach to the multiparticle scattering is valid when occupation numbers are large, i.e. E, N i , N f ∼ 1 λ and λ 1. In this study we numerically obtained the classically allowed regions of these processes in the (E, N f ) plane at several values of N i (or, equivalently, in the (E, N i ) plane at several values of N f ). They are shown in figures 1, 11, 13 and 16.
Here we considered spherically symmetric field configurations only. We believe that this simplified consideration has much in common with most general setup.
The scattering processes with the parameters E, N i and N f outside the classically allowed regions are classically forbidden. The common lore here is that their probability in quantum theory is expected to be exponentially suppressed, i.e. have the form Ae 1 λ F with some negative suppression exponent F and a prefactor A. At the same time the probability of the classically allowed processes is not exponentially suppressed. In the most interesting case 2 → N f production, with λN f > ∼ 1, the initial state contains semiclassically small number of particles. Our results show that such processes (as well as N f → 2 scattering) lie deeply in the classically forbidden region and therefore their probabilities are expected to be exponentially suppressed at least for 4π λ N f < ∼ 30 and E < ∼ 10N f . Limitations of our numerical procedure do not permit us to obtain the boundary of the classically allowed region for arbitrarily large values of N f and E. However, we obtain an evidence for existence of a limiting boundary of the classically allowed region in (E/N i , N f /N i ) plane at large N i .
In summary, our results imply that the probability of the multiparticle production processes f ew → N f in the unbroken weakly coupled λφ 4 theory is suppressed even at large λN f > ∼ 1. Therefore, factorial growth of multiparticle amplitudes observed at λN f 1 should stop at its larger values leaving the probability consistent with unitarity. Based on the perturbative results at small λN f and on the experience with other non-perturbative processes such as instanton-like transitions, false vacuum decay, soliton pair production in particle collisions we expect this suppression most likely to be exponential in the sense of eq. (1.1). Let us note that the analysis described in this study can be (at a cost of a considerably larger CPU time) extended to the case of field configurations with axial symmetry which would be more appropriate for discussion of the limit of two-particle collisions. Also, it would be interesting to apply the method of classical solutions to the scalar theory with the spontaneously broken symmetry which we leave for future study.
Although our findings indicate that f ew → N f scattering processes lie in the classically forbidden region, they tell nothing about actual value of the probability of these processes. As we already mentioned in the introduction there exist semiclassical methods to calculate the suppression exponent of the probability which, however, are difficult to apply due to singular nature of corresponding classical solutions. On the one hand, one can start, as we did in this paper, with processes N i → N f where both initial and final states contain large occupation numbers and apply the same semiclassical techniques. In this case based on the results of the semiclassical studies of baryon number violating processes [48,49] and soliton pair production [50,51] in particle collisions one can expect that corresponding classical solutions will be nonsingular. The classical solutions at the boundary of the classically allowed region obtained in the present study are expected to be close to the solutions describing transitions in the forbidden region. Using the Rubakov-Son-Tinyakov conjecture [52] the probability of 2 → N f scattering in the leading semiclassical approxi-JHEP11(2018)068 mation can be obtained from that of N i → N f process by taking the limit λN i → 0 if it exists. 3 This procedure can be viewed as a regularization to the semiclassical method of singular classical solutions. We are going to pursue this idea in future work [59].