New method of accounting for interference contributions within a multipheriperal model

We consider an inelastic scattering of protons within the simplest real scalar model $\phi^3$ (phi-cubed). Although this model is being studied for a very long time, the problem of accounting for the interference contributions for all the possible particle multiplicities observed in experiment is not solved yet. We propose the method which is based on grouping of the interference contributions into sets in such a way that the sum of all interference contributions of each particular set can be calculated with Laplace's method. This approach allowed us to calculate all the interference contributions to the cross-sections for multiplicities up to $ n\sim 50 $ at the energy $ \sqrt{s} \sim 50$ GeV. The obtained models of the energy dependence of total $pp$ scattering cross-section and the inclusive rapidity distribution are in qualitative agreement with the experiment. We also consider the well known effect of the energy dependence of the shape of inclusive rapidity distribution, and propose an explanation of this dependence and consider it exactly as the interference effect.

We consider an inelastic scattering of protons within the simplest real scalar model φ 3 (phicubed). Although this model is being studied for a very long time, the problem of accounting for the interference contributions for all the possible particle multiplicities observed in experiment is not solved yet. We propose the method which is based on grouping of the interference contributions into sets in such a way that the sum of all interference contributions of each particular set can be calculated with Laplace's method. This approach allowed us to calculate all the interference contributions to the cross-sections for multiplicities up to n ∼ 50 at the energy √ s ∼ 50 GeV. The obtained models of the energy dependence of total pp scattering cross-section and the inclusive rapidity distribution are in qualitative agreement with the experiment.
We also consider the well known effect of the energy dependence of the shape of inclusive rapidity distribution, and propose an explanation of this dependence and consider it exactly as the interference effect.

I. INTRODUCTION
When calculating the cross-sections of an inelastic protons scattering within the multiperipheral model [4,33] it is common to make an assumption that the multi-Regge region [5,12,13,32,35] makes a main contribution to the multidimensional integral. According to this assumption, the produced particles are strongly ordered in rapidity. As a result, the Reggeization models [4,17,32,33] use the approximations so that the energy-momentum varies widely within the integration domain (see the last section in [50]). At the same time it is possible to specify without any assumptions the integration domain which makes a main contribution to the multidimensional integrals for the scattering cross-section. It has been shown (see [45][46][47][48][49][50]) that the squared modulus of the multipheriperal diagram contribution to the scattering amplitude has quite a distinct conditional maximum given that the energymomentum law is satisfied. This maximum point may be found either using any numerical method or analytically within some approximation [48]. Then one may use Laplace's method [11] to calculate the integrals for the cross-sections. The main contribution to the crosssection integral is made by the neighborhood region of the maximum point; according to Laplace's method this region is determined by the second derivatives of the logarithm of the integrand at the maximum point. This region is also very different from the multi-Regge one, since it does not force the produced particles to be strongly ordered in rapidity. In particular, when the energy √ s is close to the threshold energy for production of the n secondary particles, all the rapidities at the maximum point are approximately equal and close to zero; moreover, the distances between the rapidities at the maximum point increase logarithmically slowly [50] as the * frumle@ukr.net energy √ s grows, and decrease as the n grows. Consequently, the strong ordering of the rapidities for a given multiplicity n may appear only at the energies √ s such that the partial cross-section of the process with production of the n particles is almost zero. An absence of the strong ordering of rapidities makes it pointless to consider the multipheriperal diagrams without considering the interference contributions (which represent the various ways of joining of the identical particle lines to the diagram [49,50]). Although the values of the interference contributions may be small, their number is too large to be neglected. Moreover, it has been shown (see [49]) that the contribution from the ladder diagram is much smaller compared to the sum of the other contributions.
The problem of an accounting for the interference contribution is considered in the paper [51] and briefly in [54]. However, in these papers the interference contributions have been considered within the approximations which are used in the ABFST model (see [4]). Note that the assumption about the strong ordering of rapidities is used implicitly in the ABFST model (in order to simplify the integration limits in the integral equation for the n-particle contribution to the imaginary part of the scattering amplitude at zero four-momentum transfer). In addition, in the paper [51] a calculation of the crosssections is considered only for the cases of the small multiplicities of produced particles, where the simple direct summation of the all interference contributions is possible. However, the number of the interference contributions increases rapidly (as n!) as the number of the secondary particles n grows. Each interference contribution can be calculated quite easily with Laplace's method [45][46][47][48][49][50] even for the large number of the secondary particles n. However, it is problematic to account the large number of such contributions even for the numerical computations. Therefore, the objective of this work is to propose an approximate calculation method which would allow one to arXiv:1509.04329v2 [hep-ph] 15 Jan 2020 calculate the sum of all the interference contributions for the processes with a large number of secondary particles (n ∼ 50).
The Bose-Einstein correlations [20], or the correlations between the four-momenta of the identity particles [8,15], is almost the only interference effect discussed in the literature describing the experiment. This effect is typically considered in the context of the experimental determination of the dimension of the secondary particles radiation region [1,36,44]. A solution to this problem does not require a detailed description of the dynamics of inelastic process [9,29]. At the same time, it is the dynamics that makes the major difference between our model and the models [21,24,30,31,55,58,61,64] which have been used for description of the interference effects. Let us discuss this difference.
The main problem when accounting for the interference contributions within the proposed model is that the real-valued integrands of different interference contributions reach their maximum at different points. If all the interference contributions had a common single maximum point, the sum of these interference contributions would also have a single maximum at that point. In this case one could use Laplace's method to calculate the whole sum of the interference contributions rather than calculating each contribution separately. Moreover, as we show below, in this case the maximum of such sum would become more pronounced as n grows. As a result, the accounting for the interference contributions would not be a problem for an arbitrarily large number of secondary particles. Meanwhile, in our model the different interference contributions have the different maximum points. Consequently, the sum of all interference contributions may have multiple maximum points depending on the energy √ s and the number of secondary particles n. The number of these maximum points may become so large that the time required for the numerical calculations increases critically. Even when the sum of all interference contributions have only a few maxima, they can be not pronounced enough under certain √ s and n to limit to the quadratic terms in Laplace's method. It means that one should take into account the higher order terms in the Taylor expansion of the integrand logarithm around the maximum point. By contrast, the models considered in the papers [3, 21, 24, 39-41, 58, 61, 64] assume that the secondary particles are produced by the independent sources, and the Bose-Einstein correlations occur only when the particles propagate from the source to the detector. In this case different interference contributions are the Fourier transforms of the same function but at different points. As a result, the different interference contributions reach the maximum at the same point, which means that the dynamics in these papers is fundamentally different from the dynamics considered herein. The same is true for the models [3,24,30,31,61], in which the problem of the field interacting with the fixed or random (Langevin) source is considered instead of the dynamics of interacting fields. In this case the secondary particle momentum distribution for an arbitrary number n of the secondary particles is determined by the single function (i.e. by the Fourier transform of the external source); thus the different interference contributions reach the maximum at the same point. Such models also have an obvious problem associated with the violation of the energy-momentum conservation law caused by the space-time translation symmetry breaking. The latter is evidenced by the fact that, according to these models, as well as to the models based on the multi-Regge kinematics [23], the multiplicity of secondary particles has a Poisson distribution [24]; thus the particle production processes [26] at different regions of the phase space are independent of each other; consequently, the probability of the production of an arbitrarily large number of secondary particles is not zero. In addition, the models with an external source imply the source averaging with the corresponding density matrix [7,18,43,57,58]. However, the density matrix is actually postulated rather than determined with the time evolution operator obtained within some dynamical model. From the above one may conclude that interference effects in the proton-proton scattering must be evident not only in the measurement of the Bose-Einstein correlations, but in any measured quantity. In this paper we argue that the dependence of the inclusive rapidity distribution shape [6,22], [2,14,19,53], [16](p. 590) on energy may be considered as an interference effect. where the shape of the inclusive rapidity distribution [6,22], [2,14,19,53], [16](p. 590) depends on energy √ s, may be considered as the interference effect. The experimental data for the inclusive rapidity distribution is described by the different models [10,25,27,28,34,37,59,60,62,63] which however do not take into account for the interference effects. The manifestation of the interference effects in the inclusive pseudorapidity distribution is considered in the papers [40,61], though using the above mentioned approximations in which the different interference contributions reach maximum value at the same point. By contrast, in this paper we show that the behavior of the inclusive rapidity (pseudorapidity) distribution shape under the energy √ s growth is associated with the change of the distances between the maximum points of the different interference contributions. The mentioned difference between the dynamics considered in this paper and the one considered in papers [40,61] becomes apparent from the comparison of the shapes of inclusive rapidity distributions obtained within these models. In order to focus on the problem of accounting for the interference contributions let us start with the simplest model φ 3 .

II. LAPLACE'S METHOD, φ 3 MODEL
We consider the so-called ladder diagrams for the inelastic pp scattering within the real scalar φ 3 model, assuming that the masses of the secondary particles are equal to the pion mass m π . All the physical quantities used in this work (four-momenta, masses, energies...) are nondimensionalized by the pion mass m π .
Let S n be the set of all possible permutations of N n = {1, 2, . . . , n}. Each diagram of the form Fig. 1 is characterized by a corresponding permutation π ∈ S n . Assuming the one-line notation for the permutations, π (i) denotes the index of the secondary particle which is joined to the i-th vertex in the corresponding diagram. Each permutation π ∈ S n specifies the diagram of the form Fig. 1, and thus represents an analytic expressionthe additive contribution to the scattering amplitude a n P 1 , P 2 , P 3 , P 4 , p π(1) , ...p π(n) , where P 1 , P 2 are the four-momenta of the colliding particles in the initial state, P 3 , P 4 are the four-momenta of the scattered particles in the final state, p π(1) , ...p π(n) are the four-momenta of the secondary particles in the final state. According to the Feynman diagram technique, the scattering amplitude T n of the considered process is the sum of all expressions a n corresponding to the elements of S n T n (P 1 , P 2 , P 3 , P 4 , p 1 , ...p n ) = = π∈Sn a n P 1 , P 2 , P 3 , P 4 , p π(1) , ...p π(n) . (2) For convenience, we introduce the following notations. We will consider the scattering process in the center-ofmass reference frame using the right-hand coordinate system in which z axis is codirectional with P 1 . Here the expressions a n (P 1 , P 2 , P 3 , P 4 , p i1 , ...p in ) depend on the real particles four-momenta which satisfy the following equations where M i is the mass of the i-th incident particle in the initial state, and m is the mass of the secondary particle. The components of P 1 and P 2 in the chosen reference frame are (P 1 ) x = (P 2 ) x = 0, (P 1 ) y = (P 2 ) y = 0 and (P 1 ) z = −(P 2 ) z . Thus P 1 and P 2 can be uniquely determined by the single quantity √ s = (P 1 ) 0 +(P 2 ) 0 -the total energy of the colliding particles in the center-of-mass reference frame. We use the components (p i ) x , (p i ) y of four-momenta and the rapidities y i of the secondary particles as the independent variables which uniquely determine all the four-momenta p i . The rapidities y i are determined from the following expression We also introduce another two independent variables P a x and P a y using the following expressions It follows from the equations (3) and the energymomentum conservation law that all the components of P 3 and P 4 can be uniquely determined by specifying the √ s, P a x , P a y and all the y i , (p i ) x , (p i ) y . Finally, instead of 4 (n + 4) variables (the arguments of a n ) we have 3n + 3 independent variables, so the expression (1) may be rewritten as a n √ s, y 1 , . . . , y n , (p 1 ) x , . . . , (p n ) x , It is convenient to consider the introduced variables (except of √ s) as the components of the some vector Finally, it is also natural to introduce the notation for permutations of the components of X. To this end, to each permutation π ∈ S n we assign a linear operator π : R 3n+2 → R 3n+2 which can be defined aŝ or the same: (πX) in+j = X in+π(j) for i = 0..2 and j = 1..n, while (πX) i = X i for i = 3n + 1, 3n + 2. One can see that according to the definition (8) ofπ, we associate the operatorπ 2π1 with the composition of permutations whereπ i is the operator associated with π i . Hence, using the introduced notations, the expression (2) can be rewritten in the following way Now let us consider the expression for the partial crosssections σ n of inelastic scattering, and rewrite this expression using the notations introduced above.
to integrate out the energy-momentum conserving delta function in (11) and then change the variables of integration. The Jacobian of the transformation and the multipliers appearing as the result of delta function integration can be included in the scattering amplitude T n . This procedure described in detail in [50]. As a result, the integral (11) takes the following form where R (n, √ s) = (2π) 4 /4n!I. Now, let us consider the problem of calculations of multidimensional integrals of the form (12).
A. Calculation of the partial cross-sections within φ 3 model Let us substitute the expression (10) for the scattering amplitude T n ( √ s, X) in the expression of the partial cross-section (12).
Note that one of the two sums over S n in (13) can be calculated by renaming the variables of integration.
The summands in (14) are called the interference contributions to the cross-section. Each particular permutation π specifies some interference contribution which can be represented by the so-called cut diagram (Fig. 2). The way in which the vertices from the both parts of a cut diagram are linked is specified by the permutation πaccording to (14) the π (i) is the index of vertex in the left part of the diagram which is linked with the i-th vertex in the right part of the diagram. Each particu- Figure 2. Cut diagram representing the interference contribution associated with particular permutation π such that π(1) = n, π(2) = 1,..., π(n) = 2.
lar interference contribution can be calculated using the Laplace's method (see [50]). This method allows one to easily calculate a multidimensional integral provided its integrand has a single maximum point. Indeed, it has been shown (see [50]) that the function a n ( √ s, X) has the single maximum point X (0) at fixed √ s.
At the same time the function a n ( √ s,πX) has the single maximum pointπ −1 X (0) , whereπ −1 is the operator (8) associated with permutation π −1 ∈ S n inverse of π; therefore the product a * n ( √ s, X) a n ( √ s,πX) also has a single maximum point. As a result, the Laplace's method can be applied to the calculation of the integrals in expression (14). The calculation of the total cross-sections with Laplace's method even within the simplest φ 3 model allows one to obtain the qualitative description of the experimental data (see [45][46][47][48][49][50]). Moreover, we plan to use Laplace's method to calculate the cross-sections within the multiparticle fields approach which is based on QCD (see [42,56]).
However, the number of the summands in (14) is n! which increases dramatically as the energy grows. It becomes impossible to account all of them by calculating each interference contribution separately. So we propose the method which makes it possible to account all the interference contributions for the processes with up to 50 secondary particles.

III. THE MAIN IDEA OF THE PROPOSED METHOD
Let us consider a simple example to illustrate the main idea of the proposed method. We say that the maximum points of two functions, each having a single maximum point, are close if the sum of these functions also has a single maximum point. Let φ a (x) = exp − (x − a) 2 be the function parametrized by the parameter a and has single maximum point x = a. Now let us consider the sum (16) for different values of a, and the partial sums ψ − (x) = φ −a (x) + φ −a/2 (x) and ψ + (x) = φ a/2 (x) + φ a (x). It is easy to see that in the case of a = 0 the function g 0 (x) has the single maximum point x = 0. Then, the distances between the maximum points of summands of g a (x) grow with the parameter a. For the values of the parameter a up to 1 each of the functions g a (x), ψ − (x) and ψ + (x) has a single maximum point, as can be seen in Fig. 3.a and Fig. 3.b. With further growth of a the distances between the maximum points of the summands φ of g a (x) become so large that the single maximum of the function g a (x) splits into two separate maxima (Fig. 3.c). At the same time, the partial sums ψ − (x) and ψ + (x) still have the single maximum points.
The similar effect may be observed for the scattering amplitude T n ( √ s, X) with the energy growth at fixed n > 1. The function T n ( √ s, X) has a single maximum at low energies which then splits into a few separate maxima as the energy grows. To show this, let us consider the features of the maximum point X (0) of the function a n ( √ s, X). At this point the transverse components of the particles momenta are equal to zero: i ) y = 0, P a x = 0, P a y = 0, while the rapidities of the secondary particles form the following arithmetic progression (see [50]) As seen from (17), y i = −y n−i . Now consider the expression for the scattering amplitude T n ( √ s, X) i.e. the sum (10) with taking into account the mentioned features of X (0) . As we have already shown in the previous section, each summand of this sum (i.e. the function a n ( √ s,πX)) is associated with some permutation π and reaches its maximum at the single pointπ −1 X (0) . The distance ρ between the maximum points of two such summands (associated with π 1 and π 2 ) is defined by the following expression Here we have taken into account that the X i = 0 for n < i < 3n + 2. Given the fixed number of the secondary particles n, one can see from the (18) that the common difference |∆ y | of the rapidities arithmetic progression tends to zero as the energy √ s goes to the threshold value 2M + n. In this case all the summands are close so the whole sum (10) has a single maximum point; it means that one can apply Laplace's method to calculate the integral (12), thus taking into account all the interference contributions in the simple way. However, the |∆ y | slowly increases as the energy grows; and at some point the distances between the summands in (10) become so large that the whole sum acquires several maxima. In other words, the single maximum of the scattering amplitude T n ( √ s, X) splits into a few separate maxima as the energy grows. We have already faced the similar effect for the one-dimensional function in the beginning of this section. In this case Laplace's method cannot be applied to the calculation of the integral (12) since its integrand has several maxima. One could still consider the sum (14) and calculate each interference contribution with Laplace's method, but this way it is impossible to calculate all the interference contributions for the processes with a large number of secondary particles.
The main idea of the proposed method of accounting for the interference contributions is to apply Laplace's method not to the each interference contribution separately, but to the sums of the interference contributions whose integrands have the close maximum points. The integrand of such sum has a single maximum point so the whole sum can be calculated with Laplace's method. Although the scattering amplitude T n ( √ s, X) acquires several maxima as the energy grows, we can split the sum (10) into subsums each having a single maximum point. To this end we group the neighboring vertices of the diagram into k groups. Then we group all the permutations π ∈ S n into subsets I 1 , I 2 , . . . in such a way that any two permutations belong to the same subset if and only if they specify the diagrams in which the lines of the secondary particles are joined to the same groups of the vertices regardless of the order within the groups. So the diagrams specified by the permutations belonging to the same subset I i differ only in the order of joining inside the groups of vertices. Therefore, according to (17) and (18), we expect the functions associated with such diagrams to have the close maximum points. Thus, for each subset I i we consider the sum b i ( √ s, X) = π∈Ii a n ( √ s,πX). For any given values of n and √ s we select such a number and the sizes of the groups that make each of the functions b 1 , b 2 , . . . have a single maximum point. As a result, one can use Laplace's method to calculate the sum of the interference contributions associated with each subset I i ⊂ S n rather than calculating the interference contribution associated with each permutation π ∈ S n . Such approach reduces greatly the amount of the calculations and allows one to account for all the interference contributions for the processes with the high production of the secondary particles (up to 50 secondary particles at 72 GeV at present). The specific implementation of the proposed idea will be described in the next section in detail.

A. Grouping the vertices of the diagram
Let us group the neighboring vertices of the diagram Fig. 1 into the k groups (non-empty sets) each containing n i vertices. For this purpose we consider the sequence of sets G 1 , G 2 , . . . , G k each specifying a corresponding group. The set G i contains the indices of vertices that are grouped into i-th group.
for 1 i k, where n i is the number of vertices in the i-th group. The lines of the secondary particles in the connected diagram specified by a permutation π are distributed among the groups of the vertices in some way depending on the permutation. Thus for any permutation π ∈ S n we can also consider the sets π (G 1 ) , π (G 2 ) , . . . , π (G k ), where π (G i ) = {π (ν) | ν ∈ G i } is the set containing the indices of the secondary particles which are joined to the vertices of the i-th group in the diagram specified by the permutation π. It allows us to introduce the equivalence relation ∼ on the permutations set S n in the following way Which means that any two permutations π i and π j are equivalent if they specify the diagrams in which the external lines of the secondary particles are joined to the same groups of the vertices regardless of the joining order inside the groups (Fig. 5). As a result, the permutations set S n may be split into the equivalence classes [π] = {π i ∈ S n | π i ∼ π}. In other words, the introduced equivalence relation ∼ provides the partition S n / ∼ of the underlying set S n , where S n / ∼ is the quotient set of the permutations set S n by ∼ i.e. the set contains all the equivalence classes [π].
S n = [π]∈Sn/∼ [π] (21) Taking into account (21), the expression (10) may be rewritten as Now let us consider the auxiliary function A n ( √ s, X) where [ ] is the equivalence class containing the permutations equivalent to the identity permutation such that (i) = i for all i ∈ [1 . . . n]. The maximum distance between the maximum points of summands in (23) may be significantly reduced compared to the summands of (10) by increasing the number of groups k and decreasing their sizes n 1 , n 2 , . . . n k . Therefore, for any given values of n and √ s one may select such a number of the groups k and sizes of the groups n 1 , n 2 , . . . n k that make function A n ( √ s, X) have a single maximum point. It may be shown that for any permutation π i ∈ [π] there exists a permutation e ∈ [ ] such that π i = π • e, where π is the representative (any member) of the class [π]; taking this into account and also (9), let us represent the scattering amplitude (22) in terms of A n ( √ s, X) whereπ is the permutation operator (8) associated with the representative of [π]. As a result, the integral (12) may also be written in terms of A n ( √ s, X).
One of the two summations in the integral (25) may be calculated by renaming the variables of integration as it was done between (13) and (14). Note that in this case the number of the elements in the set of summation S n / ∼ is n!/n 1 !n 2 ! . . . n k !.
Taking into account that the function A n ( √ s, X) has a single maximum point, Laplace's method may be used to calculate the summands of (26). The number of the summands in (26), which is equal to |S n / ∼ | = n!/n 1 !n 2 ! · · · n k !, may be significantly reduced by selecting the number of group k such that k << n. Each summand of (26) can also be represented by the corresponding cut diagram Fig. 6. Hence, according to (26), it is sufficient to consider only a single permutation from each equivalence class [π] rather then all permutations of S n in order to calculate all the interference contributions to the cross-section σ n . Moreover, we will show that there are many similar summands in (26) which may be easily calculated by the corresponding weight-factors.
Let us consider a summand of (26) associated with a class [π i ]. Then change the integration variables: X → eX, where e is a permutation belonging to the class [ ]. Note that the function A n is symmetric with respect to such transition: A n ( √ s,êX) = A n ( √ s, X).
According to the definition of permutation operator (8), the operatorπ l =π iê is assigned to the permutation π l = e • π i . It may be shown that the permutations π i and e may be selected such that π i and π l = e • π i belong to the different equivalence classes [π i ] and [π l ]. As a result, one can see from (27) that there are two similar contributions to the sum (26) associated with the different equivalence classes.
In other words, the definition of the function A n implies that both parts of the cut diagram Fig. 6 are symmetric with respect to the order of joining inside the groups of vertices. Consequently, the permutations which establish the same number of connections between the groups of vertices differ only in the order of linking inside the groups and therefore represent the similar contributions to the sum (26). To clarify this statement, let us consider the number of connections m ij between the i-th group in the left part and the j-th group in the right part of the cut diagram (6) specified by a permutation π.
where | · | denotes the number of elements in the set (i.e. the power of the set). In this way, each equivalence class [π] may be characterized by a matrix m of size k×k whose elements are defined by (28). If any two classes [π i ] and [π j ] are characterized by the same matrix m, then they specify the equivalent cut diagrams, which mean that they correspond to the similar contributions to the sum (26) because one of these diagrams can be obtained from another by rearranging the vertices inside the groups.
In order to consider all the unique summands of (26), we consider the set M of all possible matrices m of the size k × k such that the sum of the elements of i-th row (column) is equal to the number of vertices in the i-th group.
Each matrix m represents a summand in (26). All the summands similar to this one can be accounted by weight factor W m uniquely determined by the elements of the matrix m. To calculate the summand specified by the matrix m, one needs to link the vertices of the cut diagram in a way specified by the matrix (Fig. 7). As a result, one obtains the permutation π m such that |π m (G i ) ∩ G j | = m ij , which makes it possible to calculate summand associated with the class [π m ]. The weight factor W m may be calculated by counting the number of different equivalence classes [π] characterized by the same matrix m. The number of such classes is equal to the number of all possible ways in which one can group the elements of the sets G 1 , G 2 , . . . , G n into .k]; so that each set G j contains exactly m ij elements from the set G i .
Finally, the expression for the partial cross-section may be rewritten as Note that the number of summands in (31) may be much smaller compared to the sum (11). At the same time, Laplace's method may be used to calculate each particular summand of the sum (31).
It remains now to consider in detail all the components necessary to apply Laplace's method within the proposed approach for calculation of interference contributions.
First we consider the feature of the maximum point of the function A n ( √ s, X). We have already mentioned the features of the maximum point of a n ( √ s, X) in Section III. Since the function a n ( √ s, X) grows as the (p ⊥ ) i → 0 and P a ⊥ → 0, the function (23) reaches the maximum value at a point χ (0) where (p ⊥ ) i = 0, i = 1..n and P a ⊥ = 0; so χ (0) i = 0 for n < i 3n + 2. As we mentioned above, we group the vertices of the diagram in a way for the function A n ( √ s, X) to have a single maximum point. Taking also into account that A n ( √ s,êX) = A n ( √ s, X) for any e ∈ [ ], one can conclude thatê for any e ∈ [ ]. This important feature of χ (0) allows one to calculate the value of function A n together with its 2-nd derivatives at the maximum point in a simple way.
Let us consider the first n components of the vector of derivatives ∂A n /∂X i of A n at the maximum point χ (0) assuming that χ (0) satisfies the equation (32). We denote these components by ∂A n /∂y i where g (i) is the index of the group containing the index i (so that i ∈ G g(i) ), e −1 is the inverse permutation of e, i.e. such that e −1 • e = . The derivatives ∂A n /∂ (p x ) i and ∂A n /∂ (p y ) i are calculated in the same way as ∂A n /∂y i . For the last two components of ∂A n /∂X i we have where i = 3n + 1, 3n + 2. The matrix of the second derivatives ∂ 2 A n /∂X i ∂X j √ s, χ (0) may be split into blocks with respect to the variables: y i and (p x ) x , (p x ) y . Then each block may be considered separately. Let us start with the block containing the derivatives ∂ 2 A n /∂y i ∂y j √ s, χ (0) which we denote briefly by D ij .
In order to calculate the sum (36), we denote the derivatives ∂ 2 a n /∂y i ∂y j √ s, χ (0) by d ij and consider the following cases keeping in mind that χ (0) satisfies the equations (32): 1. g(i) = g(j) and i = j 2. g(i) = g(j) and i = j The rest of the blocks containing the derivatives may be calculated in the same way as (37) - (39).

IV. RESULTS
We applied the proposed method within the phi 3 model to calculate the energy dependence of the pp scattering total cross-section and the inclusive rapidity distribution for inelastic pp scattering. The models presented in Fig. 8 and Fig. 9 are obtained at the value of effective coupling constant L = 29 (the expression for L is the same as in [50], page 876).
As we can see from Fig. 8 and Fig. 9, comparison with the experimental data shows that the proposed approach allows us to obtain the theoretical predictions that are in qualitative agreement with the experimental data even within the simplest model. The following expression has been used for the calculations of inclusive rapidity distribution where N ( √ s) is for the maximum available number of the secondary particles (each having mass m) which may be produced as the result of inelastic scattering at a fixed energy value √ s.
It is of interest to consider the peaks behavior in Fig. 9. At low energies one can observe the single peak at rapidities (pseudorapidities) close to zero. However, this peak becomes less pronounced and transforms into the so-called rapidity plateau as the energy √ s grows. With further energy growth the plateau splits into the two separate peaks at non-zero points located symmetrically about zero; The effect described above may be explained by analyzing the behavior of the maximum points of the scattering amplitude absolute value |T n ( √ s, X) | 2 with the energy growth. Note that the function |T n ( √ s, X) | has the same maximum points as T n ( √ s, X) since the scattering amplitude (10) within φ 3 model is the real-valued function, except for the constant complex factor which may be neglected.
T n √ s, X = π∈Sn a n √ s,πX Let us consider (41) at the values of energy close to the inelastic threshold energy √ s 2M + m. In this case, according to (17), for any n = 1, .., N ( √ s) the rapidities of the secondary particles y Hence, the most probable inelastic processes at low energies are the processes in which all the secondary particles have near-zero-rapidities. I.e. the mean value of the vector of secondary particle rapidities at low energies is equal to zero. As a result, the inclusive rapidity distribution has the single pronounced maximum at zero. The continuous growth of the energy √ s increases the distances between the maximum points of the summands in (41). As a result, the single maximum of T n ( √ s, X) becomes less pronounced, i.e. the absolute values of eigenvalues |λ i | of the Hessian matrix (with respect to y) of the function T n at the maximum point ζ (0) decrease and tend to zero. It means that the variances of the rapidities of the secondary particles are increasing with energy growth while the mean values remain zero. Accordingly, the single maximum of the inclusive rapidity distribution transforms smoothly into the so-called rapidity plateau. The further energy growth makes the maximum points of the summands in (41) move away from each other so far that the single maximum of the T n ( √ s, X) observed at lower energies splits into a few separate maxima. As a result, the eigenvalues λ i become positive while the first derivatives are still equal to zero, which indicates that there is a local minimum of the function T n ( √ s, X) at the point ζ (0) . It means that as the energy grows, the processes in which the secondary particles have nonzero rapidities become more probable than the processes in which all the secondary particles have zero rapidities. Hence, the maxima of the inclusive rapidity distribution move away from zero. Moreover, the maximum points of the function T n are symmetrically distributed about zero due to the symmetry of the considered physical system with respect to the inversion of the collision axis. As a result, the maxima of the inclusive rapidity distribution at high energies are also symmetrically distributed about zero, which is in agreement with the experiment.

V. CONCLUSIONS
We developed the method of accounting for the interference contributions to an inelastic scattering crosssections. Using this method, we obtained the models of energy dependence of the total cross-section and rapidity distribution for the proton-proton scattering at energies up to 72 GeV and multiplicities up to 50. This models reproduce the experimental data only qualitatively, which can be explained by the fact that the modeling has been performed within the φ 3 (phi-cubed) modelthe simplest dynamical model with the scalar field.
The obtained results indicate that it would be appropriate to use the proposed method with a more complex model, for instance, with the multi-particle-fields model based on QCD [56].
In addition, the idea of the proposed method allowed us to analyze the energy dependence of the shape of inclusive rapidity distribution and to propose the physical explanation for this effect.