Gauge-invariant approach to quark dynamics

The main aspects of a gauge-invariant approach to the description of quark dynamics in the nonperturbative regime of QCD are first reviewed. In particular, the role of the parallel transport operation in constructing gauge-invariant Green's functions is presented, and the relevance of Wilson loops for the representation of the interaction is emphasized. Recent developments, based on the use of polygonal lines for the parallel transport operation, are then presented. An integro-differential equation is obtained for the quark Green's function defined with a phase factor along a single, straight line segment. It is solved exactly and analytically in the case of two-dimensional QCD in the large $N_c$ limit. The solution displays the dynamical mass generation phenomenon for quarks, with an infinite number of branch-cut singularities that are stronger than simple poles.


Introduction
Quantum chromodynamics (QCD), the theory of the strong interaction in which the fundamental fields are the quarks (the matter fields) and the gluons (the gauge fields), has the primary property of being asymptotically free [1,2]: The interaction weakens at short distances between particles or sources, a feature that allows the use of perturbation theory in that domain. This theoretical prediction has been widely verified by many experimental processes and data and represents one of the major justifications for the advent of QCD.
The counterpart of asymptotic freedom is that the effective coupling constant increases at large distances, where perturbation theory ceases to be valid. The failure of perturbation theory in the infrared domain actually has a deep origin, confirmed by nature: Quarks and gluons, the fundamental constituents of QCD, are not observed in nature as asymptotic free particles; at low energies or at large-distance separations, they are confined into color singlet bound states, which are the hadrons (mesons and baryons). Their presence there can only be detected indirectly through spectroscopic analyses. This phenomenon is a new feature that had not been addressed in earlier known field theories.
The precise description of the confinement mechanism requires a nonperturbative approach and still remains a challenge on analytic grounds.
The only systematic tool that succeeds in enabling the treatment of the large-distance regime of QCD is lattice theory, first introduced by Wilson [3,4,5]. Lattice theory is based on the discretization of spacetime and a numerical calculation of the path integral. Its limitations come from the facts that the theory is considered in Euclidean space, rather than Minkowski space, and the continuous Poincaré invariance is explicitly broken by the discretization mechanism and replaced by discrete symmetries. Furthermore, limitations of computational power also place restrictions on the precision of the calculation. It is worthwhile noticing, to the credit of lattice QCD, that it predicts analytically, in the strong-coupling approximation, the confinement of quarks.
Another systematic tool is provided by the Dyson-Schwinger equations [6,7,8,9]. These are integral equations for Green's functions of the theory, which, however, are infinite in number and coupled to each other. For a practical resolution, one considers mainly the equations relative to the simplest Green's functions with a truncation of the infinite series of coupled sectors. Numerical solutions compatible with confinement have been found. Also, solutions verifying chiral symmetry breaking have been obtained [10]. This line of approach is still under investigation. Inherent difficulties might be related with the gauge-invariance problem in nonperturbative approaches, where truncation of exact equations should be done in a consistent way in order not to break gauge invariance of physical quantities to the order of approximation that is considered. In QCD, because of the masslessness of the gluon field, artificial infrared divergences or singularities easily appear in intermediate quantities. Ensuring infrared finiteness of physical quantities is another difficult task in nonperturbative approaches.
Because of the infrared instabilities of noninvariant quantities under gauge transformations, a line of approach based on the use from the start of gauge-invariant quantities has been under development for a long time [11,12,13,14,15]. However, because gaugeinvariant Green's functions are extended composite objects, equations derived for them have more complicated structure and this fact has inhibited rapid and straightforward progress. Nevertheless, it turns out that the confinement problem has a sounder formulation in this approach than in more conventional approaches and provides simple criteria for its verification [3,16].
The remaining part of this article is devoted to a presentation of characteristic features of the corresponding formalism and of new developments obtained in recent years.

Gauge transformations
The gauge group of QCD is SU(N c ), acting in the internal space of color degrees of freedom; the physical value of N c is 3, but leaving N c as a free parameter allows one to study in more generality the properties of the theory.
The quark fields belong to the defining (N c -dimensional) fundamental representation of the gauge group. They are denoted ψ a α (x), where a, the color index, runs from 1 to N c ; α is the spinor index; and x is the spacetime coordinate, which we shall consider in Minkowski space. Antiquark fields are denoted ψ b,β (x); they belong to the conjugate (N c -dimensional) fundamental representation. Actually, there are six different types of quark, with different masses and appropriate electric charges, classified according to a global index, called flavor; however, for simplicity, we shall consider in the following only one type of quark, without specifying its flavor. (In QCD, the interaction does not couple different flavor states.) Under a local gauge transformation, the quark and antiquark fields transform as where Ω is an element of the group SU(N c ), whose parameters are in general x dependent (unless specified otherwise, repeated indices are assumed to be summed over).
Because of the x dependence of the parameters of the group transformations, the derivatives of the quark and antiquark fields do not transform as members of the fundamental representations, as in Eqs. (1). Fields or functions of fields that transform as being members of irreducible representations will be said to be transforming covariantly, although here the term covariance is not related to the position of the color index (upper or lower).
To define a covariant derivative operator one introduces a connection, represented by the gluon gauge field, that belongs, for x-independent gauge transformations, to the adjoint representation of the group; since the latter can be obtained from the composition of the fundamental representation with its conjugate, the gluon field can be represented as depending on two color indices, related to the two previous representations: A a b,µ [17]. One also has the properties A †a b,µ = A b a,µ and A a a,µ = 0. The more conventional notation represents the gluon field with a single index B running from 1 to (N 2 c − 1). The relationship between the two notations can be obtained in the following way. If T B are the generators of the group transformations in the fundamental representation, then we (The factor 1/ √ 2 comes from a normalization convention.) The covariant derivative acting on the quark field is where g is the (dimensionless) coupling constant of the theory.
As is the case of connections, the gluon field does not transform covariantly under gauge transformations; it obtains an inhomogeneous part that breaks covariance. It is possible, however, to construct a new field that has covariance properties. This is the gluon field strength F (the analog of the curvature tensor of differential geometry), given by Its transformation law is With the aid of the covariant fields, one can easily construct the gauge-invariant Lagrangian density of the theory: where m is the mass parameter of the quark fields and γ µ are the Dirac matrices (the spinor indices being omitted).

Parallel transport
Let us consider the ordinary two-point Green's function of a quark field: where the first expression in the right-hand-side represents the vacuum expectation value of the chronological product of the two quark fields, while the second expression corresponds to the notation used in the path-integral formalism; we recall that, in the latter formalism, ordinary products of functions actually represent chronological products of the operator formalism.
Considering the gauge transformation law of the quark fields as given in Eqs.
(1), we immediately deduce that the above Green's function is not invariant under local gauge transformations, since the vacuum state and the Lagrangian density (5) are gauge invariant. The consequence of this property can only be the vanishing of the Green's function for general values of the coordinates x and y [18]. This would mean the impossibility of treating perturbatively noninvariant quantities. The resolution of this difficulty comes from adding into the Lagrangian density (5) gauge-fixing terms and ghost fields [19]. The classical local gauge invariance is then replaced by a quantized global invariance, the BRST symmetry [20,21]. However, in gauge-invariant Green's functions the gauge-fixing and ghost-field contributions generally cancel each other and one may continue reasoning, at a primary level, with the classical gauge transformations.
To construct gauge-invariant Green's functions, one should relate in some way gauge transformations at the point x with those at the point y. This is generally done by devising an operation of parallel transport, for instance from point x to point y. This operation is usually realized with the aid of the connection. By considering the parallel transport along an oriented curve C yx , its representative function, which we denote U(C yx ; y, x), takes the form of a path-ordered phase factor of the gluon field: where P represents the path-ordering operation, meaning that the gluon fields are ordered according to their position on the curve C yx . The above exponential can also be decomposed into products of other exponentials, each defined successively on a smaller part of the curve C yx . A more explicit expression of U can be obtained by a series expansion of the exponential in terms of the coupling constant, in which path ordering is taken into account. Parametrizing the function z on the curve with a parameter λ varying between 0 and 1, such that z(0) = x and z(1) = y, and defining z ′ (λ) = ∂z(λ) ∂λ , one obtains where the abbreviation A(z(λ)) = A(λ) has been used.
Under a gauge transformation, U behaves as a covariant bilocal field, the transformation depending only on its end points x and y, independently of the form of the curve C yx [22]: By applying the operator U on a quark field at point x, its variation at x cancels the one coming from the quark field and the whole object transforms as belonging to the fundamental representation at point y: Finally, by applying the previous object from the left on an antiquark field, the ydependent variation of the latter is canceled and the whole quantity becomes invariant under gauge transformations: Considering this object in the path-integral formalism, or taking the chronological product of the fields once their color indices are fixed, allows one to define, up to a normalization factor, a gauge-invariant quark Green's function. The latter depends, in addition to the coordinates of the quark fields, on the line C yx that appears in the definition of the parallel transport operation from point x to point y. We shall denote the resulting Green's function S(x, y; C yx ); it no longer depends on color indices and is given by

Wilson loops
Another gauge-invariant useful object is constructed with the sole parallel transport operator (7). Considering in the latter the case of a closed line C xx , obtained from C yx by the limit y → x along a nontrivial curve, and then taking the trace of the color indices at the matching point, one clearly obtains from Eq. (9) a gauge-invariant result: The line C xx being closed and the trace having been taken, the point x no longer plays a particular role in that object and hence can be removed from the representation of the operator, where only the closed contour may appear. The resulting object is called a Wilson loop: Its vacuum expectation value [3] will be denoted by The Wilson loop provides a simple criterion for the confinement of quarks. In the case of a static pair of a quark and an antiquark, considered as the extreme limit of heavy quarks, which remain fixed in position space and which are observed at equal times, the whole dynamics of the system in the color singlet state can be described by means of a Wilson loop average (15) along a rectangular contour in a plane, where one of the sides represents the fixed distance R between the quark and the antiquark and the other side the interval T of the evolved time (see Fig. 1). For large separation distances and large Figure 1: Rectangular contour of a Wilson loop representing the evolution of a static quark-antiquark system. time intervals, the following relationship can be shown [16]: where V (R) is the static potential energy corresponding to the force that is exerted between the quark and the antiquark. However, W (C rect. ) can be analytically calculated in lattice QCD [3,4,5] in the strong-coupling approximation, yielding where σ is a positive constant with dimension energy/length. Comparing the two expressions (16) and (17), one deduces which shows that the interquark static potential energy is a linearly increasing function of the distance. This is sufficient to bind the quark and the antiquark together and to prevent their separation to infinite distances. Since the product RT represents the area of the rectangular contour, the above result is called the area law of confinement. Assuming that the previous result remains in general true in the continuum limit of lattice QCD, one ends up with the conclusion that QCD in the continuum should confine quarks, in accordance with the observational facts. The constant σ is called the string tension; the numerical value of √ σ lies within the bounds 420-480 MeV (with units whereh = c = 1) [23]. The large-time behavior in the static limit of the Wilson loop average determines therefore the interquark potential energy. Because the latter is a quantity defined mainly in nonrelativistic theories, it appears that the Wilson loop itself should be the quantity that would play the role of potentials in the general case. Being gauge invariant, the information coming from it about the dynamics of quarks would have an unambiguous interpretation. In the case of general simple contours, the area law obtained for the rectangular contour would be replaced by the area of the surface, supported by the contour, having a minimal value (called a minimal surface) [3]. Wilson loops are often used in relativistic approaches to bound state problems [24,25,26,27,28].
Wilson loop averages can be considered as functionals of the contour C from which they are defined. They may therefore receive a geometric interpretation or be subjected to a geometric analysis. Properties of Wilson loops have been thoroughly studied, starting with works by Polyakov [29] and Makeenko and Migdal [30,31,32,33].
The analysis is based upon functional derivations acting on the path-ordered phase factor (the parallel transport operator), which may be submitted to local deformations of its line. Two main equations are obtained, the Bianchi identity and the loop equation, the latter resulting from the equation of motion operator of the gluon field. Owing to the complexity of the loop equation, obtaining an exact nonperturbative solution of it displaying confinement has not been possible. However, Makeenko and Migdal showed that, for large simple contours, the loop equations have asymptotic solutions that satisfy the area law of confinement [31], obtained independently in lattice QCD [3,4,5]. In perturbation theory, renormalization of Wilson loops was shown in [34,35].
The loop equation was also analyzed in two spacetime dimensions by Kazakov and Kostov [36,37], who could exactly solve the equation for many kinds of contour. In particular, for simple contours, the solution is given by the exponential of the area of the surface enclosed by them: The area law is thus explicitly satisfied. The same results were also obtained by Bralić [38], who used the non-Abelian version of the Stokes theorem for the calculation of the Wilson loops. Actually, the main simplification in two dimensions arises from the fact that the QCD coupling constant g has now a dimension of mass (h = c = 1) and its square has the dimension of the string tension that appears in the area law; this feature naturally drags the solutions to the area law, which is also obtained in perturbation theory in the axial gauge.
A general survey of questions related to loop equations and gauge theories can be found in [39]. Properties of Wilson loops saturated with minimal surfaces are studied in [40]. It has been shown that the minimal surface is the only type of surface that satisfies the Bianchi identity. This implies that fluctuations of surfaces around a minimal surface violate in general the Bianchi identity [32,40]. The violation of the Bianchi identity is a signal of the presence of a color monopole current. However, the Bianchi identity plays an important role in the verification of some consistency conditions of the loop formalism, in particular for ensuring the commutativity of two successive functional derivatives [40]. It has also been shown that, by neglecting the short-distance (i.e., perturbative) interactions and renormalizing the coupling constant g 2 N c quadratically with respect to a short-distance regulator, one ends up, for simple contours, with a solution to the loop equation that is identical to a minimal surface [40]. This provides another verification of the asymptotic solution obtained in [31]. (In perturbation theory, g 2 N c vanishes with the short-distance regulator as the inverse of a logarithm.)

Polygonal lines for parallel transport
Although loop equations have not been solved exactly or analytically in four spacetime dimensions, the past investigations provide us with a general idea about the main features of Wilson loop vacuum averages. For simple large contours, minimal surfaces are asymptotic solutions, in accordance with the area law obtained in lattice QCD in the strong-coupling regime. For small contours, ordinary perturbation theory should be applicable for their evaluation. Adopting this qualitative basis, one can go further, using Wilson loops in dynamical equations where they play the role of gauge-invariant potentials or kernels. Since most experimental observations related with QCD concern hadrons, which involve quarks as main degrees of freedom, the study of the properties of gauge-invariant quark Green's functions appears as a natural step for such an investigation. We already defined in Eq. (12) the two-point gauge-invariant quark Green's function (2PGIQGF), which is a functional of the line defining the parallel transport from the quark field to the antiquark field.
At this point the question arises as to the impact of the type of phase factor line on physical quantities: Do the latter depend on the choice made of the line followed by the parallel transport operation? A general answer to this question has been given in Ref. [40] (end of Appendix A), in the simplified case of a static quark-antiquark system evolving in time, within the framework of the saturation of Wilson loop averages by minimal surfaces. One compares the evolution law of a system defined with a straight line segment for the phase factor line with that of a system defined with an arbitrary (continuous) line. The physical observable is represented here by the energy content of the system. The latter is explicitly exhibited through an exponential factor when the evolution time interval T tends to infinity (see Eq. (16)). The equations of the minimal surfaces allow us to study in detail this limit: Both minimal surfaces yield the same energy dependence, governed actually by that of the straight line segment and the rectangular Wilson loop (Fig. 1). The line dependencies are factored out within the wave functional of each system. The above analysis displays the effect the choice of the phase factor line may have on a physical system under study: It affects the wave functional of the system, but it has no influence on physical quantities. The case of nonstatic quarks does not seem to modify in general the previous conclusion and we adopt the latter in the forthcoming part of our investigation.
Since the simplest choice for a phase factor line is the straight line segment, one would be tempted to stick to Green's functions defined with such lines. This is, however, not possible as a single stage, because the equations of motion couple, through interaction terms, Green's functions with different lines. A single choice of a rigid line does not lead to a closed system of equations. The situation here is very similar to what happens with the Dyson-Schwinger equations, which consist of an infinite number of coupled equations. In the present case, in view of the complexity of all types of line, it is natural to seek classes of line that might form a closed system under the equations of motion. This problem actually has a positive answer: Polygonal lines, made of a succession of straight line segments, joined to each other, do form a closed set, at least for studies related to quark Green's functions [41].
Polygonal lines have several advantages worth mentioning: (a) They can be classified according to the number of segments they contain; accordingly, the same classification appears also at the level of the corresponding Green's functions. (b) Their basic ingredients, which are the straight line segments, are Lorentz invariant in form, which also reflects itself on the previous classification scheme. (c) The straight line segment has a well-defined, unambiguous limit when its two end points approach each other: It shrinks to a point. (d) Polygonal lines form a complete set for the present problem, in the sense that all equations that are derived will be closed within the category of polygonal line; no other types of line will be needed to complete the description. We shall henceforth consider for the parallel transport operator the class of line having a polygonal structure in spacetime.
One of the key ingredients in our study is a functional variation formula, first established by Mandelstam [14], which results when local deformations of the line C yx are considered. If the line C yx is deformed at all its points by infinitesimal variations δz(λ), including displacements of its end points, then the parallel transport operator undergoes the variation [14,15,22,29,42] where the abbreviation U(λ 2 , λ 1 ) has been used to designate the operator U corresponding to the line going from z(λ 1 ) to z(λ 2 ) along the curve C yx , with z(0) = x and z(1) = y. This formula exhibits the relationship between dynamical operations, such as the insertion of the gluon field strength inside the phase factor, and geometric operations, such as line deformations. We now consider for the line C yx an oriented straight line segment going from x to y and designate by U(y, x) the corresponding parallel transport operator. A displacement of one end point of the rigid segment in form, while the other end point remains fixed, generates a displacement of the interior points of the segment. This defines a rigid path displacement. By parametrizing the interior points of the segment with a linear parameter λ varying between 0 and 1, such that z(λ) = λy + (1 − λ)x, the rigid path derivative operations with respect to y or x yield where the second terms on the right-hand sides represent the contributions of the interior points of the segment and are given by the following expressions: The superscript + or − on the derivative variable takes account of the orientation on the segment and specifies, in the case of joined segments, the segment on which the derivative acts.
In the case of a polygonal contour C n , with n segments and n junction points x 1 , x 2 , . . ., x n , the vacuum average of the Wilson loop will be designated by W n = W (x n , x n−1 , . . . , x 1 ) (see Fig. 2). and (n − 1) junction points, the corresponding Green's function will be designated by

Quark Green's functions
where each phase factor U is along a straight line segment indicated by its two end-point coordinates. (Spinor indices are omitted and the color indices are implicitly summed.) The simplest such function is S (1) , having a phase factor along a single straight line segment: (We shall generally omit the index 1 from that function.) Figure 3 represents pictorially two different Green's functions. For the internal parts of rigid path derivatives, we have definitions of the typē The Green's functions satisfy the following equations of motion concerning the quark field variables: For n = 1, we have Similar equations also hold with the variable x ′ . A graphical representation of the equations of motion of S (1) and S (3) is displayed in Fig. 4. Figure 4: Graphical representation of the equations of motion of S (1) and S (3) . The cross on the dashed lines represents the rigid path derivation; it is placed near the end point of the segment that is submitted to the derivation.
Multiplying Eq. (27) with S(t 1 , x) and integrating with respect to x, one obtains functional relations between various 2PGIQGFs. A typical such relation is (Integrations on intermediate variables are implicit. Here, y 1 is an integration variable.) A graphical representation of this equation for n = 3 is given in Fig. 5.
Since this equation is valid for any n ≥ 1, one can use it again on its right-hand side for S (n+1) . One therefore generates an iterative procedure that eliminates successively the higher index 2PGIQGFs in favor of the lowest index one, S (1) . Assuming that the Figure 5: Graphical representation of the functional relation among S (3) , S (1) , and S (4) .
terms rejected to infinity are negligible, one ends up with a series where only S (1) appears together with Wilson loop averages along polygonal contours with an increasing number of sides and rigid path derivatives along the segments. This result shows that, among the set of the 2PGIQGFs S (n) , n = 1, 2, . . ., it is only S (1) , having a phase factor along one straight line segment, that is a genuine dynamical independent quantity. Higher index 2PGIQGFs could in principle be eliminated in terms of S (1) and polygonal Wilson loops and their rigid path derivatives.
The above procedure is also repeated on the right-hand sides of the equations of motion (27) and (28) to express the rigid path derivatives in terms of Wilson loop averages. The equation of S takes at the end the following form: K nµ− (x ′ , x, y 1 , . . . , y n−1 ) S (n) (y n−1 , x ′ ; x, y 1 , . . . , y n−2 ) , where the kernels K n (n = 1, 2, . . .) contain Wilson loop averages along polygonal contours, which are at most (n + 1)-sided, and the 2PGIQGF S and its derivative. The total number of derivatives contained in K n is n. Once the Wilson loop averages and the various derivatives have been evaluated and the high-index S (n) s have been expressed in terms of S, Eq. (30) becomes an integro-differential equation in S, which is the primary unknown quantity to be solved. This equation is the analog of the self-energy Dyson-Schwinger equation for ordinary Green's functions. Furthermore, the fact that it has been obtained with the sole aid of polygonal lines for the parallel transport operator, without the need of other types of line, is an indication that polygonal lines form a complete set for the study of the 2PGIQGFs.
The various kernels that are present on the right-hand side of Eq. (30) can be analyzed in terms of their behaviors at large and short distances. It seems that the series is globally perturbative with respect to the inverse of the number of sides of the polygonal contours, the first terms, having the least number of sides, being the dominant ones. The kernel K 1µ− is actually null for symmetry reasons; it is therefore the kernel K 2µ− , corresponding to a triangular contour with two rigid path derivatives, that is expected to be the dominant term of the series.
One of the mathematical difficulties of Eq. (30) comes from the property that the kernels do not have a convolutive structure; this is related to the fact that Wilson loops act in each term on all the points that are present in the corresponding expressions. This is actually reminiscent of the difficulty of obtaining functional inverses of gauge-invariant Green's functions because of their nonlocal structure.

Two-dimensional QCD
The equations obtained in the previous sections also remain valid in two spacetime dimensions and could be analyzed more easily in this case. Two-dimensional QCD in the large-N c limit [17,43] provides a simplified framework for the study of the confinement properties, which are expected to be qualitatively similar to those of four dimensions. Wilson loop averages were explicitly calculated in two dimensions [36,37,38]: For simple contours they satisfy the area law. In that case, the second-order derivative of the logarithm of the Wilson loop average reduces to a two-dimensional delta function. Higher order derivatives give zero in Eq. (30), since they act there on different segments of the polygonal contours. Overlapping self-intersecting surfaces, which give more complicated expressions, are assumed to give negligible contributions, for the residual terms they produce are probably of zero weight under the integrations that are involved.
In the series of terms of Eq. (30) it is only the kernel K 2 that survives and the integro-differential equation takes the following expression [44]: where σ is the string tension. The above equation can be analyzed by first passing to momentum space. Designating by S(p) the Fourier transform of S(x), one can decompose it into Lorentz-invariant components: The solution of Eq. (31) can be searched for by using the analyticity properties of the 2PGIQGF. It turns out that the equation can be solved exactly and in analytic form. The functions F 1 and F 0 are found to have an infinite number of branch cuts located on the positive real axis of p 2 (timelike region), starting at thresholds M 2 1 , M 2 2 , . . ., M 2 n , . . ., with fractional power singularities equal to −3/2. Their expressions, for complex p 2 , are The Green's function S [Eq. (32)] then takes the form The masses M n (n = 1, 2, . . .) are positive, greater than the free quark mass m, and ordered according to increasing values. For massless quarks they remain positive. The masses M n and the coefficients b n , the latter also being positive, satisfy, for general m, an infinite set of algebraic equations that are solved numerically. Their asymptotic values, for large values of n such that n ≫ m 2 /(πσ), are The functions (M 2 n − p 2 ) −3/2 are defined with cuts starting from their branch points and going to +∞ on the real axis; they are real below their branch points on the real axis down to −∞.
Expressions (33) and (34) represent weakly converging series. The high-energy behavior of the functions F 1 and F 0 is obtained with a detailed study of the asymptotic tails of the series and the use of the asymptotic behaviors of the parameters M n and b n [Eqs. (36)]. One finds that they behave asymptotically as in free field theories, which is here a trivial manifestation of asymptotic freedom [45], since in two dimensions and in the large-N c limit the coupling constant does not undergo renormalization.
In summary, the solution of Eq. (31) is nonperturbative and infrared finite. The masses M n are dynamically generated, since they do not exist in the Lagrangian of the theory. They could be interpreted as dynamical masses of quarks with, however, the following particular features. First, they are infinite in number. Second, they do not appear as poles in the Green's function but rather with stronger singularities. In x space the latter do not produce finite plane waves at large distances and therefore quarks could not be observed as free asymptotic states. Nevertheless, the above singularities, being gauge invariant, should have physical significance and would show up in the infrared regions of physical processes involving quarks.
The fact that the singularities of the Green's function appear only in the timelike region of real p 2 is an indication that the quark and gluon fields satisfy, even in the nonperturbative regime, the usual spectral properties of quantum field theory [46,47,48]. Expression (35) can be interpreted as fitting a generalized form of the Källén-Lehmann representation [41,44,49,50], where the denominator of the dispersive integral now has a fractional power, while the spectral functions are saturated by an infinite series of dynamically generated single quark states with alternating parities. The latter still satisfy Lehmann's positivity conditions [50].
Finally, the question may arise as to the dependence of the predictions so far obtained on the choice of polygonal lines for the phase factor paths. According to our discussion at the beginning of Sec. 5, physical quantities, such as masses and momentum space singularities, should be insensitive to changes of line type and therefore should be line independent.

Conclusion
Gauge-invariant Green's functions are adequate tools for a systematic investigation of the nonperturbative properties of QCD. Polygonal lines, supporting parallel transport in quark Green's functions, form a complete set of paths for the study of the dynamical properties of the latter. An equation playing the same role as the self-energy Dyson-Schwinger equation for ordinary Green's functions has been obtained, in which the kernels are represented by Wilson loop averages along polygonal contours with rigid path derivatives on their segments. It is expected that the leading term in the series of kernels is provided by the Wilson loop along a triangular contour with two derivatives.
The application of this equation to two-dimensional QCD in the large-N c limit provides an exact nonperturbative analytic solution, not known from conventional approaches, that displays a dynamical generation of a series of massive quark states, with the characteristic feature that their singularities in momentum space are stronger than simple poles.
The consistency of the results obtained in two-dimensional QCD is a positive test for the general approach presently developed for investigations in four dimensions.