Schwinger–Dyson equation for quarks in a QCD inspired model
Abstract
We discuss formulation of QCD in Minkowski–spacetime and effect of an operator product expansion by means of normal ordering of fields in the QCD Lagrangian. The formulation of QCD in the Minkowski–spacetime allows us to solve a constraint equation and decompose the gauge field propagator in the sum of an instantaneous part, which forms a bound state, and a retarded part, which contains the relativistic corrections. In Quantum Field Theory, for a Lagrangian with unordered operator fields, one can make normal ordering by means of the operator product expansion, then the gluon condensate appear. This gives us a natural way of obtaining a dimensional parameter in QCD, which is missing in the QCD Lagrangian. We derive a Schwinger–Dyson equation for a quark, which is studied both numerically and analytically. The critical value of the strong coupling constant , above which a nontrivial solution appears and a spontaneous chiral symmetry breaking occurs, is found. For the sake of simplicity, the considered model describes only one flavor massless quark, but the methods can be used in more general case. The Fouriersine transform of a function with logpower asymptotic was performed.
I Introduction.
Strong interaction physics should be described by Quantum Field Theory (QFT) with the Quantum Chromodynamic (QCD) Lagrangian Yang and Mills (1954); Faddeev and Popov (1967); Faddeev (1969); ’tHooft (1971); Gross and Wilczek (1973); Politzer (1973). As it shown Gross and Wilczek (1973); Politzer (1973), the running coupling constant is strong enough at small energies, so that perturbation expansion is not applicable. This is a big problem due to the lack of general methods of nonperturbative calculations.
In order to describe the strong interaction, the phenomenological models were developed that are based directly on experimental data and use partly the QCD knowledge: the QCD sum rules Shifman et al. (1979a, b); Ioffe (1981), the Chiral Perturbation Theory Volkov and Pervushin (1978, 1975, 1976); Gasser and Leutwyler (1984), the Nambu–JonaLasinio model and its generalizations Nambu and JonaLasinio (1961); Ebert and Volkov (1983); Volkov (1984, 1986); Reinhardt and Alkofer (1988); Bernard et al. (1988); Le Yaouanc et al. (1985), bag models Bogolubov (1968); Chodos et al. (1974a, b) and others Haymaker (1991); Pawlowski (2007). These models can relatively easily reproduce experimental data. However, they have a number of disadvantages: each of these models works in a certain application area but fail in others, the accuracy of theoretical calculations are limited and often less than the accuracy of modern experimental data. And these models are not true theory of strong interactions. This gives impetus to construct models based directly on QCD, for instance: instanton liquid model Shuryak (1982a, b, c); Dyakonov and Petrov (1984), domain wall network Efimov and Nedelko (1995); Burdanov et al. (1996); Kalloniatis and Nedelko (2001, 2004); Galilo and Nedelko (2011); Nedelko and Voronin (2015), various estimations from Schwinger–Dyson equations Roberts and Williams (1994); von Smekal et al. (1998); Alkofer and von Smekal (2001); Maris and Roberts (2003); Fischer (2006); Epple et al. (2007); Aguilar et al. (2008).
Thus, there exist various approximations to the theory of strong interactions with their specific simplifications of the QCD. In the QCD researches, one should find answers to the key questions, which are the description of QCD vacuum, spontaneous breaking of chiral symmetry, the absence of color particles (confinement problem), the description of bound states, their masses and decay widths.
We consider the theory of strong interaction at low energy. Our aim is to emphasize the importance of formulation in Minkowskispacetime and effect of an operator product expansion by means of normal ordering of fields in Lagrangian, and to discuss some consequences of this novel approach.
The formulation of QCD in the Minkowskispacetime allows us to solve a constraint equation and decompose the gauge field propagator in the sum of an instantaneous part, which forms a bound state, and a retarded part, which contains the relativistic corrections. At the first stage, we should neglect the retarded part and use the instantaneous part to construct the bound state. Then the retarded part gives corrections to the already existing bound state. This idea cames form QED Shilin and Pervushin (2013) (see also Polubarinov (2003); Pervushin (2003, 1990)), where any attempts of working with the entire propagator do not lead to satisfactory results or the decomposition occurs implicitly. Our approach enable us to cover the both high and lowenergy ranges and find the relation between fundamental QCD parameters and low energy constants.
In QFT, for a Lagrangian with unordered operator fields, one can make normal ordering by means of the operator product expansion. Then the gluon condensate and a low energy effective gluon mass appear. This mechanism gives us a natural and fundamental way of obtaining a dimensional parameter in QCD, which is missing in the QCD Lagrangian. The existence of nonzero condensates directly linked to the conformal anomaly of QCD.
In the next section, we start from QCD Lagrangian and derive an effective action of strong interaction. Then in section III, by using this effective action we obtain the Schwinger–Dyson equation for a quark, which is solved both numerically (in section IV) and analytically (in the subsequent sections). In conclusion, we summarize the obtained results and discuss the prospects of the developed methods. Here, for the sake of simplicity, we intentionally neglect some effects, for example, the considered model describes only one flavor massless quark. While investigating of the Schwinger–Dyson equation, we focus mainly on the question of the spontaneous symmetry breaking. Nevertheless, all the assumptions made to derive the equation are transparent and wellcontrolled.
Ii Effective Action for the Strong Interaction.
Let us start with the Quantum Chromodynamics Lagrangian, with number of colors and number of flavors :
(1) 
where , , , , and are the gluon field, gluon field strength tensor, quark field, quark current mass, and color current of quark, respectively.
An effective action for the mesonlike bound state can be derived from the Lagrangian (1). To this end, some restrictions and assumptions are needed. Below the symbol • is introduced for convenience when we discuss another one assumption or resrtiction. Some of the restrictions are not principle but imposed in order to not overload the reader by technical calculations. Anyway, in the developed model, we outline main ideas that may be important for correct description of mesonlike bound state rather than give a complete description of strong interaction, which certainly remains a tremendous problem.

First, we choose the frame of reference where the bound state, which we obtain and discuss below, is as whole at rest. Therefore only the static problems are considered. We emphasize that the proper choice of the reference frame should be done in Minkowskispacetime rather than in Euclidianspacetime. Note that the generalization of this theory to one bound state moving on mass shell Pervushin (2003); Cherny et al. (2013); Pervushin (1990); Pervushin et al. (2012, 2015) can easily be done, it is sufficient to rewrite various quantities in the comoving frame of reference.
We fix the gauge
(2) 
where and are the space and gluon color indexes, respectively.
The gluon term in the Lagrangian takes the form
(3) 
The third term on the right hand side contains the time derivative and thus can be neglected, because only the static problems are considered, as is noted above.
After quantization, the gluon field becomes an operator field. One can consider the vacuum 2point correlator
(4) 

We assume that and . Actually depends on the energy, but for simplicity we suppose to be a constant. The constant can be determined from a phenomenology.
The fields in Eq. (4) obeys the condition (2). A question at once arises: “How should the formula (4) be rewritten in other gauges?” The answer is to make a gauge transform to (2), then impose the condition (4), and then make the inverse gauge transform. As is said above, when solving the onshell bound state problem, we always have a privileged frame of reference, in which this bound state as whole is at rest; therefore, we always have the privileged gauge (2), and thus we define (4) in a gaugecovariant manner in this way. Note that physically privileged reference frame is absent for any scattering problem of quarks and gluons, and one cannot define (4) the same way.
Usually in Quantum Field Theory, a Lagrangian contains only normally ordered operator fields. This is a result of normal ordering of an initial Lagrangian where the above correlatorlike terms, arising due to the ordering, are omitted, because they are considered as (infinite) vacuum energy contributions. Keeping these terms, we have after the normal ordering
(5) 
The term with comes from the last term of formula (3) and . Here we use the relation . The quantity might be interpreted as an effective gluon mass in the gauge (2). This is essentially a modeldependent quantity. In this approach, the gluon mass appears before a perturbation expansion. Phenomenological models in which gluons have nonzero effective mass at small energies have been considered earlier by some authors (see Aguilar et al. (2008); Mandula and Ogilvie (1987); Amemiya and Suganuma (1999); Cornwall and Soni (1983); Bogolyubskaya et al. (1990); Halzen et al. (1993); Larin (2016a, 2015, b) and references therein).

Let us consider the dotted terms in (5) as a perturbation and neglect them. This assumption means that we suggest that bound states are formed by only some of the terms which explicitly written in expression (5), while the other terms merely give some corrections to the already existing bound states. In the basic model developed in this paper, these terms are neglected. The neglected terms can influence on quantitative characteristics of the bound states, but not their presence, and numerical amount of corrections might be not small due to large value of strong coupling constant.
Substituting (5) without dotted terms into Lagrangian (1), we arrive at the generating functional
The source is not involved, since the field is not dynamical degree of freedom with the gauge (2). This is owing to the fact that the corresponding equation of motion is a constraint Dirac (1955).
Making integration over yields
The term
includes a combination of the GellMann matrices, which may be rewritten in the form

We restrict ourselves to the colorless mesons and so neglect the second term. This term is the diquark channel, which plays a role when baryons are taken into account (”baryon = diquark + quark”).
Thus within this approximation, the above term can be rewritten in the form
where the above formula is the definition of the operator , and was shifted to the left. One can see that color indexes and have been summed inside pairs , so the pair as whole is colorless.

Let us treat as a real bilocal field.
The operator is symmetrical and has an inverse operator that can be defined by:
This allows us to introduce new bilocal field and make a bosonization transform (HabbardStratanovich transform) Ebert and Pervushin (2004, 1976); Kleinert (1976); Pervushin et al. (1979); Stratonovich (1957); Hubbard (1959):
Finally the generating functional for effective action of strong interaction takes the form
(6) 
With the help of this generating functional, one can write down any diagrams for processes of interest.
Iii Schwinger–Dyson equation.
In what follows we restrict ourselves only to the question of spontaneous symmetry breaking in the theory described by the functional (6). For this purpose, it is convenient to derive and investigate the Schwinger–Dyson (Gap) equation for the quark.
It is difficult to examine the Schwinger–Dyson equation in the general form.

For the sake of simplicity, we use the Stationary Phase method (that is the Semiclassical approximation). This method simplify the Schwinger–Dyson equation but retain its main properties.
According to this method, we should integrate out the Fermion variables and in (6), thus deriving the functional for the action
(7) 
Then we arrive at the Schwinger–Dyson equation
(8) 
which gives us the Fermion spectrum inside the bound state Pervushin et al. (1979); Ebert and Pervushin (2004, 1976); Pervushin et al. (1979, 1990); Kalinovsky et al. (1989a, 1990, b); Pervushin (2003).
We introduce the operator
and define its inverse as
In this notations, formula (7) reads
Inserting the corresponding into equation (8) we arrive at
(9) 
Below in this article, the solution of this equation is denoted by
It is convenient to introduce the operator
which, on the stationary solutions obeying Eq. (9), coincides with the earlier introduced operator
The inverse operator is defined in the standard manner
Acting with the operator on the both sides of Eq. (9) and using the above notations, we obtain
(10) 

We are looking for a simplest solution of this equation and adopt the following ansatz
Due to the isotropy, is radially symmetric and depends only on .
Making the Fourier transform of equation (10), we have
(11) 
In momentum space, the operator can easily be reversed
where we put by definition , and: . One can see that the timecomponent appears in Eq. (11) only through and, hence, can be integrated out.
After integrating over the solid angle in 3D momentum space, finally the Schwinger–Dyson equation takes the form
(12) 
with and being the absolute values of and , respectively.
As discussed above, Eq. (8) describes a fermion spectrum inside the bound state. Thus, the physical meaning of is a running quark mass, and, hence, it should be positive for any momentum. At , the value corresponds to a constituent quark mass, while the current quark mass is . One can introduce instead of quark charge a strong coupling constant . It is well known in QCD is a running coupling, whose value strongly dependent of the energy scale. Moreover, at low energies, the dependence of momentum can not be calculated from the perturbation theory, which is inapplicable due to the large value of . In the literature, there exist various predictions about the shape of (see, e.g., Pawlowski (2007); von Smekal et al. (1998); Alkofer and von Smekal (2001); Maris and Roberts (2003); Fischer (2006); Epple et al. (2007); Bogolubsky et al. (2009) and references therein). Nevertheless in this paper, we assume that is a constant; which is consistent with, as we mention above, neglecting corrections to the bound states; this means in particular neglecting all the loop corrections to , and is really a constant in the framework of this approach. So in a way, the used in this article constant can be understood as an average of the strong coupling over within a lowmomentum range.

We solve the equation (12) only for , which can be justified by the phenomenology. Indeed, , because the current mass of light and quarks is about . On the other hand, the constituent mass of the same quarks is of order for different models.

One can demand when . Due to the asymptotic freedom at large momenta, the running quark mass tends to the current mass. Although the existence of the asymptotic freedom in our model is questionable, we do not want to violate it explicitly. In addition, if this restriction is fulfilled then the equation (12) does not need any renormalization.
It is convenient to introduce the dimensionless variables , , and . In this variables, the Schwinger–Dyson equation (12) takes the form
(13) 
It is obvious that there always exists the solution . Of course, we are looking for a nontrivial solution of this equation.
An attractive feature of the Schwinger–Dyson equation (13) is that it is controlled only by one external parameter , which should be fixed from the phenomenology. In particular, it follows from the definition of the dimensionless variables and Eq. (13) that the constituent quark mass is a linear function of and the coefficient of proportionality depends only on : .
Iv Numerical solution of the Schwinger–Dyson equation.
To solve equation (13) numerically, we use the following algorithm. Let us take a zerothorder approximation function , it is desirable that differs from the solution not much. Then substituting into the integral in the righthand side (13) we get in the lefthand side. Then is substituted again, and so on. After a certain number of steps we get, up to the errors of computer calculations, the exact solution for which the substitution into the righthand side (13) gives itself. Strictly speaking, we should prove that this algorithm is convergent. We did not try to prove this because in all cases that we calculated this algorithm turned out to be convergent. Moreover, there is no difference in the choice of (see below for details).
In some sense, the convergence of the algorithm can be explained by the stability of the solution under small perturbations. Namely, substituting the function , where , into integral (13) we have up to terms
The last term is smaller than and has a minus sign. That is why the obtained expression is closer to the solution than .
After some attempts to solve equation (13) numerically, we have found that there are some things that should be avoided at numerical computation:

The upper limit of integration must be and cannot be replaced by finite quantity ; otherwise a strong dependence of the solution form appears.

, otherwise the integral diverges.

It is better to avoid replacing the continuous function by a discrete table with fixed numbers of points . That is because the value of at the penultimate point (at the last point , as it is noted above) depends mainly on the behavior of between this point and the end point and has a weak dependence on the values of a the other points; the value of at the next to penultimate point depends on the value of at the penultimate point and the behavior of between these three points, and so on. One cannot approximate the behavior of the function between two points by a linear segment, otherwise this leads to very low accuracy of numerical calculations. Preferably, expands in a series of known functions. One may also add that maybe we have more accurate results than in paper Puzynin et al. (1999), where a similar equation was considered numerically and such replacing by the table was done.
Put by definition , then equation (13) can be rewritten in the form
(14) 
Let us define the new function
(15) 
One can easily see that has the properties
(16)  
(17)  
New variables can be introduced (where – is some parameter)
In this variables the SchwingerDyson equation takes form
(18) 
On there is a convenient system of the Fourier series functions:
Using the Fourier series expansion we can avoid all the numerical difficulties which were discussed above. As the Fourier harmonics are periodical functions, it would be better if the area of the fastest change of the function lay closer to the center of the interval. The point corresponds to , so it dictates the choice of . Of course, before the calculation we do not know what value should be taken; fortunately, the incorrect leads only to hight inaccuracy and low speed of calculation. Equation (18) now takes the matrix form
(19) 
where , where