Application of Methods of Computer Algebra for Doing Physical Sums

The methods of computer algebra are considered as modern trend in doing physical sums. It has been given a general classification of computer algebra systems. The system of analytical calculations and their application for solving various problems of the theoretical physics are stated in details.


Introduction
In solving the broad class of physical sums, especially the sums of statistic physics, physical kinetics, quantum chemistry, astrophysics, etc. by methods of perturbation theory appears the necessity of performing large amounts of monotonous lengthy analytical calculations.Formally, most of these problems can be solved with any preassigned precision.In fact, in the vast majority of cases, accounting of the third and higher orders of perturbation theory leads to a need of performing such a large amount of analytical calculations that it is not possible to do them manually for such a foreseeable time.In addition, by a high professional skill he possesses, is not immune from basic errors while performing such calculations, (no matter how skilled he is).Therefore, quite long time ago the methods and the systems of analytical solution of such tasks have begun to develop by using electronic computers symbolic and graphic information with the help of a computer is for a long period of time (approximately since 50th of the last century).Necessity in their development appeared with demand of solving bulky mathematical sums in different branches of science and technics.First of all they are the problems of the heavenly mechanics, quantum mechanics, static physics, quantum chemistry and other branches of science.The detailed classification of system, related to the 70th, is represented in the survey (Gerdt, V. P. et al., 1980).In the following years in connection with great increasing of possibilities of computer technics the working out the system of the analytical calculations and methods of CA obtained important development.
The aim of the present work is to give a short survey of the modern state of researchers in the development of systems of the analytical calculations (including the works of the authors) and their using in solving different sums of the theoretical physics.

Analytical Systems in the Sums of Mechanics and Heavenly Mechanics
In mechanics as in many branches of the applied researches more often the complex mathematical models of the investigated objects are considered.The researcher has to clash with necessity to operate by bulky analytical expressions when formating a model and then under its analyses.These difficulties may be so great that get narrow the model class, that may be compared and analysed at the first steps of research.
In connection (Pochtarenko, M. V., 1984) with it some more often met in practice the sums of mechanics are made algorithmic and created a package of the applied programs DYNAMICS of processing of symbol information for solving such sums.
For solving problem-oriented sums of mechanics such complexes of CA are used as NASRTAN, ANSYS, COS-MOS, MSC.ADAMS, APM WINMACHINE.They are convenient for the most complex calculations.The shortage of these packages is their orientation for discreet conception of incoming and outgoing information about the object of investigation.
In the packages of the object-oriented language of programming and modern surroundings of visual projecting program there are fitted systems of analytical calculation, realizing the symbol calculations.The representative of such packages is the program complex (PC) KIDIM developed at the department of the theoretical mechanics NTU HPI (Practical work in theoretical and analytical mechanics with using electronic computer, 2002, Andreev, Yu. M. et al., 2004;Lavinsky, D. V. et al., 2007).In the scientific work (Yermilov, A. N., 1989) the analytical solving of the Potts model is started in details including its simplex presentation.
In the heavenly mechanics (HM) one of the most important problem is the construction of the analytical theory of movement of heavenly bodies: planets, satellites, artificial Earth satellites (Sputniks), space stations, etc.For its construction it is necessary to solve the system of differential equations, describing moving of a body.Only in trivial cases the decision can be found in closed aspect.The most general approach to the construction of the analytical theory of heavenly bodies is the theory of perturbation.Considering that movement of bodies in HM, as a rule, almost periodical, the decision more often find in the row by trigonometrically functions, the arguments of which depend linearly from time, and coefficients of row are the polynomials by small parameters of the given tusk (sum).Such rows are called the Puason's rows.The high exactness, requiring in HM, leads to the necessity to take into account hundreds, thousands and even a dozen of thousands of members in Puason's rows and connected with them expansions.Therefore, as early as in 1985 (Davis, M. S., 1958) astronomers were interested in application of electronic computers for the analytical calculations and in a year the first program (Herget, P. & Musen, P., 1959) was published by IBM-650 for the analytical transformations of the Puason's rows.Later on quite a number of programs as well as CAB were worked out especially for calculations in HM (Gerdt, V. P. et al., 1980;Barton, D. & Fitch, J. P., 1972;Jeffrys, W. H., 1971).

The Analytical Systems in Quantum Mechanics and Quantum Chemistry
The overwhelming majority of the quantum mechanics problem, as a rule, has not exact analytical solution.The part of these sums is solved by either using simplified mathematical models which do not reflect' the whole completeness of the problem, or with using slowly coinciding rows of the theory of perturbation.With the development of CA it was possible to solve considerably more wide class of sums (Banshchikov, A. V. et al., 2008).
The most natural and exact method of researching the evolutional sums of the quantum mechanics is the solution of the non-stationary Schrodinger's equation because the calculation of the wave function gives the most possible description of the quantum system.
For example, for quantum chemistry the systems of the CA are the powerful tool of the standard analytical trans-

The System of Analytical Calculations in Statistical Mechanics, Physical Kinetics and Physics of Plasma
The use of perturbation theory allows analytical treatment of a number of problems in statistical physics, physical kinetics and the theory of non-ideal plasma.Computer algebra systems in their turn provide means of automation of the lengthy steps in deriving analytical expressions for the physical quantities one is interested in.
Here we will discuss one particular problem arising in the theory of non-ideal plasma physics, that is the calculation of electronic equilibrium distribution function of weakly non-ideal plasma.To effectively solve the problem at hand we also provide the details of computer algebra system specifically designed for the solution of this class of problems.
The expressions for equilibrium distribution functions of plasma species in weakly non-ideal (here, we restrict ourselves to the pair-collision approximation) plasma could be readily obtained as the solutions of Boltzmann kinetic equation.Of particular interest is electronic distribution function (Chapman, S. & Cowling, G., 1960).The approximate analytical solution of the corresponding Boltzmann equation could be conveniently obtained with the use of Chapman-Enskog method (Klimontovich, Yu. L., 1975;Polak, L. S., 1971;Ferziger, J. H. & Kaper, H. G., 1972;Lifshitz, E. M. & Pitaevskii, L. P., 1960).The non-linear integro-differential equation: is solved with the help of Chapman-Enskog's procedure, consisting in the expansion of the first order correction (correction of the order l/L 1, where l is the electronic mean free path and L is the characteristic length scale of the problem) to leading-order Maxwellian electron distribution function in a series of orthogonal Sonin-Lager's polynomials.The substitution of this ansatz in the above integro-differential equation and subsequent integration over the particle velocities leads to the infinite system of linear algebraic equations (Ferziger, J. H. & Kaper, H. G., 1972): Here, in (1,2) e k , m k , v k , f (v k ) stand for electrical charge, mass, velocity and distribution function for particles of sort k, g k j is the relative velocity of the colliding particles k and j, b stands for impact parameter, ϕ-azimuthal scattering angle and E is the intensity of electric field.Velocities marked by primes correspond to the velocities of particles after collision, T denotes plasma temperature, n e -electron concentration and δ kn stands for Kroneker's symbol.
In general, the solution of this infinite system of linear algebraic equations is obtained by truncating the above series and taking into account only a number of first terms of the expansion.Careful examination of the integrals entering expressions for the coefficients of our infinite system of linear algebraic equations shows, that one needs to know the transport cross-sections (discussed below) at high orders.But the most difficult part consists in calculation of multiple bracket integrals (1).Nevertheless, this truncated system of equations could be solved analytically -the only problem is that the solution procedure should be definitely automated.The use of general-purpose computer algebra systems for this task (we are having in mind to account for as more terms in the above mentioned expansion as possible) is not feasible, so here we will present the details of custom computer algebra system written in Visual Basic, specifically designed to solve this kind of problems.

The Analytical Calculation of the Transport Cross-sections
The description of collision processes is provided by corresponding collision cross-sections.According to the classical mechanics the differential cross-section for electron-electron (e − e) and electron-ion (e − i) scattering is given by: where b-is the impact parameter, which in the case of Coulomb and short-range Coulomb (here we are accounting for the Debye screening) potentials (srCP) is given by: (5) Here ξ = r cp /b 0 ± 1, b 0 = e 2 /μv 2 , χ-scattering angle, μ-reduced mass, v-the relative velocity of particles and b 0 is the impact parameter, corresponding to the scattering angle equal to 90 degrees.The minus sign in the expressions for ξ corresponds to the scattering of opposite charges, while the plus sign to the scattering of similar charges.Differentiation of expressions ( 4) and ( 5) with respect to scattering angle χ and substitution of the resulting expressions in (3) leads to -for Coulomb potential and -for short-ranged Coulomb potential.
In the limit ξ 1, the main contribution into cross-section comes from the small angle scattering and Equation (7) transforms into Rutherford formula (Kudryavtsev, E. V. & Ionov, V. P., 1961).The transport cross-section at n's order is given by the following expression (Vyzhol, Yu. A. & Fedoseeva, E. V., 1998): The order of the transport cross-section used is influenced by the mass correlation of the colliding particles.Momentum transmission under collision of light and heavy particles is determined by the transport section of the first order, while for the collision of particles with the same mass, such as electron-electron scattering, one needs to account for high-order cross-sections.Integrating (8) with n = 1 in the case of (e − i) and (e − a) collisions we get: -for Coulomb potential and -for short-ranged Coulomb potential.
The difficulties appear in the case of the transport cross-sections for the colliding particles with the same mass (electron-electron (e − e) collisions).In this case, the transport cross-section of the first order σ 1 t ee vanishes, as the average momentum transfer is zero here.Non-vanishing contributions are provided by transport cross-sections of even orders n = 2, 4, 6, ....The analytical expressions for these transport cross-sections could be easily obtained with the help of any general-purpose computer algebra system.Let us illustrate how this could be done with MathCad.First, one needs to substitute Equations ( 6) or (7) into Equation (8).As a result one obtains the following expressions for the transport cross-sections in the case of Coulomb and short-range Coulomb potentials: Now, with the use of MathCad functions expand and f actor, the transport cross-sections of 2 nd and 4 th orders in the case of Coulomb potential take the form: Similar, the expressions for transport cross-sections of 6 th and 8 th orders in the case of short-range Coulomb potential are given by: Both numerical and analytical solutions of the kinetic equation, with the exception of so called τ-approximation, when the collision integral is expressed as S t f e = f e − f 0 e τ , turn out to be difficult to solve problems.On one hand we are always dealing with the functions of several variables, depending from a number of macroscopic parameters, such as plasma temperature, composition and so on.On the other hand, approximate analytical solutions for distribution functions always involve a large number of routine analytical transformations for expressions having dozens and hundred thousands of terms.It is precisely this problem of huge number of terms, which one encounters in the calculation of bracket integrals.To deal with this problem we first worked out the complete algorithm for the whole process of obtaining approximate analytical solution to the Boltzmann kinetic equation, starting from original bracket integrals and ending with the expressions, which could be calculated numerically, and second-we designed a custom computer algebra system implementing it.Proceeding this way allows us considerably reduce the dimension of integration space, such that the remaining integrals could be easily calculated numerically.

Analytical Calculation of Bracket Integrals
In procedure (2), the most laborious task is the calculation of the bracket integrals To integrate with respect to the angles in (17), it is convenient to change from the laboratory frame of reference to the frame of reference related to the center of mass of the colliding particles; then, the bracket integral takes the form where G and g are the velocity of the center of mass and the relative velocity of the colliding particles and S 3/2 n (x) is the Sonin polynomial of degree n and order 3/2.The further calculation of integrals (18) in the general form for an arbitrary ratio of the masses of the colliding particles is a tedious task.For that reason, we first perform the calculations for two limiting cases when the masses of the colliding particles are very different from each other and when the colliding particles have the same mass.Then, we generalize the results for the case of the arbitrary ratio of the masses.

Calculation of the Collision Integral for Light and Heavy Particles
For the specific cases of electron atomic (e − a) and electron-ionic (e − i) collisions, the calculation is considerably simpler because the ratio of the masses is small: 1. Taking into account the fact that weak external perturbations do not distort the Maxwell distribution of the ionic and atomic components, we obtain from ( 18), upon integrating with respect to the angles θ , ϕ and the absolute values of the center of mass velocity G, one can show that the bracket collision integral for the light and heavy particles has the general form The Sonin polynomial where r is a real number and s is a nonnegative integer.For example, The condition for orthogonality of these polynomials for the given index r and various indexes ss is It is seen from ( 19) that the coefficients L e j rn form a symmetric matrix; that is L e j rn = L e j nr .In the general case, integrals in (19) depend on the temperature and concentration of particles as parameters.

Calculation of the Collision Integrals for Particles of the Same Mass
Consider the bracket integrals in the case when the colliding particles have the same mass using electron-electron collisions as an example.Since the momentum exchanged between the particles is zero in this case, the terms containing cross sections of odd orders 2k − 1 in the collision integral vanish.The bracket integral characterizing the collision of identical particles of other kinds has exactly the same form; only the mass of the electron and the transport cross sections for electron-electron scattering of various orders are replaced by the mass and the transport cross sections for scattering of particles of the corresponding kind.The coefficients L ee rn also form a symmetric matrix; in addition, L ee 0n = L ee n0 , (n = 0, 1, 2, ...), which is also characteristic of the bracket integrals of any two particles of the same mass because these integrals are obtained in the general form without taking into account the specific interaction potential.We calculate the coefficients L ee rn using L ee 11 as an example.Substituting the Sonin polynomial of the first degree S 3/2 1 (x) = 5/2 − x, the expressions for the cross sections for electron-electron collisions (see Vyzhol, Yu. A. et al., 1997), and taking into account the equality of masses of the colliding particles, we obtain Making the change of variables in ( 14) and performing elementary algebraic transformations of the integrand, we obtain Here, η is the angle between the velocities g and G, g is the relative velocity of the colliding particles after the collision, and G is the velocity of the center of mass of the two-particle system.Since the collisions are elastic, we have ||g | = |g||, where g is the relative velocity of the particles before the collision.For the scalar product, we have g G = gGcosη, where cosη = cosχcosθ + sinχsinθcosϕ.Integrating with respect to ϕ and θ and transforming the integrand, we obtain The last integral in ( 25) is easily calculated, and we finally obtain where σ 2 t ee (y 2 ) is the second-order transport cross section for electron-electron scattering.Thus, using the coefficient L ee 11 as an example, we see that the integrals for electron-electron collisions can be reduced to a double integral.
For electron-electron interactions, no general formula of form ( 19) can be obtained; therefore, one has to calculate each bracket integral individually.As the indexes r and n in the bracket integral L ee rn grow, the number of required analytical transformations grows as a factorial.For r, n ≤ 4, the results of such calculations are presented in (Kudryavtsev, E. V. & Ionov, V. P., 1961).When r, n > 4, it is impossible to carry out these calculations manually.For that purpose, we developed a system of analytical computations; the main characteristics of this system are discussed in the next section.

Calculation of the Collision Integrals for Particles of Arbitrary Masses
Consider the general case when the ratio of the masses of the colliding particles is arbitrary.Assume that two particles of the kinds k and j of masses m k and m j , respectively, collide.Define the quantity μ = m j /m k .Then, from the relations M j = m j /(m e + m j ); M e = m e /(m e + m j ), we obtain The quantity μ varies in the range [1, ∞); the lower limit of this range corresponds to the collision of particles of the same mass; the upper limit corresponds to the collision of light particles with heavy ones.Using the results of paragraphs 4.3 and 4.4 and applying method of mathematical induction we will get the general expression of the bracket integral of arbitrary order L k j rn ; the calculation of this integral is reduced to the cases considered in the preceding sections:

Analytical Calculation of Bracket Integrals on a Computer
Now, we consider an implementation of a computer system for analytical calculation of bracket integrals for colliding particles of the same mass.To calculate such integrals, one has to calculate the expressions (29) where It is seen from ( 21) that the initial, intermediate, and final results can be written using polynomials of the form where K a i is a constant coefficient, a a i j is exponent of a factor in the ith term of the polynomial A, and a is the degree of this polynomial.
The program was written in Visual Basic 6.0.To store each term in computer memory, eight bytes for the coefficient K and nine bytes for the exponents are needed.In addition, to speed up the search for similar terms, the control sum (four bytes) α a i = 8 j=0 α a i j d j , where d ≥ α a i j for i, j ∈ R was calculated and stored in memory for each term.In the course of intermediate transformations, the number of terms in the polynomials was as high as several hundreds of thousands.The solution of this problem using the integer representation of numerators and denominators of the coefficients in (31) in computer memory and calculations based on rational fractions yields exact results for the elements up to L ee 2020 inclusive.This algorithm was implemented in (Vyzhol, Yu. A. & Fedoseeva, E. V., 1998) and used for the calculation of the electron kinetic coefficients of non-ideal plasma in (Vyzhol, Yu. A. et al., 1997;Zaika, E. V. et al., 2000;Vyzhol, Yu. A. et al., 1999).In this algorithm, the greater part of the execution time is spent on finding the greatest common divisor for reducing fractions.The calculation of bracket integrals of higher orders than L ee 2020 is needed because the Chapman-Enskog process converges nonuniformly for various types of interaction potentials (transport sections) of colliding particles.The use of floating point representation for the coefficients considerably (by several orders of magnitude) reduces the execution time but results in roundoff error accumulation, which reduces the accuracy of the coefficients by four to five orders (compared with machine precision) (Vyzhol, Yu. A. et al., 2011).Upon integration with respect to each variable, similar terms appeared that canceled one another.When the coefficients were represented by floating point numbers, a considerable roundoff error was accumulated.For that reason, to estimate the probability that a coefficient is equal to zero, the computational error was simultaneously estimated.If the error became of the same order as the coefficient itself, the corresponding term of the polynomial was eliminated.
To estimate the accuracy of computations, the absolute error of the coefficient K of the term was stored in memory (four bytes); this error was calculated using the well-known formulas Polynomials were stored in dynamic arrays, and each element of such an array occupied 25 bytes of memory.

Implementation of Elementary Operations in the Program
To carry out the required transformations, the following elementary operations on polynomials were defined: addition, multiplication, and integration with respect to the variables ϕ, θ, and x.To reduce the length of the polynomials and, correspondingly, the amount of memory, similar terms were collected after each elementary operation.
Below, we describe the algorithms using a pseudocode, which is similar to popular programming languages (C or Pascal).The indentation indicates the nesting level of the corresponding operator.Every line includes comments after the symbol .Sometimes, the operation of the algorithm is described informally when this is clearer.In addition, some technical details (for example, error handling and input-output operations) that are necessary in real life programs but can obscure the essence are omitted (see Vyzhol, Yu. A. et al., 1997).
Polynomials are added (C = A + B) using the following algorithm that includes the operation of collecting similar terms.
C ← A the polynomial A is copied to the polynomial C for i ← 1 to Lb Lb is the length of the polynomial B find ← false find is the flag indicated whether there is a similar term in the result if the control sums of the terms are equal to each other Vol.4, No. 3; 2012 then find ← true the result contains a similar term After running this program, we obtain the resulting polynomial in the form This is a part of the expression obtained at the intermediate phase of calculations after collecting similar terms copied from the computer screen (the ellipses symbols at the beginning and the end of the expression denote the set of all the preceding and succeeding terms).Depending on the order of the bracket integrals, such expressions cover from several tens to several thousands of screen pages.
Polynomials are multiplied (C = AB) using the following algorithm that includes the operation of collecting similar terms.
for i ← 1 to La La is the length of the polynomial A The coefficient of the Sonin polynomials were calculated using the formula where m!!! is the product of all odd integers that are less than or equal to m.
Calculation of the Sonin polynomials (C = S 3/2 n (A)) of the polynomial A and the calculation of C was performed using the following algorithm.
Integration with respect to ϕ was carried out using the formulas 2π 0 where m!! is the product of all even integers that are less than or equal to m.To integrate the polynomial A and form the polynomial C, the following algorithm was used.
for i ← 1 to La for each term of the polynomial Integration with respect to x was carried out for even powers (the terms containing odd powers were canceled out) using the formula Integration with respect to θ was carried out using the formulas The integration algorithm with respect to and θ is analogous to the integration algorithm with respect to ϕ.The integrals calculated using this algorithm have the form Here, the ellipses symbol also denotes the set of all the preceding or all the succeeding terms in the expression.
The factors sin m χ were eliminated from the polynomials using the formula sinχ (odd powers of sin m χcos n χ = sin m−2 χcos n χ − sin m−2 χcos n+2 χ did not occur in the calculations).To form the polynomial C from the polynomial A, the following algorithm was used.

Repeat
FlagSin ← False FlagSin indicates whether there are terms including the factor sinχ The terms containing 1 − cos m χ were distinguished in the polynomials using the formula cos m χ = (1 − cos m χ) + 1; an algorithm similar to that used for eliminating sin m χ was used.

Elimination of Zero Terms
In the course of calculations, many terms cancelled out, and, for small values of the indexes (r, n < 3), the result contained terms with zero coefficients.Such terms were eliminated from the result.For large values of the indexes r and n, the result contained terms with the coefficients that did not include the factors 1 − cos m χ.This indicates that such terms appeared due to roundoff errors.A small absolute value of the coefficient cannot be the reason for deleting the corresponding term from the resulting polynomial because the coefficients of high-order Sonin polynomials can be small.For that reason, the absolute error was estimated after each elementary operation on the coefficients of the polynomial, and the term was removed from the result only if the relative error was ≈ 1.
In addition, the error estimates showed that the accuracy of the coefficient computation can be considerably reduced in the course of computations (up to five orders of magnitude compared with the machine precision).In this paper, we give only correct significant digits (with account for roundoff errors) in the results.

Algorithm Underlying the Program
The algorithm for calculating the expression It is easy to verify that expression (34) can be represented in the form Calculations by formula (35) do not require the calculation of the polynomial W and considerable effort for finding similar terms in the resulting polynomial because the polynomial F rn (y, χ) is fairly short.However, in this case, the multiplication, integration with respect to three variables, and collection of similar terms must be performed within a subroutine that must replace the block of subroutines enclosed in a frame in algorithm (34).The algorithm underlying this subroutine is as follows.
for i← 1 to Lv Lv is the length of the polynomial V

Validation of Operations
The validity of performing operations was checked by comparing the results of numerical and analytical computations.The integration results for the coefficients up to the fourth order were compared with the coefficients calculated in (Kudryavtsev, E. V. & Ionov, V. P., 1961); for the coefficients of the fifth to seventh orders, they were compared with the coefficients calculated in (Vyzhol, Yu. A. & Fedoseeva, E. V., 1998;Vyzhol, Yu. A. et al., 1997;Zaika, E. V. et al., 2000;Vyzhol, Yu. A. et al., 1999;Vyzhol, Yu. A. et al., 2009;Vyzhol, Yu. A. et al., 2011).The comparison was performed using a graphical representation of the formulas in the conventional form.
The layout of the start-up screen of the program is shown in the figure.Depending on the value of the degree and order of the collision integral, the general result or its individual terms are obtained.

Calculation Results
Using the program described above, we calculated the bracket integral up to the 20th order.The final result of the analytical computations was obtained in the form Double integral (36) can be used for the numerical integration with respect to y and χ for the computation of the electron kinetic plasma coefficients.
As a result, an array of numerical values of the coefficients K i k in the integrand in (36) up to the 20th order inclusive was obtained.This array is huge, and we presented in the Appendix of (Vyzhol, Yu. A. et al., 2011) only the coefficients for the ninth approximation (in this table, y is the degree and D is the exponent of cosχ).

Applications of Asymptotic Expansions to the Problems of Quantum Field Theory (QFT)
To illustrate the power of modern CAS applied to the problems of quantum field theory we chose three particular problems related to the calculation of high order radiative corrections to the different observables in Standard Model (SM) and Minimal Supersymmetric Standard Model (MSSM) solved with the use of large-mass asymptotic expansions.Thus, after a general description of asymptotic expansion procedure, we present results of calculations of two-loop bosonic electroweak corrections to r, two-loop O(α 2 s ) MSSM corrections to the pole masses of heavy quarks and light-by-light scattering contribution to the anomalous magnetic moment of the muon.

Asymptotic Expansion
One of the approaches to the computation of multi-loop diagrams (at least in certain kinematical limits) is based on asymptotic expansion (Note 1).An asymptotic expansion may be viewed as a generalization of a Taylor expansion.In both cases we have an expansion in powers of a small quantity.However, in case of the asymptotic expansion contrary to the case of Taylor expansion the corresponding coefficients are not constant but contain non-analytical functions of the small parameter.As an example let us consider the function f (x) = Li 2 (1 − x) for which the Taylor expansion for x → 0 does not exist beyond the leading order.Nevertheless there is an asymptotic expansion which reads with the non-analytic (for x → 0) function ln x in the coefficient of the linear term in x.
In the case of Feynman diagrams the situation is similar.There are systematic procedures for the cases of a large external momentum (large-momentum expansion) and a large internal mass (large-mass expansion) (Smirnov, V. A., 1991;Smirnov, V. A., 1995).Both procedures apply to problems which can be formulated in Euclidean space.
In contrast to that there are also phenomena which are tightly connected to the Minkowskian space-time.Also in such cases there are rules for an asymptotic expansion.In particular, two-loop on-shell two-point diagrams (Smirnov, V. A., 1997a;Czarnecki, A. & Smirnov, V. A., 1997) and two-loop vertex diagrams in the Sudakov limit (Smirnov, V. A., 1997b) have been considered.Also the threshold expansion (Beneke, M. & Smirnov, V. A., 1998) belongs to this class of phenomena.
It should be noted, that both large-momentum and large-mass expansion procedures could be treated similarly.So, here we only present the general formulae for the case of large external momenta-the transition to the large-mass expansion procedure is straightforward.The prescription for the large-momentum procedure is summarized by the following formula (Note 2): Here, Γ is the Feynman diagram under consideration, {Q} ({m, q}) is the collection of the large (small) parameters, and the sum goes over all subgraphs γ of Γ with masses m γ and external momenta q γ (the conditions, to which subgraphs should satisfy, will be given below).T {q,m} is an operator performing a Taylor expansion in {q, m} before any integration is carried out.The notation Γ/γ T {q,m} γ indicates that the subgraph γ of Γ is replaced by its Taylor expansion which should be performed in all masses and external momenta of γ that do not belong to the set {Q}.
In particular, also those external momenta of γ that appear to be integration momenta in Γ have to be considered as small.Only after the Taylor expansions have been carried out, the loop integrations are performed.In general, the set {γ} is referred as hard subgraphs or simply subgraphs and {Γ/γ} as co-subgraphs.
The conditions for the subgraphs γ are different for the large-mass and large-momentum procedures (Note 3).For the large-momentum procedure, γ must • contain all vertices where a large momentum enters or leaves the graph • be one-particle irreducible after identifying these vertices.
From these requirements it is clear that the hard subgraphs after Taylor expansion become massless integrals where the scales are given by the large momenta.In the simplest case of one large momentum one ends up with propagator-type integrals.The co-subgraph, on the other hand, may still contain small external momenta and masses.However, the resulting integrals are typically much simpler than the original one.
In the case of large-mass expansion procedure, γ has to • contain all the propagators carrying a large mass • be one-particle irreducible in its connected parts after contracting the heavy lines.
Here, the hard subgraphs reduce to tadpole integrals with the large masses setting the scales.The co-subgraphs are again simpler to evaluate than the initial diagram.
As an example of application of large mass expansion procedure let us consider one-loop diagram depicted in Figure 1 and given by the following expression: The large-mass expansion applied to F(M 2 , m 2 , q 2 , d) could be written as It should be noted, that the main simplification, which is common to all kinds of asymptotic expansions, comes from the fact that the expansions in the small parameters are done before any momentum integration is performed.The proof that this leads to correct results is based on the so-called strategy of regions (Smirnov, V. A., 1999).There different regions of each loop momentum are selected and in each of them Taylor expansions with respect to the small parameters are performed.In the limits of the large-mass and large-momentum procedures an interpretation of the different regions in terms of subgraphs and cosub-graphs is possible (see above).This is, however, different in the case of the threshold expansion (Beneke, M. & Smirnov, V. A., 1998) where a graphical representation becomes much less transparent.However, the application of the strategy of regions (Smirnov, V. A., 1999) leads to correct results.

Two Loop Electroweak Corrections to r
Let us consider an application of asymptotic expansions to the calculation of two-loop bosonic electroweak corrections to r.
Fermi constant is defined as the coupling constant in the low energy four fermion effective Lagrangian describing the decay of the muon where e and μ are electron and muon fields, ν e and ν μ are the corresponding neutrinos and G F is the Fermi constant.
From the Lagrangian (42) one gets the following value for the muon lifetime where the factor Δq describes all the quantum corrections in the low energy effective theory (i.e.QED corrections).At one-loop order these corrections have been computed a long time ago (Berman, S. M., 1958;Kinoshita, T. & Sirlin, A., 1959).Recently also the two-loop result for Δq has been obtained (van Ritbergen, T. & Stuart, R. G., 1999;van Ritbergen, T. & Stuart, R. G., 2000).Taking it into account, the error of G F is nowadays dominated by the experimental error of τ μ measurement (Note 5).
In order to relate G F to M W one can use the matching condition between effective theory (42) and the Standard Model, which requires that the value of τ μ does not depend on whether it is evaluated in the Fermi theory or in the full Standard Model.This brings us to the relation where e is the electric charge, s W = sin θ W is the SM mixing angle and M W is the W-boson mass.The factor e 2 /8s 2 W m 2 W corresponds to the matching at the tree level, while the quantity Δr = Δr 1−loop + Δr 2−loop + . . . is defined to absorb quantum corrections coming from the matching procedure at higher loop orders.At the one-loop level Δr was first evaluated in (Sirlin, A., 1980).Further improvement of Δr was achieved in (Degrassi, G. et al., 1996;Degrassi, G. et al., 1997) where the leading and subleading large t-quark corrections were computed at O(α 2 ) order.Other two-loop fermionic corrections were considered in a series of papers (Freitas, A. et al., 2000a;Freitas, A. et al., 2000b;Freitas, A. et al., 2002).The calculation at the two loop level has been completed recently with the evaluation of the electroweak bosonic corrections (Awramik, M. & Czakon, M., 2002;Awramik, M. & Czakon, M., 2003a;Onishchenko, A. & Veretin, O., 2003a;Onishchenko, A. & Veretin, O., 2003b;Awramik, M. et al., 2003).Subsequently, the fermionic part, which was first obtained in (Freitas, A. et al., 2000;Freitas, A. et al., 2002) has been recalculated (Awramik, M. & Czakon, M., 2003b).A difference between the two results has been found and the original result (Freitas, A. et al., 2000;Freitas, A. et al., 2002) corrected.There is now agreement on all parts of the two loop contributions.
Let us say a few words about factorization procedure.In order to separate the short and long distance dynamics at the level of a single Feynmans diagram we use the large mass expansion procedure (Tkachov, F. V., 1983;Chetyrkin, K. G., 1988;Smirnov, V. A., 1990).Namely, given the graf F, then asymptotically where the sum runs over all "hard" subgraphs H of the diagram F; S is a "soft" subgraph obtained from F by shrinking H to a point and T stands for the Taylor expansion (before integration!) of H with respect to all "soft" parameters.The exact rules for construction of hard subgraphs are discussed in details in (Tkachov, F. V., 1983;Chetyrkin, K. G., 1988;Smirnov, V. A., 1990).
At the level of Feynman amplitude the situation is more involved, since each diagram has its own subgraphs.However it is possible to rearrange the action of asymptotic expansion (45) in the sum of all diagrams such that it exponentiates.Such rearrangement is based on the combinatorial properties of R-operation and the special structure of the expansion (45).The rigorous prove of factorization theorem can be found in (Gorishnii, S. G., 1989).

<Figure 3>
In order to demonstrate the idea we consider a simple example.In Figure 3 the contibution of ladder topologies to muon decay amplitude A SM is depicted, where A SM is given by One can notice however that this sum can be writen as a product of two factors.One corresponds to short distance Wilson coefficient function (G F ) and another to long distance matrix element.All other topologies can be considered similary.As a result of this procedure, the evaluation of G F is reduced to computaion of bubble diagrams only.
To perform the calculation we used the methods of asymptotic expansion (Tkachov, F. V., 1983;Chetyrkin, K. G., 1988;Smirnov, V. A., 1990) in two different regimes: 1) large Higgs mass expansion and 2) expansion in difference of masses M 2 H − M 2 Z .In addition we also considered s 2 W = sin 2 θ W as a small parameter in the spirit of (Jegerlehner, F., 2002).
• G F is nothing but the Wilson coefficient function of the corresponding operator in the low energy four fermion effective theory.Therefore it depends only on the "hard" physics and is insensitive to low energy field modes in loops.The simplest way to compute such a quantity is the use of the factorization theorem (if it exists), which allows one to separate "hard" and "soft" physics (Note 6).This theorem for the construction of low energy effective theories is known for a long time (see e.g.Gorishnii, S. G., 1989) and is based on the euclidean structure of large mass expansion.
• Using the factorization theorem the problem is reduced to the evaluation of two-loop bubble diagrams.
• All the calculations have been performed in a general R ξ gauge with three arbitrary gauge parameters.By explicit calculation we checked that contribution to Δr is gauge invariant.
• In order to have explicit gauge invariant expressions also at the level of bare quantities we explicitely include tadpole diagrams.We have checked that applying this prescription all the bare parameters and all the counterterms for the masses and the coupling constant are explicitely gauge invariant.
• The renormalization was performed in two different schemes: MS and on-shell scheme.For the transformation from one scheme to another one needs the results of (Jegerlehner, F., 2002) where the transformation formulae for the gauge boson masses were given to two loops.Vol. 4, No. 3;2012 To obtain the two-loop bosonic contribution to r a number of two-loop diagrams of the order of 10 4 should be evaluated.In our calculation the computer algebra system FORM (Vermaseren, J. A. M., 1991) was used.In order to obtain the FORM readable input we make use of the input generator DIANA (Tentyukov, M. & Fleisher, J., 2000).
The results for the on-shell renormalised Δr (2)  bos were given in a twofold expansion, in the large Higgs boson mass and in the mass difference between the W and the Z boson.The leading behaviour both in the Higgs mass and in the sine of the Weinberg angle has been factorized out (Note 7).
Note that the leading term in the Higgs boson mass can be resumed in sin 2 θ W to give the behavior The expansion coefficients read (these coefficients expanded to the order O(z 5 H ) could be found in (Awramik, M. & Czakon, M., 2002;Awramik, M. & Czakon, M., 2003;Onishchenko, A. & Veretin, O., 2003;Awramik, M. et al., 2003) ) MSSM Corrections to the Pole Masses of Heavy Quarks The particle pole mass is defined by the singularity of the corresponding two-point function.As explicit perturbative calculations showed, the pole mass of a quark is an infrared finite and gauge invariant quantity.For this reason it is considered as a physically meaningful quantity (Tarrach, R., 1981).Here we employ the definition of the pole quark mass, where it equals the value of p, at which the inverse quark propagator turns to zero.
The pole mass M of a quark is defined as a formal solution for p (in the Minkowski metric) at which the reciprocal of the connected full propagator equals zero where Σ( p, m) = mΣ 1 (p 2 , m) + ( p − m)Σ 2 (p 2 , m) is the one-particle-irreducible two-point function, m may stand for the bare or renormalized mass, m bare or m ren , depending on the prescription used in evaluating Σ.The solution to equation ( 50) is sought order by order in perturbation theory.To two loops we have where Σ (L) is the L-loop contribution to Σ, and the prime denotes the derivative with respect to the first argument.
In what follows, we will be interested in the relation between the pole quark mass M and the running DR mass m computed in Minimal Supersymmetric Standard Model (MSSM) up to the O(α 2 s ) order (Bednyakov, A. et al., 2003).Technically, to solve this problem, we need to evaluate ≈ 200 two-loop propagator type diagrams involving many different mass scales (see Figure 4).In general, it is quite a complex problem to be solved exactly.However, in our case all mass scales can be divided into two groups denoted by m soft and m hard , such that scales from m soft are much less than each of the scales from the m hard group Here m g is the gluino mass, m q stands for different squark masses.We also explicitly denoted which mass scales belong to soft or hard groups in the cases of b-and t-quarks.
Given this mass hierarchy, one can employ the large mass expansion procedure to reduce evaluation of multi scale two-loop integrals to the calculation of single scale on-shell two-loop integrals, two-loop tadpole integrals with two scales and products of one-loop on-shell integrals and one-loop tadpole diagrams.To compute two-loop single scale on-shell integrals, we made use of the ONSHELL2 package (Fleischer, J. & Kalmykov, M. Y., 2000) .For evaluation of two-loop tadpole diagrams with two scales, the recurrence relations of (Davydychev, A. I. & Tausk, J. B., 1993) were used.Calculation of one-loop self-energies, including their derivatives with respect to momentum, is an easy task and we will not make further comments on it.

<Figure 4>
The two-loop QCD results in the DR-scheme for the b-quark (with four light quarks)  , 1997).
We find very a simple answer to one particular limit.Namely, in the case, when mt where a q = A q − μ cot β Here A q -soft trilinear coupling in MSSM lagrangian, μ-Higgs mass parameter and tan β = v 2 v 1 -ratio of vacuum expectation values of two MSSM Higgs fields.
Note, that we put equal to each other the left and right squark masses m 2 qR = m 2 qL but not its physical mass eigenstates m 2 q1 = m 2 q2 .The complete formulae (Note 8) are too large to be presented here and could be found in (Bednyakov, A. et al., 2003).

Light-by-light Scattering and the Anomalous Magnetic Moment of the Muon
In the standard model (SM) the theoretical value of the muon anomalous magnetic moment a μ = (g μ − 2)/2 is given by a sum of three contributions: a μ (SM) = a μ (QED) + a μ (weak) + a μ (had).The contribution from QED, a μ (QED), which includes those from virtual leptons, and the one from weak interactions, a μ (weak), can be uniquely evaluated in perturbation theory with the results: a μ (QED) = 11658470.57(0.29)× 10 −10 (Mohr, P. J. & Taylor, B. N., 2000) and a μ (weak) = 15.1(0.4)× 10 −10 (Czarnecki, A. & Marciano, W. J., 2001).The hadronic piece a μ (had), however, is sensitive to long distance effects and cannot be evaluated in a perturbative framework.An important NLO contribution to a μ (had), which will be considered as example here, originates from light-by-light scattering.One important part of light-by-light-scattering amplitudes comes from neutral, low mass intermediate resonances, dominantly π 0 and, less important, η, another from charged pion loops.
To evaluate the contribution of charged pion loops we consider the pion as an elementary pointlike particle and use the Lagrangian to describe its interaction with photons.This is well justified since the integral is convergent in the high energy region, and the internal structure of the pion is not yet resolved for momenta of order m μ or m π .In view of the low mass of the muon and pion, as compared to the characteristic scale of the pion form factor, m 2 ρ ≈ 0.6 GeV 2 , pions can be treated as pointlike.
The 21 diagrams for the light-by-light contribution as derived from scalar QED are displayed in Figure 5 (permutations of legs omitted).Numerical results for this contribution were already reported in refs.(Kinoshita, T. et al., 1985;Hayakawa, M. et al., 1996).Analytical expressions, however, were obtained for the first time in (Kuhn, J. H. et al., 2003).

<Figure 5>
To compute the contributions of different diagrams we first project onto the relevant form factor of the anomalous magnetic moment (Aldins, J. et al., 1970;Kukhto, T. V. et al., 1992).Subsequently we employ the well-known method of asymptotic expansion in the small mass ratio m/M (Tkachov, F. V., 1983;Chetyrkin, K. G., 1988;Smirnov, V. A., 1990).As an example let us consider the large mass expansion of the generic diagram depicted in Figure 6.To obtain the asymptotic expansion for a given diagram one has to compute the sum of different a μ (γγ) . (57) The expansion in scalar QED, relevant for the charged pion contribution, has the form  These results where obtained in a general covariant gauge, thus providing additional checks at different steps of the calculation.

Conclusion
The using of the system of CA in physics does not exhaust by the examined examples here, it is perspective and effective sphere of investigation.In this work only a short survey of modern state of investigation in the sphere of working out the system of the analytical calculations given and their using to solution of different sums of mechanics, physical kinetics and physics of plasma, in quantum theory of field.
The system of analytical calculations on electronic computers in the kinetics of plasma is determined in details, the realization of the program of the analytical transformations of the bracket integrals is shown, the analytical calculation of the transport sections on electronic computers and the system of analytical calculations, worked out for calculation of Feynman's diagram in the quantum theory of field.
Using the system of CA in the scientific research lets to automatize and, consequently, to accelerate the analytical calculations.Vol. 4, No. 3;2012 Notes Note 1.Here, for the presentation of asymptotic expansion procedure we are closely following (Steinhauser, M., 2002).
Note 2. In the case of the large-mass expansion procedure one essentially has to replace Q by the large mass M.
Note 3. Actually they are very similar and it is certainly possible to merge them into one condition using a more abstract language.
Note 5.The present value is G F = 1.16637(1) × 10 −5 GeV −2 and possible future experiments could reduce the error by an order of magnitude.
Note 6.In our case by "soft" we understand energies and momenta of order of muon mass m μ , while "hard" means scales m μ .
Note 8.They could be found in the source archive file for SUSY paper submitted to the HEP database.
of the polynomial Computation of the absolute error C[ j] break end of the loop on the polynomial C if not find if no similar terms are found Adding a new term to the polynomial C Lc ← Lc + 1 incrementing the length of the polynomial C C[Lc] ← B[i] copying the term B[i] into the last term of the C for j ← 1 to Lb Lb is the length of the polynomial B M ← A[i] B[j] M is the intermediate term Calculation of the control sum α m Calculation of the absolute error M Search for a similar term and adding M to the polynomial C ! zeroing the exponent of cosϕ α m [i, 7] ← 0 multiplying by π Calculation of the control sum α m Calculation of the absolute error M Search for a similar term and adding M to the polynomial C for i ← 1 to La La is the length of the polynomial A if α a [i, 3] > 0 then if the exponent of sinχ is greater than zero 0 FlagSin True the polynomial includes a term with the factor sinχM ← A[i] M an the intermediate term α m [3] ← α m [3] − 2 calculation ofthe term containing sin m−2 χcos n χ calculation of the control α m search for a similar term and addition of the term M to the polynomial C K m ← −K m α m [3] ← α m [4] + 2 calculation of the term containing-sin m−2 χcos n+2 χ calculation of the control sum α m search for a similar term and addition of the term M to the polynomial C A ← C the polynomial C is copied to the polynomial A Until FlagSin completion of the procedure if the polynomial does not contain terms with the factor sinχ

F
Elimination of zero terms Integration with respect θ for j ← 1 to Lu Lu is the length of the polynomial U M ← A[i] B[j] M is an intermediate term Integration of M with respect to ϕ Integration of M with respect to θ Integration of M with respect to x Calculation of the control sum α m Calculation of the absolute error M Search for a similar term and addition of M B to the polynomial F rn (y, χ).

HFigure 2 .Figure 4 .
Figure 1.One loop diagram with two mass scales m and M (mM) and small external momentum q (q 2 M 2 )

Table 1 .
Coefficients of the collision integral