Numerical study of multiparticle scattering in $\lambda\phi^4$ theory

We study numerically classical collisions of waves in $\lambda\phi^4$ theory. These processes correspond to multiparticle scattering in the semiclassical regime. Parametrizing initial and final wavepackets by energy $E$ and particle numbers $N_{i}$, $N_{f}$ 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). Tree level calculations performed in (3 + 1)dimensional scalar λφ 4 theory show [3][4][5][6][7] 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 [8] and violates unitarity at sufficiently large multiplicities N f ∼ 1/λ. It was observed [9], 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 multiparticle cross sections were developed in [10][11][12][13][14]. They involved singular solutions of classical equations of motion and allowed to extend the observations obtained with perturbative calculations to higher energies. Still the most reliable results were obtained in the regime of small λN f where corresponding cross sections turns out to be exponentially suppressed. Unitarity bounds [2,15,16] 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 [17][18][19], where it was found that the perturbative factorial growth is replaced by exponential suppression for λN f > ∼ 1. But still behaviour of the probability (1.1) in λφ 4 field theory at large λN f is unknown 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 number 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 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 are not exponentially suppressed at E > E * [16]. Otherwise the existence of nontrivial minimum at N i ∼ 1/λ would indicate on exponential suppression of the 2 → N f scattering for E < ∼ E * . This idea of exploring classical region was used previously for study of several different processes induced by collisions of particles. Among of them are false vacuum decay [20], the baryon number violating processes in the Standard Model [21,22] and soliton-antisoliton pair production in (1 + 1)-dimensional scalar field theory [23].
In this paper we study classical solutions describing scattering of wave packets in the unbroken λφ 4 theory. In particular, we try to approach boundary of the corresponding classically allowed region. Classical scattering of waves in relation to multiparticle production was previously studied to some extent in (1 + 3)-dimensional φ 4 model [24], (1 + 1)-dimensional abelian Higgs model [25] and non-abelian gauge theories [26][27][28]. In particular, authors of Ref. [24] constructed 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 production of many "quanta". In our study we use stochastic sampling technique to scan numerically over classical solutions describing 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 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 non-trivial classically allowed region for 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 studied energy range. We examine properties of the classical solutions corresponding to the boundaries. Obtained results give indirect indication for exponential suppression of the probability of f ew → N f processes for considered energies.
Let us note that recently an interest to the problem of multiparticle production has been renewed [29][30][31][32][33][34][35][36][37][38][39]. In particular, calculations of Refs. [35,39] performed in λφ 4 theory with spontaneously broken Z 2 symmetry, exhibited that the probability of 2 → N f scattering process near threshold is unsuppressed at λN f 1. In view of this result it would be interesting to apply the method of classical solutions used in our study to explore the classically allowed region in (E, N i , N f ) parameter space in the spontaneously broken λφ 4 theory.
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 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 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 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 position dependent interaction. Corresponding equation of motion reads and the energy related to the field configuration χ(t, r) is In what follows we solve 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 field configuration as follows 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 multiparticle scattering. In linear regime the time-dependent 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) 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 are initial and final particle numbers, respectively. For convenience we introduce the following shorthand notations For 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. [40]. Some of the results have been verified with 4-th order Runge-Kutta method with a very small time step. For convenience, we introduce 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 [41] to compute field configuration χ(t, r) as (2.7) and interaction term (2.9).

The method
In this Section we describe a method used to explore classical transitionsÑ i →Ñ f at fixed energyẼ. In particular, we are interested in solutions of 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 . 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 very mild dependence on 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 e.o.m. (2.8) (i.e. 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Ñ 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 particle number changes as much as possible at a given energyẼ. Since we solve initial value problem it is more convenient to fix initial particle numberÑ i and find 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 respectful normalization of the Fourier amplitudesf n .
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 stochastic sampling technique in combination with the simulated annealing method, see [40,42]. Namely, we generate an ensemble of the classical solutions with fixed initial particle numberÑ i weighted by probability where If large positive β and ξ are taken the ensemble will be dominated by solutions having small F thus driving their distribution towards lower boundaryÑ min f (Ẽ * ) of classically allowed region. To reach the upper boundaryÑ max f (Ẽ * ) signs of β and ξ should be negative.
To generate the ensemble (3.4) we use Metropolis Monte Carlo algorithm. We start with a solution specifying by a randomly chosen set of f n , see (3.1), which have initial particle numberÑ i . This solution has an energyẼ and final particle number N f which is found after solving of classical equations of motion. Next, we choose a few randomly chosen amplitude numbers and change corresponding amplitudes f n by small amounts, 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. Modified set of amplitudes f n is rescaled in a proper way and used to construct initial wave packet having the same initial particle numberÑ i and some value of energyẼ , see Section 3.1. Then we evolve the system forward in time far enough until it reaches linear regime and we can find final particle numberÑ f of new solution. We compute ∆F ≡ F − F using (3.5). New set of amplitudes f n is accepted with probability and used for 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 [43] to where i is 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 variation of a and β 0 . Let us note, that the stochastic sampling technique was previously used to study classically allowed sphaleron transitions [21,22] and classical soliton-antisoliton pair production in particle collisions in (1+1)-dimensional scalar field theory [23], where its efficiency in determination of the boundary of classically allowed region was confirmed. which determines the relevant energy interval for our study. We comment on this choice below. The solutions which seed our Monte Carlo search are randomly chosen with the the only condition that they haveÑ i = 1 and linearize at initial and final times. These solutions typically have very small change in particle number are situated near the linẽ N f = 1 on (Ẽ,Ñ f ) plane. Initially, we take 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 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 Fig. 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 classical evolution increases with the energy of solutions but does not exceed 12% in the chosen energy range. Let us stress that due to time reversal symmetry one can always interchange the initialÑ i and finalÑ f particle numbers in the discussion of classically allowed region in the  Fig. 3 where χ(t, r) = φ(t,r) r is shown for convenience. 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 Fig. 3. The initial wave packet propagates towards 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.
According to the used extremization procedure which involves stochastic sampling, the found classical solutions describe scattering of waves withÑ i = 1 and maximal 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 initial and final wave packets. However, we can calculate 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 evolution of the field is linear. On the right panel of the Fig. 4 we plot the instantaneous particle number for the boundary solution shown in Fig. 3. Comparing it with the time evolution of the field presented on Fig. 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 Fig. 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 smearing is less than 10 −3 . In a way similar toÑ (t) we define 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 Fig. 4 (left panel) for the solution from Fig. 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.
Nontrivial change in particle number implies a redistribution of energy between low and high frequency modes. In Fig. 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 particle number. Also we see that with increase of colliding energyẼ 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 a given lattice size and spacing. This is 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 consists mostly of nonrelativistic modes and larger space interval R is needed to obtain solutions which are linear 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 (Ẽ,Ñ f ) plane forÑ i = 1 and for the boundary solutions on parameters of the lattice, namely on the spacial cutoff R and on 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 classically allowed region for the case R = 20, N r = 600 with smaller lattice spacing and for the case R = 30, N r = 600 with larger space interval but the same lattice spacing. The results are presented in 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Ẽ > ∼ 11 which reflects appearance of new high frequency modes. Examples of the upper boundary solutions for the same set of energies as in Fig. 2 but for the case N r = 600, R = 30 are shown in Fig. 7. Comparing them with the solutions presented in Fig. 2 we see that the difference lies in 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 Fig. 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  to time translations. This symmetry is explicitly broken in our numerical setup by the choice of the initial configuration (3.1) and (3.3) at t = t i ≡ 0, which in particular implies χ(0, r 1 ) = χ(0, r 2 ) = 0 .
Still the oscillating form of the initial wave packets (see Figs. 2 and 7) and the time evolution of the particle number indicate that an approximate discrete time-shift symmetry should exists in the lattice version of the model. In terms of initial wave packets this symmetry connects the configurations which are related by time evolution and which satisfy conditions (4.3) as long as the spiky part of the initial wave train lies inside the interval [r 1 , r 2 ] . Indeed, with our numerical technique (making several runs at the same energyẼ * with different seeds) we actually find several branches of boundary solutions which are related with each other by such shifts in time. In fact, for constructing the boundaries in Fig. 1 we selected the solutions whose initial wave packets have the longest soft oscillating part while their sharp edges are placed as far as possible (for given R and initial space interval [r 1 , r 2 ]) from the origin. In Fig. 9 we plot initial field configurations (red solid lines) corresponding to another branch of the upper boundary solutions at the same values of energy as in Fig. 2. These wavepackets have spiky parts shifted closer to the origin as compared to those in Fig. 2. Also in the Fig. 9 in thin (blue) lines we show the initial field configurations from Fig. 2 evolved forward in time until the positions spiky part of the different solutions coincide. We observe perfect coincidence of the form of that part. As the spiky part of the wavepackets is responsible for the change of the particle number we expect that time-shifted solutions would have similar values ofÑ f . This is indeed the case and in Fig. 6  have 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 initial particle numberÑ i . In Fig. 11 we show the upperÑ f (E) and lowerÑ f (E) boundaries of the classically allowed region forÑ i = 0.1 and relevant energy interval 1.5Ñ i , 15.0Ñ i with the values ofẼ * defined by (4.1). We see that maximal relative difference of initial , 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 Fig. 12 we show examples of the initial wave packets and instantaneous particle number evolution for upper and lower boundary solutions forẼ ≈ 0.6. The space interval for initial wavepacket [r 1 , r 2 ] is shown by 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 interaction region results in actual change in particle number and a soft oscillating part. As in the previous case we obtain several branches of 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Ñ i = 0.1, 1.0. In Fig. 13 boundary of the classically allowed region is shown with 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Ñ i = 0.  Figure 12: Initial wave packets (left panels) and evolution of the instantaneous particle number (right panels) for boundary solutions withÑ i = 0.1 andẼ ≈ 0.6. Upper and lower panels correspond to upper and lower boundary solutions respectively.  andÑ min f (E), consist of two different branches of classical solutions. At energies lower than about 70 the boundary solutions have two distinct spiky parts. An example of initial field configuration φ(t = 0, r) corresponding to lower boundary is presented in Fig. 14 (upper left panel) forẼ ≈ 65. The evolution of instantaneous particle num-berÑ (t) as well as distribution of energy over the modes k for initial and final field configurations are also shown in this figure 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 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 Fig. 14 has δr ≈ 2.8 and δk ≈ 2.3. Comparing 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 different branch of solutions whose initial (and final) space field configuration contains already three spikes. On the right panels in Fig. 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 the two-spike solutions atẼ > ∼ 70 and three-spike solutions atẼ < ∼ 70 are 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 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 Fig. 13 by thin red (solid) line. In similar way we find continuation of three-spike branch of the solutions atẼ < ∼ 70 which is shown by this blue (dashed) line in Fig. 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   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 initial field configuration; d) the same as in c) but for final field configuration. Left and right panels correspond to the lower boundary solutions withẼ ≈ 65 and 80, respectively. extrema toÑ f . In the linear regime these solutions look as a sum of two or three single-spike wavetrains which differ from those already described by somewhat different distances between the wavetrains. In Fig. 15 Figure 15: Initial wave packets (left panels) and evolution of the instantaneous particle number (right panels) for solutions withÑ i = 10 andẼ ≈ 80 corresponding to local minima ofÑ f . Upper and lower panels correspond toÑ f >Ñ i andÑ f <Ñ i .
N f >Ñ i . The distance between the wavetrains in the initial field configuration is larger than that of for the boundary two-spike solution. Similarly, right panels show an example of solution near the lower boundary delivering a local minima toÑ f with three spike. Again we observe somewhat larger distance between wavetrains in the initial field configuration, c.f. Fig. 14 (upper right panel). These two branches of solutions are shown in Fig. 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 N 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 seed some initial field configuration which was obtained forÑ i = 1 as a boundary solution which is of single-spike type and increasing N 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 Fig. 16 we show our numerical results for the boundary of 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 Fig. 17. At the same time we find that 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 Fig. 16 we present only the envelopesÑ min  be considered as an indication on existence of some limiting boundary of 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 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 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 considered energy range E < ∼ 10N f m. In 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 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 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 Figs. 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 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 at classically allowed region for arbitrarily large values of N f and E. However, we obtain an evidence for existence of a limiting boundary the classically allowed region in (E/N i , N f /N i ) plane at large N i .
Although our results 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 [44,45] and soliton pair production [46,47] 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 [48] the probability of 2 → N f scattering in the leading semiclassical approximation 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 [55].