Laplace's method for elastic scattering diagrams within multi-particle fields model

We apply the multi-particle fields model to calculate the differential cross-section d{\sigma}/dt of elastic proton-proton scattering. This problem includes the calculation of multidimensional integrals arising from the loop Feynman diagrams. We demonstrated how these integrals can be reduced with Laplace's method to one- and two-dimensional integrals which can be calculated numerically. The obtained result qualitatively describe the minimum in differential cross-section dependency d{\sigma}/dt(t).


I. INTRODUCTION
The process of proton-proton elastic scattering is the subject of both theoretical and experimental research for the last few decades. The humanity have accumulated a lot of experimental data [1], but we still not have a dynamic theory from the first principles [2] yet. The descriptions of such processes are usually built with either the phenomenological approaches [3] which are different variations of Regge theory [4,5], or with the geometrical models [6][7][8]. All these approaches are based on some assumptions that are not the corollaries of fundamental physical principles, but may even be in conflict with these principles. In particular, the Reggeized models [9][10][11] use the assumption that multi-Regge region [12][13][14][15] provides the dominant contribution to the integrals for observables. This region of the phase space contains the points corresponding to the significantly different values of energy-momentum of the secondary particles in the final state of the scattering [16], which violates the energymomentum conservation law.
Among the various experimental data regarding the proton-proton(antiproton) scattering the most theoretical investigations have centred on the elastic scattering [17]. The later is justified by expectation for the description of elastic processes to be significantly simpler compared to the description of inelastic ones. In particular, the scattering amplitude for elastic process is the function of just two Lorentz-invariants [4]. The description of inelastic processes is then obtained from the description of elastic processes with the Abramovsky-Gribov-Kancheli (AGK) cutting rules [18][19][20][21]. However, this approach requires additional assumptions regarding the form of multi-Regge vertex [22,23] and introduction of unobservable quantity such as cross-section of scattering with an exchange of a certain number of Reggeons. All these assumptions are directed to establish the elastic scattering amplitude dependence on Mandelstam variable s. At the same time, the dependence on another variable t remains undefined. In the simplest one-Reggeon model tdependence is included in both Regge pole trejectory and residue. While the trajectory is defined by the masses of t-channel resonances, the residue remains completely undefined. In the context of Regge diagrams this residue * nata.podolyan@gmail.com is explained as the product of vertices corresponding to the interaction of Reggeon with the hadrons undergoing elastic scattering [24][25][26]. The certain t-dependence of such vertices is postulated in such a way that allows one to describe the experimentally measured differential cross-section dσ/dt at low values of t. This dependence is usually postulated to be exp −R 2 |t| , where R is a fitted parameter also known as the Regge radius of hadron [24]. As being said, such assumption provides the desirable t-dependence of the differential cross-section dσ/dt of elastic scattering at low t which is close to linear one in logarithmic scale applied to the cross-section axis [27]. However, this dependence is not linear on the whole range of measurements, but is non-monotonic and has maxima and minima [28][29][30][31]. The models with the multi-Regge exchange of simple poles and the model of quasi-eikonal multi-Regge vertices [27,32,33] are failed to explain these features. Indeed, this non-monotonic behavior has been described within the additive quark model [25] due to the interference contributions from processes with an exchange of various number of Reggeons. As for the phenomenological approaches, this dependence is also reproduced within the models with a Regge multiple poles [34][35][36]. However, these models have high uncertainty in t-dependence, since the increasing in the order of pole requires increasing in the number of terms in Laurent series whose forms are also need postulating.
The non-monotonic dependency of differential crosssection dσ/dt(t) was also reproduced in the models with simple Regge poles of various signatures and with the two-Reggeon cuts [37].
The present research is devoted to the description of effect of non-monotonic t-dependence of the elastic scattering differential cross-section dσ/dt built upon a purely dynamical model. In other words, we use the model that is based on the fundamental physical principles, beginning from the Lagrangian and its corresponding dynamical equations, quantization, and solution of the equations describing the dynamics of corresponding relativistic quantum system in Fock space. We do that in the framework of the multi-particle fields model.

II. PROBLEM STATEMENT AND LITERATURE REVIEW
We described the model of multi-particle fields in [38][39][40] and showed that hadrons can be considered as the arXiv:2101.04755v3 [hep-ph] 6 Jan 2022 quanta of field defined on the subsets of simultaneous events. These subsets are extracted from the tensor product of two (for mesons) and three (for barions) Minkowski spaces. The codomains of these fields are invariant subspaces of the tensor product of bispinors, where the scalar (for mesons) and bispinor (for protons) representations of Lorentz group act. The trivial representation of SU c (3) group acts also on these subsets, which corresponds to the fact that hadrons have no color. As usually, we build a Lagrangian for these fields and apply gauge principle. The local SU c (3)-invariance of Lagrangian is provided by the gauge fields transforming by the tensor representation of the Lorentz group and internal symmetry groups. In this case, the common way of providing the gauge invariance through the introduction of covariant derivatives is the special case of the multi-particle approach for obtaining the local invariance. We have shown [38][39][40] that these tensor fields can be used for description of the creation and annihilation processes of glueballs, i.e. the bound states of confined gluons. The same fields provide both the confinement of quarks within the protons and mesons, and interaction between the quarks of different protons through the exchange of glueballs. As the result, we have the dynamical model for the threeparticle bispinor fields that interact via the two-particle glueball fields. The calculation of observable quantities within this model can be performed with Feynman diagram technique, since such nonperturbative effects as the confinement of quarks and gluons are already accounted for in the internal dynamics of multi-particle field quanta. The non-zero masses of glueballs are also obtained from the dynamical equations in a natural way. This leads to the finite value of the elastic scattering amplitude at t = 0 due to the strong interaction. It is known that the total proton-(anti)proton scattering cross-section is finite after exclusion of electromagnetic interaction [8,41]. As the result of optical theorem [4], the scattering amplitude must have a finite value at t = 0. This finiteness is simply postulated in the mentioned phenomenological models. So it is necessary to develop a model that would describe the mentioned features of scattering amplitude without simply postulating them. In QCD perturbation theory one have an infrared singularity for the scattering with a massless gluon exchange. It is clearly problematic to use perturbative QCD at t = 0, where finiteness of the elastic scattering amplitude may be caused by nonperturbative effects. However, it is still unclear which nonperturbative effects and how they lead to the finiteness of the scattering amplitude. Meanwhile, in the model of multi-particle fields it is the non-zero mass of the glueball (which is actually the consequence of nonperturbative effects) that leads to the finite value of the scattering amplitude.
The qualitative description of the inelastic scattering may be done considering tree-level diagrams only [42]. However, the tree-level diagrams (see Fig. 1) are not enough to describe the elastic scattering. Nevertheless, the non-monotonicity of the cross-section dependence appears even in such a simple model, except that it reproduces only the minimum, but not maximum. Pole tree-level diagrams of elastic proton-proton scattering. P1, P2 are the four-momenta of initial protons, P3, P4 are the four-momenta of outgoing protons. Double lines correspond to the glueballs (bound state of gluons). The minus between the diagrams reflects the fact that protons are described by Fermi-Dirac statistics.
Therefore, the description of experimental data requires the calculation of more complex loop diagrams.
The calculation of such diagrams is reduced to the calculation of multidimensional integral over virtual fourmomenta. Similar integrals have been considered in [43] with an application of Feynman parametrization, where, however, only the location of the singularities was studied but not the calculation method. The Feynman identity simplifies the integrand due to the fact that instead of the product of multiple Feynman denominators we obtain the power of single denominator, while, unfortunately, it also leads to the essentially more complex integration domain. For this reason, here we consider an approximate method of calculation starting from the expression for the loop diagram contribution to the elastic scattering amplitude. There is also another approach to the similar integrals [44] based on the eikonal approximation. However, such approximation significantly changes the pole structure of the integrand. While the eikonal approximation allows one to reduce the dimension of the integral for the eikonal, such approach does not solve the problem of calculating the integrals themselves. The main problem here is to calculate the limit of multidimensional integral as the parameters that shift the integrand poles tend to zero. Note that taking the limit inside the integral is not allowed in this case, because the poles will move inside the integration domain so the integral will diverge. This puts restrictions on the application of numerical integration methods, because taking the limit in the end requires the calculation at small parameters, which makes the poles close to the integration domain and thus complicates the numerical calculations.
To solve the outlined problem, we use the Laplace's method [45] which worked well for the description of inelastic scattering processes [46][47][48][49]. This method allowed us to effectively calculate the integrals with dimension up to 100. However, in those works the limit could be taken before the integration, because the poles of integrands were outside of the integration domain, which is not the case here. Now let us figure out how to apply Laplace's method in the present situation. We know that the problem arises from the Feynman denominators that correspond to the lines of a diagram. Let us number these lines in an arbitrary order and put the expression of each denominator in the form (z a − iε) −1 , where a is the line number, z a is the expression corresponding to the a-th line, ε is the parameter which should be made zero after the integral is calculated. Selecting a subset of k lines we can equate the corresponding expressions z a to zero and consider the obtained system of equations. If the obtained system of equations is consistent, it defines the subset of integration domain where the integrand is equal to ε −k . As a result, the major contribution to the integral comes from the region in which the greatest number of denominators is equal to zero. Let us denote this number by l and the total number of integration variables by n. From the corresponding system of equations we can express the l variables through the rest n−l ones. Then it is convenient to change the first l integration variables. These new variables are the deviations from the values that satisfy the equation system for l denominators. It means that the absolute value of the integrand now has a distinct maximum at zero values of the first l variables regardless of the next n−l variables. Then we can apply the Laplace's method to integrate over the first l variables. As we have already mentioned, the absolute value of the integrand at the maximum point is equal to ε −l . At the same time, when applying Laplace's method, there comes the gaussian integral which leads to the factor ε −l . It removes the ε from expression and allows one to turn ε to zero before the integration. Then the obtained integral can be calculated using numerical methods.
In the present paper we apply the described idea to calculate the sum of simplest single-loop diagrams for the elastic scattering of protons (Fig. 2).
We calculate the elastic proton-proton scattering differential cross-section within the multi-particle fields approach. We consider the contributions of the diagrams in Figs where g is the effective coupling, M p and M G are the masses of the proton and the glueball respectively -are the solutions of the Dirac equations, γ a s4s2 , γ b s3s1 -are the elements of the Dirac matrices, and t ab -is the tensor whose components are defined by where a, b = 0..3.
We consider the problem in the center of mass reference frame with the right hand coordinate system whose z axis is oriented along P 1 , and the x axis is perpendicular to z and lies the plane containing vectors P 1 and P 3 . All quantities are expressed in the units of the proton mass, M G and g are considered as the model parameters.
First, we need to calculate the components of the tensor t ab . After the tensor components t ab have been calculated, we use the squared absolute value of the scattering amplitude to calculate the elastic scattering differential cross-section dσ/dt taking into account the contributions from the pole diagram from Fig. 1.
To calculate the integral in (2) we use Laplace's method. According to the calculation method described in the previous section, for a small value of ε we determine which region makes the major contribution to the integral. If the ε is small and non-zero, the absolute value of the integrand reaches its maximum in a region where the maximal number of real parts of the denominators (i.e. the parts that do not contain ε) turn to zero. Let's denote the real parts of the denominators in (2) as Taking into account (3), expression (2) can be rewritten as FIG. 3. The coordinates k 1 and k 3 of the vector k in initial coodrinate system and its coordinates k 1 and k 3 after transformation of the coordinate system.
The squared absolute value of the denominator in (2) is the sum of the real part squared and the ε 2 . When the real part in (3) is equal to zero and ε is non-zero, then the integral converges and has the maximum. Assume that one of the denominators (3) is zero, which yields some subset of the integration domain where the real part of the denominator is zero. Again, considering the integration over this subset, the main contribution will be provided by that region where the real parts of some of the other denominators turn to zero. So a natural question to ask is how many expressions in (3) can be turned to zero at the same time? It has been shown [50] that either the first pair of expressions z 1 and z 2 corresponding to the horizontal lines in Feynman diagram (Fig. 2a), or another pair z 3 and z 4 corresponding to the vertical lines, can be turned to zero simultaneously. Thus, both the horizontal and vertical lines cannot be turned to zero at the same time. Next we present the calculations of each denominator in more detail.
Let us consider the tensor (2) and denote the numer-ator as f ab k 0 , k = (k a + 2P 2a ) (2P 1b − k b ). Taking into account that in the center of mass reference frame P 2 = − P 1 , the denominator z 2 in the tensor (2) may be represented as follows Tensor (4) contains now 7 non-zero terms, and each of them has to be calculated separately. Note that they differ only in numerators and have the same denominators, which means they all can be calculated in the same way. It is then convenient to change the coordinate system as shown in Fig. 3.
The vector k has the coordinates k 1 along x axis and k 3 along z axis. We transform the axes so that x becomes parallel to the vector P 1 − P 3 . The coordinates of vector k expressed through the new coordinates k 1 and k 3 in the transformed system are The in (4) is also expressed through the coordinates (6). For a further calculation of (4) we also change the variables as follows so the tensor (4) can be rewritten as where the denominators (3) take the following form Note that z 1 and z 2 differ only in sign before q 0 . At the same time the only difference between z 3 and z 4 is the sign before q 1 . It allows us to shorten the calculations by introducing the following notations so that z 1 = z − 1 , z 2 = z + 1 , and z 3 = z − 3 , z 4 = z + 3 . Since the differential cross-section depends on the transmitted four-momentum t, we express the vectors P 1 and P 3 through Mandelstam variable t where P = | P 1 | = | P 3 |.
Applying one more change of variable and taking into account (11), we obtain the new expressions for (10): Finally, introducing the spherical coordinates (q, θ, φ): where 0 ≤ θ ≤ π and 0 ≤ φ < 2π, we can rewrite (8) as follows where As mentioned in the problem statement, we are going to use the Laplace's method to calculate the multidimensional integral (14).
The essential idea behind this method is that the integral whose integrand has a single maximum point in the integration domain, can be approximated nicely by the corresponding Gaussian integral. We have already seen above that either the first pair of denominators (associated with the horizontal lines of the diagram in Fig. 2a) or the second pair (two vertical lines) can be turned to zero simultaneously. Thus, we examine the first and the second pairs of denominators separately. We find the corresponding regions of the integration domain where either the first or the second pair of denominators turn into zero, and calculate the contribution of these regions to the integral.

III.1. The first pair of denominators
Let us consider the system of equations for the first two denominators z + 1 − iε and z − 1 − iε in (14). If we set the real parts of these denominators to zero, their product reduces to ε 2 , which provides the maximal contribution to the integrand in (14).
It is easy to verify that q 0 = 0, q = P satisfies the equation system (16). This solution describes the 2-dimensional sphere and thus provides some subset of the 4-dimensional integration domain where the absolute value of the first two denominators in (14) has a single minimum.
Next, let us apply Laplace's method to the first pair of denominators Considering the second order Taylor approximation of the exponent in (17) and changing the variable q = P + x in (15), we obtain The convenience of variable x = q − P is that together with q 0 they are the offset from the point (q 0 = 0, q = P ) where the real parts of denominators in the left-hand side of (18) take minimum value. Consequently, the absolute value of exponent in (18) has maximum at q 0 = 0, x = 0 .
Changing the variable q → x = q − P in the integral (14) and substituting (18), we get the approximation for the integral (14) where the same variable change q = P + x is also performed for f ab , z + 3 and z − 3 . We next turn to analyzing the integrand in (19). We denote it shortly by a ε q 0 , x, θ, φ . The exponent (18) in the integrand depends only on the integration variables q 0 and x. The absolute value of this exponent has maximum at q 0 = 0, x = 0 , and the width of its peak tends to zero as ε → 0. Since the other part of the integrand (19) containing the second pair of denominators does not have singularities at q 0 = 0, x = 0 , we can use such the exponential suppression of the integrand to cut the integration domain where q 0 cut > 0 and 0 < x cut < P . This idea is demonstrated in Fig. 4a. First, we found that the function a ε q 0 = 0, x = 0, θ, φ has maximum at (θ = 0, φ = 0). Then we considered a restriction of a ε q 0 , x, θ, φ to a ε q 0 , x, θ = 0, φ = 0 which we denoted by a ε q 0 , x and plotted its absolute value in the q 0 x plane.
In this way we consider the contribution to the original integral (14) supplied by the first pair of denominators z 1 and z 2 . Note that the actual value of this contribution does not depend on the selection of q 0 cut and x cut , which becomes clear as soon as we make the following change of variables It makes the integration (20) to take the form (22) and also cancels all the ε in denominators in (19). Since there is no ε in the denominator now, we can finally pass to the limit when ε → 0. At this point it is clear that the resulting integration limits will not depend on the particular selection of q 0 cut and x cut because both q 0 cut /ε and x cut /ε tend to infinity as ε → 0 (see Fig. 4b), while the width of the peak of |a (E, y) | remains unchanged. Now (19) can be rewritten as where f ab (θ, ϕ) now denotes f ab q 0 = 0, x = 0, θ, ϕ . The integrals over E and y in (23) are now reduced to Poisson integrals where e is the Euler's number (i.e. the base of the natural logarithm).
Combining (23) with (24) and substituting expressions for z ± 3 we obtain Finally, there remains only two-dimensional integral over θ and ϕ for each component t ab . Moreover, we can calculate the integral over ϕ for each pair of indices (a, b) analytically (these calculations are quite long but yet straightforward). So we end up with one-dimensional integral over θ which can be calculated numerically.

III.2. The second pair of denominators
Let us now consider the region of the integration domain in (4), where the divergence of the integrand arises due to the second pair of denominators If we let both the energy component k 0 of fourmomentum k and the length of its spatial part | k| tend to infinity simultaneously, the difference of their squares k 0 2 − | k| 2 may still be finite. In contrast, the sum k 0 2 +| k| 2 takes infinite values even when either k 0 or | k| tends to infinity. Note that the first pair of denominators contains k 0 in first power, which provides the decreasing of the integrand as k 0 tends to infinity. At the same time, in the center of mass system the energy component of P 1 − P 3 is equal to zero, and the denominators (26) do not contain the first power of k 0 , but only the difference k 0 2 − | k| 2 which, as we have already mentioned, may take finite values as k 0 → ∞ and | k| → ∞. In this case the whole integrand goes to zero due to the first pair of denominators containing the first power of k 0 . If we consider the component t 00 of the tensor (4), the numerator f ab k 0 , k = s − k 0 2 , where s is the Mandelstam invariant, tends to infinity as the denominator of the integrand. In this case the whole integral diverges. However, this divergence of (4) is only delusive. Let us show that. As mentioned above, if the expression contained the sum instead of the difference q 0 2 − | q| 2 , there would be no problem with the divergence at all. It is possible to turn this difference into a sum if we succeeded by transforming the integration contour in the complex plane from the real axis to the imaginary one, which can be done by the Wick rotation. The possibility of this transformation depends on the location of the integrand poles. Les us take a closer look at them.
First, we rewrite the expression (8) in the form Let us analyze the placement of quantities η and −η on the complex plane. It depends on the variables q 1 , q 2 , q 3 , so it varies across different subsets of the q 1 , q 2 , q 3 integration domain. Thus the major question here is which of the two conditions is met for the particular subset. Whether the first or the second condition is met, determines the distribution of the poles among the quadrants of the complex plane.
In case the second inequality (28) holds, each quadrant of the complex plane contains a pole of the integrand (27). It means that one cannot do Wick rotation, because one would cross the poles during rotation.
In the case of the first inequality, the poles are located in the second and fourth quadrants, so it is possible to rotate the integration path in the first and third quadrants. Thus the integration over q 0 can be transformed into the integration along the imaginary axis (Wick rotation), and the expression q 0 2 − | q| 2 transforms into − q 0 2 − | q| 2 as the integration variable changes q 0 → iq 0 . If one of these components tends to infinity, their sum also tends to infinity, which solves the problem. The integrand tends to zero as q 0 2 + | q| , which guarantees the convergence of the integral.
The possibility to apply Wick rotation arises in the region of high values of | q|, i.e. in that part of the integration domain where the integral can diverge. Therefore, selecting this region and applying corresponding Wick rotation, we obtain the convergent integral. The rest of the integration domain is finite and does not impact the convergence of the integral, since there are no singularities of the integrand in this region at ε = 0. This was actually the reason for calling the divergence of the integral "delusive" earlier in this section.

III.3. Calculation of the second contribution with the Wick rotation
Let us implement the idea explained in the previous section. We split the integral (8) into two integrals. The first integral over the finite region, where the Wick rota-tion cannot be applied, we denote as where c > 1. The value of c has no impact on the result of integration. However, if c = 1 and | q| < P , then the poles of integrand (27) lie on the integration path and the integral diverges, which imposes the lower bound on c.
The second contribution to the tensor t ab , i.e. the integral over the region where the Wick rotation can be applied, we denote as Here we also use the notations (10) for denominators with P 1 and P 3 expressed through t (11) Next we consider the system of equations for the second pair of denominators similar to what we did for the first pair of denominators. First, we set the real parts equal to zero and then find the second conditional maxima Let It is easily seen that q 1 = 0, q 0 = ±χ is the solution of (32). Note that q 0 = ±χ can take both positive and negative values, but we consider only positive solutions, since the integrand is an even function of q 0 . We can now proceed the calculation of the integral (29) analogously to what we did in Section III.1. First, we represent the second pair of denominators in exponential form In order to cancel the ε in denominators in (33), we make the following variable change (analogously to (21)) so after passing to the limit when ε → 0, the right-hand side of the expression (33) takes the form where ω ± = −2Eχ ± |t|x. Then we can rewrite the integral (29) as where D < = { q 2 , q 3 ∈ R 2 | q 2 2 + q 3 2 < c 2 P 2 } and f ab q 2 , q 3 denotes f ab q 0 = χ, q 1 = 0, q 2 , q 3 . According to Laplace's method [45], we transform the square roots in the denominators in (36) containing ω ± as 1 and take the second order Taylor approximation of the exponent at E = 0, x = 0 Similarly, for the other two terms in these denominators exp   i arccos Substituting (38) and (39) into (36), the integration with respect to E and x is now reduced to the Poisson integrals After these transformations we can rewrite the expression (36) as It is convenient to use the polar coordinates for further calculations q 2 = q cos (α) , Then the expression (41) takes the form As a result, we obtain the two-dimensional integral with respect to the variables q and α. This integral can be calculated numerically and allows us to calculate each component of the tensor t ab separately.
Let us get back to the expression (30). Applying Wick rotation, we change the variable q 0 → iq 0 in the integral (30), so it can be rewritten as where Note that ε in denominators can be set to zero now, because the integration is performed along the imaginary axis and we obtain the convergent integral. Rewriting the denominators in (44), we obtain where ξ 2 = q 0 2 + q 1 2 + q 2 2 . It can be seen that the product of two fractions in (45) is even function with regard to q 0 and q 2 and has the maximum at q 0 = 0, q 2 = 0. Increasing the absolute values of q 0 and q 2 results in increasing the denominators, i.e. in decreasing the whole fraction. Therefore we return back to (44) and separate the integration with respect to q 0 , q 2 and q 1 , q 3 . For the convenience and simplicity of the further calculations, we introduce the following polar coordinates The expression (44) then can be rewritten as where The expression (47) contains the integration with respect to q from cP to +∞. We can change the variable q → x = 1/q in order to make the integration limits finite. At this point the integrand in (47) becomes large enough if we substitute the expressions for (48) for z ± 1 and z ± 3 which depend on every integration variable x, α, q 0 , q 2 . However, introducing the following notations and changing the integration variables again q 0 → E = xq 0 and q 2 → y = xq 2 , we can finally rewrite (47) using notations (49) t > ab = i where , The further calculation of (50) depends on the values of indices a and b. We demonstrate the calculation workflow for the component t > 00 only, while the calculation of other components is completely similar. In this case (a = 0, b = 0) the function f ab in the integrand of (50) is equal to f 00 = E 2 + x 2 s , so the expression for t > 00 has the form The integration with respect to E and y can be reduced to the calculation of the Poisson integrals where W (x, α) = W 1 (x, α) + W 2 (x, α). Substituting (53) into (52) yields Finally, we obtain the two-dimensional integral with finite limits which can be calculated numerically. The other non zero components of t > ab have to be calculated in the same way.
After the components t < ab and t > ab get calculated, we can finally calculate their sum t ab = t < ab + t > ab and return to the analytical expression (1) for the diagram Fig. 2a.
So far, we have calculated the analytical expression for the simplest one-loop diagram shown in Fig. 2a. However, we should also consider all possible contractions of bispinor and glueball field operators. In the simplest case of the elastic scattering of two protons there is another diagram Fig. 2b added to the one-loop diagram Fig. 2a. The analytical expression for this diagram is pretty similar to (1) and the integral in this expression can be calculated in the same way as for t ab .

IV. RESULTS AND DISCUSSION
The analytic calculations presented in the previous sections allowed us to calculate the differential cross section of elastic proton-proton scattering dσ el /dt (t). This calculation includes the contributions from the tree-level (pole) (Fig. 1) and one-loop (Fig. 2) diagrams, as well as contributions from the diagrams with the P 3 and P 4 interchanged. The model used for calculation contains two parameters -M G -the mass of glueball, and G -the effective coupling of proton-glueball interaction. All the quantities were expressed in the units of the proton mass M P = 0.938 GeV. The obtained dependency dσ el /dt (t) is presented in Fig.5-8 at different energies √ s and different values of M G and G.
As can be seen in Fig.5, the obtained dependency is non-monotonic and qualitatively describes the first minimum of the experimental curve. The values of M and G were chosen in each case (energy) differently to better reproduce the experiment. However, the observed non-monotonic behavior of the theoretical dependency is  FIG. 6. The differential cross-section dσ el /dt (t) of elastic pp scattering at √ s = 30.5 GeV. The solid curve is the calculated dependency (MG = 0.13 and G = 2.2), the dots with error bars are experimental data [28]. preserved in all three cases even for the same values of M and G.
Although the results obtained are only in qualitative agreement with the experimental data, we expect that inclusion of the diagrams with the higher number of loops may help to achieve the quantitative description of the experiment.
In order to demonstrate the qualitative behavior of the dependency dσ el /dt (t) (including the second fall), we plotted it separately in a wider range of t up to 30 GeV 2 (Fig. 8).
We also performed the analogous calculations within another model φ 3 (phi-cubed) with the scalar field. However, within that model we have not obtained the effects of non-monotonic dσ el /dt (t) dependencies similar to the presented here. It suggests that the physical mechanisms  responsible for this non-monotonicity are associated with the spin effects. The spin flip effects observed in the protons scattering have been calculated within the Regge theory [51]. However, our results are different in that within the Regge approach this non-monotonicity arises due to an exchange of reggeons having different signatures, i.e. the reggeon contributions to the scattering amplitude with different signs. In [51] the spin effects provides only quantitative changes in the behavior of differential cross-section dependency dσ el /dt (t). In our approach it is the spin effects that are responsible for the appearance of non-monotonicity in the differential crosssection dependency dσ el /dt (t) .

V. SUMMARY
We demonstrated that the multi-particle fields approach can be used for a development of dynamic models which describe the scattering of multi-quark systems. Using this approach, we considered the elastic protonproton scattering and calculated the differential crosssection. As a result, we succeeded to obtain the qualitative agreement with the experiment.
Laplace's method allowed us to approximately calculate the loop diagrams. Unfortunately, the calculation is lengthy enough and requires a further improvement.
We found that the experimentally observed non-monotonicity of the differential cross-section dependency dσ el /dt (t) is caused in our model by the spin effects. Taking these effects into account in the tree-level and simplest loop diagrams allowed us to obtain the qualitative description of experiment. The obtained results suggest that inclusion of the more complex loop diagrams to the calculation will help achieve the quantitative description of experimental data.