{'fraction_non_alphanumeric': 0.03449746940991681, 'fraction_numerical': 0.012177858791134208, 'mean_word_length': 4.391531625718766, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hadronic unquenching effects in the quark propagator

We investigate hadronic unquenching effects in light quarks and mesons. Within the non-perturbative continuum framework of Schwinger-Dyson and Bethe-Salpeter equations we quantify the strength of the back reaction of the pion onto the quark-gluon interaction. To this end we add a Yang-Mills part of the interaction such that unquenched lattice results for various current quark masses are reproduced. We find considerable effects in the quark mass function at low momenta as well as for the chiral condensate. The quark wave function is less affected. The Gell--Mann-Oakes-Renner relation is valid to good accuracy up to pion masses of 400-500 MeV. As a byproduct of our investigation we verify the Coleman theorem, that chiral symmetry cannot be broken spontaneously when QCD is reduced to 1+1 dimensions.


I. INTRODUCTION
Dynamical chiral symmetry breaking is one of the most important properties of lowenergy QCD. The breaking pattern has profound impact for phenomenological quantities, as e.g. the appearance of the pseudoscalar Goldstone bosons in the chiral limit of QCD and the non-degeneracy of chiral partners. Chiral perturbation theory [1,2] describes these effects very efficiently on the level of hadrons but has nothing to say about the underlying structure of the full theory. In this work we investigate dynamical chiral symmetry breaking on the level of (non-perturbative) correlation functions between quarks and gluons.
We are particularly interested in the interplay between the fundamental quark and gluon degrees of freedom and the resulting bound states. In full QCD there are hadronic contributions to the fully dressed quark-gluon interaction. These effects are generated by the inclusion of dynamical sea quarks in the quark-gluon interaction and are therefore only present in unquenched QCD. In this work we concentrate on the back reaction of the pions on the quark propagator and investigate the impact of the unquenched interaction on the quark-pion system.
Pion effects on the quark propagation are important for several reasons. They account for (at least part of the) pion cloud effects in baryons and mesons. Thus they need to be incorporated in bound state calculations aiming at an adequate description of phenomenological properties of these objects. Furthermore they allow for the possibility of hadronic intermediate states in bound state calculations and therefore generate the finite width of meson spectral functions. On a more fundamental level Gribov argued that pion effects in the quark propagator may be responsible for quark confinement [3]. In addition, the inclusion of these effects in a continuum framework allows for an extrapolation of unquenched lattice results for large quark masses towards the physical up/down quark values. Finally, it will become apparent that Goldstone effects are essential in smaller dimensions. In two dimensions they prohibit the dynamical breaking of a continuous symmetry [4].
We study the pion back reaction on the quarks in the Green's function approach to Landau gauge QCD using Schwinger-Dyson and Bethe-Salpeter equations (SDE/BSE) [5,6,7,8]. In the past these effects have been investigated to some extent within the NJL-model [9,10,11]. Here we are extending these model studies to full QCD. Unquenching effects in the quark-antiquark system have also been investigated by Watson and Cassing within the SDE/BSE approach in [12]. There a coupled system of equations for the Bethe-Salpeter kernel and the connected quark four-point function has been solved in a certain approximation, which neglects the back reaction of the hadronic resonances on the quark propagator. Our study here is therefore in some sense complementary with their work. It is also complementary to the investigations reported in [13,14], where unquenching effects in the gluon polarization have been considered.
In order to make our results as concise as possible we use a combination of two methods. We devise a truncation scheme for the combined Schwinger-Dyson equations for the gluon, ghost and quark propagator as well as the quark-gluon vertex, which is combined with the corresponding Bethe-Salpeter equation for light mesons. The details of the quark-gluon interaction are fixed such that the available unquenched lattice results for the quark propagator [15] are reproduced. This combination of methods allows us to study the relative impact of the pion contributions compared to the pure Yang-Mills part of the quark-gluon interaction.
The paper is organized as follows. In subsection II A we discuss the hadronic contributions in the quark-gluon vertex and detail our approximation scheme of this vertex and the quark Schwinger-Dyson equation (SDE). In subsection II B we explain the extraction of pion contributions in the quark-antiquark scattering kernel of the corresponding Bethe-Salpeter equation (BSE). The resulting interaction is supplemented by the familiar one-gluon interaction from the rainbow-ladder approximation. Both types of interaction together represent an approximation to the full unquenched quark-gluon interaction, which retains important properties of the mechanism of dynamical chiral symmetry breaking. A different perspective to this approximation in the framework of 2PI effective actions is discussed in subsection II C. In subsection II D we explain our procedure to fit the parameters of the one-gluon exchange part of the interaction such that results from unquenched lattice simulations at large quark masses are recovered. Some technical details of our numerical method are presented in III A. In subsection III B we discuss our numerical results for the quark propagator and the pion at physical quark masses. We then investigate the chiral limit in subsection III C. We show that the Gell-Mann-Oakes-Renner relation is satisfied in our approach and give a value for the chiral condensate. In subsection III D a preliminary analysis for the analytical structure of the resulting quarkpropagator is presented. Our results are summarized and discussed further in section IV. Throughout the paper we work in Euclidean space and in the isospin symmetric limit of equal up and down quark masses.

II. THE QUARK-GLUON INTERACTION
As emphasized above there is a distinct imprint of dynamical chiral symmetry breaking on the low-energy meson spectrum: the lowest-lying pseudoscalar meson states are identified with the (pseudo-)Goldstone bosons of this mechanism. In the chiral limit, these states are exactly massless although they can also be described as bound states of constituent quarks with masses of the order of 300-400 MeV. This dichotonuous nature of the pseudoscalar states has been discussed from first principles in [5]. It is the axialvector Ward-Takahashi identity (axWTI), relating the quark self energy and the quark-antiquark scattering kernel, that enforces a binding energy of the pseudoscalar meson system which exactly cancels the quark and antiquark masses. It is therefore mandatory for every meaningful computational scheme of QCD to respect the axWTI.
A widely used practical truncation of the SDE/BSE framework with this property is the rainbow-ladder approximation for the quark interaction. It has a history of remarkable successes (summarized e.g. in [7]). However, there are also shortcomings that limit the credibility of such an approximation. Consequently, several efforts have been made to extend this scheme, see e.g. [14,16,17,18,19,20]. In the following we will extend the rainbow-ladder scheme by taking into account additional hadronic contributions to the quark-gluon interaction.
A. SDEs for the quark-gluon vertex and the quark propagator Let us first briefly sketch the so-called renormalization-group-improved rainbow-ladder truncation. The starting point is the SDE of the quark propagator where S −1 0 (p) = ip · γ + m denotes the inverse bare quark-propagator, while S −1 (p) = ip · γA(p 2 ) + B(p 2 ) is the dressed propagator and Z 2 the renormalization factor of the quark field. The quark self energies A(p 2 ) and B(p 2 ) can be recombined into the quark mass M(p 2 ) = B(p 2 )/A(p 2 ) and the quark wave function Z f (p 2 ) = 1/A(p 2 ). The quark self energy is given by with k = p − q, the Casimir C F = (N 2 c − 1)/(2N c ) and the renormalization factor Z 1F of the quark gluon vertex. The self energy depends on the fully dressed quark-gluon vertex Γ ν (q, k) and the gluon propagator with the gluon dressing function Z(k 2 ). Throughout this paper we will work in the Landau gauge, which guarantees the transversality of (3). The rainbow-ladder approximation amounts to the replacement where Γ(k 2 ) can be viewed as a combination of the gluon dressing function and a purely k 2 -dependent dressing of the γ ν -part of the quark-gluon vertex. Often the function Γ(k 2 ) is referred to as an effective running coupling. Indeed it can be shown, that the ultraviolet behavior of Γ(k 2 ) has to resemble that of the strong running coupling in order to reproduce the correct ultraviolet behavior of the quark propagator as known from resummed perturbation theory and the operator product expansion [22]. The rainbow approximation then gives reliable results for the quark propagator in the ultraviolet. In the infrared, however, one needs to go beyond simple approximations like (4). The infrared shape of the quark propagator and also its analytical structure do depend on the details of the quark-gluon interaction [13,23]. In particular, tensor structures other than γ ν seem to be important there [24] The Schwinger-Dyson equation for the quark-gluon vertex is given diagrammatically in Fig. 1. In the first line we show the full, untruncated equation. For very small momenta, a selfconsistent solution to this equation has been given in ref. [24]. Here we are primarily interested in the mid-momentum behavior of the vertex and in particular in hadronic contributions. To lowest order in a skeleton expansion such contributions can only occur in the diagram with the bare quark-gluon vertex at the external gluon line. In the second line of Fig. 1 we expand the quark-antiquark scattering kernel of this diagram in terms of resonance contributions to the kernel and one-particle irreducible Green's functions ('skeleton expansion'). Amongst other terms one obtains one-meson and one-gluon exchange between the quark and anti-quark lines. Also diquark exchange contributions arise which are not explicitly shown in Fig. 1. The first baryon exchange diagram shows up as a 'two-loop' diagram, which also involves diquarks and the baryon-quark-diquark vertices. Note that double counting is trivially avoided in this combined expansion due to different quantum numbers in the exchange channel.
The computation of the hadronic diagrams in Fig. 1 is rather involved. The mesonexchange diagram requires in addition the solution of a coupled system of the Schwinger-Dyson equation for the quark propagator and a corresponding Bethe-Salpeter equation (BSE) for the meson-quark vertex. Even more complex is the baryon-exchange diagram, which involves the computation of the diquark-BSE and a Faddeev-type equation for the baryon bound state. In this work we therefore wish to concentrate solely on contributions from the one-pion exchange diagram in addition to the non-hadronic contributions. In fact this choice should be regarded as a first approximation to the fully unquenched system, since all other hadronic diagrams are at least suppressed by factors of Λ 2 QCD /m 2 H with H ∈ {K, ρ, N, ...}. The same is true for diquark exchange contributions.
On the quark level, the one-pion exchange diagram in Fig. 1 involves a closed quark loop and therefore can only appear in unquenched QCD. Our investigation therefore complements previous studies of unquenching effects in the SDE framework, where quark- loops in the gluon polarization have been investigated [13,14]. The Yang-Mills part of the vertex is conveniently approximated by where we denote the quark momenta by p 1 and p 2 and the gluon momentum by p 3 . The explicit expression for Γ Y M (p 2 3 ) is detailed in subsection II D. The ansatz (5) is an example for the rainbow-ladder approximation, cf. Eq. (4). A similar form has been used in [25,26], where quenched lattice results for the quark propagator have been reproduced. Note, however, that (5) involves only the γ ν -part of the tensor structure of the full vertex.
In general one can decompose the quark-gluon vertex Γ ν (p 1 , p 2 , p 3 ) into twelve different tensor structures, given by The details of this basis are given in ref. [27]. Here we only note that L 1 ν = γ ν ; the explicit forms of the other structures are not needed. In principle, all twelve tensor structures can be important in the intermediate momentum regime. This has been explored to some extent in lattice simulations [28]. We have verified by projection methods that the one-pion exchange diagram of Fig. 1 contributes to all these structures.
These considerations lead to the following approximation scheme for the quark-gluon vertex: we subsume the Yang-Mills part of the quark-gluon interaction into a form as given in Eq. (5) and add the one-pion exchange diagram of Fig. 1. The resulting vertex is given in Fig. 2. It contains pion contributions in all twelve dressing functions λ i and contributions from the Yang-Mills sector in λ 1 . A conceptually different justification of this approximation scheme in terms of a 1/N c -expansion is given in subsection II C. Certainly, an explicit diagrammatic calculation of the Yang-Mills part of the vertex would be preferable to an ansatz of the form of Eq. (5). Such a calculation is under way and will be detailed elsewhere.
Inserting the vertex of corresponding kernel for the Bethe-Salpeter equation of the pion we need to perform an additional approximation as detailed in the next subsection.
B. The quark-antiquark scattering kernel We start with the general expression for the homogeneous Bethe-Salpeter equation (BSE) for quark-antiquark bound states, which can be written as Here Γ is the Bethe-Salpeter vertex function of a quark-antiquark bound state and K is the Bethe-Salpeter kernel. The momenta k + = k + ξP and k − = k + (ξ − 1)P are such that the total momentum is given by P = k + − k − . The momentum partitioning parameter ξ reflects the arbitrariness in the relative momenta of the quark-antiquark pair and can be set to ξ = 1/2 without loss of generality. Latin indices (t, u, r, s) refer to color, flavor and Dirac structure. The BSE is a parametric eigenvalue equation with discrete solutions P 2 = −M 2 n where M n is the mass of the resonance. The lowest mass solution corresponds to the physical ground state. Since P 2 is negative, the momenta k ± are necessarily complex in Euclidean space and thus the quark propagator functions are evaluated with complex argument. We will come back to this issue below.
Chiral symmetry constrains the Bethe-Salpeter kernel K tu;rs via the axialvector Ward-Takahashi identity (axWTI), which relates the kernel to the quark self energy In the rainbow-ladder approximation this relation can be satisfied easily. For our Yang-Mills part of the interaction, Eq. (5), one then obtains the kernel The resulting BSE includes an effective one-gluon exchange between the quark-antiquark pair as shown in Fig. 4. The construction of a Bethe-Salpeter kernel corresponding to the pion-exchange part of our interaction is more complicated. We were not able to find an exact solution of the axWTI for the interaction given in Fig. 2. In principle one then has two choices: either one works with an approximate solution of the axWTI or one approximates the interaction further such that an exact solution of the axWTI to the approximated interaction is possible. Since the first option destroys the Goldstone character of the pion, only the second choice is viable and shall be followed here.
To this end we reconsider the quark-SDE given in Fig. 3. For diagrammatical reasons the corresponding Bethe-Salpeter kernel should look like a superposition of the effective one-gluon exchange and an interaction involving pions. Let us assume for the moment that only the Yang-Mills part of the interaction would be present in the BSE. Then the quark-SDE could be rewritten by inserting the BSE into the second diagram on the right hand side. The result is given diagrammatically in Fig. 5 and represents an interaction given by an effective one-gluon and a one-pion exchange diagram in the SDE.
The only remaining question concerns the description of the one-pion exchange, which will be needed in the t-channel (see Fig. 3 and Fig. 5). In the s-channel the pion contribution to the quark-antiquark scattering kernel is given by [5] M pion tu;rs (q, p; P ) = [Γ i π ] rs (q; −P ) where Γ i π is the Bethe-Salpeter vertex function,Γ i π (q; −P ) = C −1 Γ i π (−q; −P ) T C and C = γ 2 γ 4 . As already stated we neglect other contributions in the quark-antiquark scattering kernel as we expect the Goldstone modes to be dominant. The corresponding expression for the t-channel should be related by crossing symmetry (if we take incoming and outgoing particles on-shell). As the latter will however be broken when only considering the leading term in Eq. (10), a transformation of the leading term in Eq. (10) to the t-channel is not unique and depends on the choice for incoming and outgoing momentum. In order to find a suited kernel for the BSE which fulfills the axWTI in Eq. (8) later on, we take the arithmetic mean of the two possible contributions and obtain for the self-energy with k = p − q and D π (k) = (k 2 + M 2 π ) −1 being the pion propagator. The factor 3 in front of the pion contribution stems from the flavor trace and represents contributions from π + , π − and π 0 . These are treated on equal footing in the isospin-symmetric limit adopted here. The general structure of the Bethe-Salpeter vertex of the pion can be represented by For the case H π = 0 it is then possible to respect the axWTI in Eq. (8) for P 2 → 0 in the chiral limit by choosing the kernel K pion tu;rs (q, p; P ) = − Neglecting H π would therefore be a suited approximation, especially since it is known to have very small impact on hadronic observables (although the vertex function is then, strictly speaking, no longer Poincare invariant). We will however follow another strategy here by symmetrizing our kernel in (q, p, P ) ↔ (−q, −p, P ) for which we obtain This is the only kernel we have found to analytically respect the axWTI for P 2 → 0 in the chiral limit for the general structure of the Bethe-Salpeter vertex given in Eq. (12). From a theoretical viewpoint it is less well grounded than the choice Eq. (13): whereas (13) can be represented by proper Feynman diagrams the symmetrized version Eq. (14) does not share this property. It is however motivated from a pragmatic point of view since it allows one to check the influence of H π explicitly. We will therefore prefer the version Eq. (13) in this work. As a result we will see that neglecting contributions from H π introduces errors on the level of a view percent into observables such as the pion decay constant. For many applications it may therefore be possible to use the simpler kernel (13). Both kernels have the correct charge conjugation properties and respect multiplicative renormalizability of the BSE by construction.
The resulting diagrammatic expression of the Bethe-Salpeter equation is shown in Fig. 6. The sum over all contribution in Eq. (14) is presented by the one-pion exchange diagram. Our approximation of the quark-DSE in Fig. (3) by the one given in Fig. (5) is justified provided the Bethe-Salpeter vertex functions of the pion with or without the pion interaction term do not differ strongly. We will show in section III B that this is indeed the case.
With the kernel (14) our final expression for the quark-SDE is then given by with k = p − q. The corresponding expression for the BSE reads with the kernels K Y M tu;rs and K pion tu;rs given in Eqs. (9) and (14). The Bethe-Salpeter vertex of the pion, i.e. the dressing functions E π , F π , G π and H π are normalized according to the condition with χ j π (k, P ) = S(q + P/2)Γ j π (k, P )S(q − P/2). The condition is written for the momentum partitioning ξ = 1/2 without loss of generality. The trace is over Dirac matrices and the conjugate vertex function Γ is defined as with the charge conjugation matrix C = −γ 2 γ 4 . The normalization condition (17) guarantees that the corresponding residue of the pion pole in the four-point quark-antiquark Green's function is unity [29]. The leptonic decay constant characterizing the pion coupling to the point axial field is then given by [30] where again the trace is over Dirac matrices. The numerical details needed for solving Eqs. (15) and (16) simultaneously are discussed in subsection III A. We wish to emphasize that the approximations made in the course of this section respect multiplicative renormalizability. This is shown explicitly in appendix A. Since also the axWTI is preserved we may hope that our approximation scheme captures essential properties of full QCD. In this context there are some interesting points to note: • A crudely simplified version of our interaction was considered by Gribov in his investigations of quark confinement due to supercritical charges [3]. He conjectured that the inclusion of the pion back reaction to the quark induces important changes in the analytic structure of the quark propagator. We will come back to this point in subsection III D.
• We observe the quark self energy in Eq. (15) to be infrared divergent in two dimensions for any non-vanishing Bethe-Salpeter amplitude in the chiral limit. More precisely the axWTI requires Γ i π (p, 0) = τ i γ 5 B(p)/f π [5] and we have in good approximation for d dimensions. For any finite value of the dressing function B(p 2 ) at any momentum p, the integral of the pion contribution is infrared singular for d = 2 and would shift the value B(p 2 ) → −∞, i.e. would give an infinite repulsive back reaction. Thus only the trivial solution B(p 2 ) = 0 is left in two dimensions. This supports Coleman's theorem [4], which states that continuous symmetries in two dimensions cannot be broken dynamically.
• Using the approximation Γ i π (p, P ) ≈ τ i γ 5 B(p)/f π [5] in the s-channel contribution of the pion to the quark-antiquark scattering kernel (see Eq. (10)) a simpler choice for K pion tu;rs (q, p; P ) could be We consider this as an approximation to our construction for the kernel.
• Finally we note that the inclusion of a back reaction by Goldstone bosons is possible for any case of dynamical symmetry breaking on an equal footing, e.g. it can similarly be applied to color-superconducting phases elaborated in [31] for our approach.
In the following subsection we will reconsider our approximation scheme from a different perspective.

C. The 1/Nc-expansion
The approximation specified in the last subsection may also be viewed in the light of the 1/N c -expansion. The pion exchange contributions are then the (approximate) next-to-leading order correction to the leading one-gluon exchange of the pure Yang-Mills interaction. Such a view has been advocated already in model frameworks such as the NJL-model [9,10,11]. In the following we discuss this aspect also for QCD in the framework of two-particle irreducible (2PI) effective actions, i.e. the CJT-formalism, which also allows for a derivation of the SDEs for propagators. The CJT action [32] for truncations with effective one-gluon exchange is a functional of the quark propagator and given by Here we already neglected the dependence on the expectation value of quark fields ψ and ψ as they are vanishing for the ground state.
For the diagrammatic expansion of the 2PI functional Γ 2 [S] we can then use a 1/N cordering. The 1/N c -expansion of vacuum-vacuum contributions needed for Γ 2 [S] can nicely be arranged by the topology of the contributing diagrams [33]. Since we do not consider pure gluonic contributions at order O(N 2 c ), the leading contribution is given by planar gluonic diagrams with a quark line as a boundary, being O(N c ). The corresponding contribution to Γ 2 [S] is the left diagram of Fig. 7. To next-to-leading order, i.e. O(1) , we have the topology of a cylinder with two quark lines as boundaries. Those contributions to Γ 2 [S] are subsumed in the right diagram of Fig. 7.
Therefore the truncated quark-SDE to next-to-leading order in a 1/N c -expansion, given by δΓ[S] δS = 0, turns out to be where Σ 0 (p) = Σ (rainbow) (p) is the rainbow self energy from Eqs. (2)(3)(4) and with M tu;rs (q, p; P ) = K Y M tu;rs (q, p; P ) is the next-to-leading order contribution.
Since we do not want to perform the ladder resummation Eq. (25), we furthermore want to approximate M tu;rs (q, p; P ) in order to reduce the numerical complexity. First we note that Eq. (24) is still correct up to next-to-leading order in 1/N c when using the propagators from the rainbow-ladder approximation. The matrix M tu;rs (q, p; P ) then contains Goldstone bosons in the s-channel. As we again presume the Goldstone bosons to be dominant, we will focus on the leading term in Eq. (10). Following our discussion [Γ a π ] tu ( q; P ) for the ladder resummed fully-amputated quark-antiquark scattering amplitude. This is again a pragmatic approach as we need to symmetrize the arguments in M tu;rs (q, p; P ) in order to find a manageable Bethe-Salpeter kernel. We then recover the quark-SDE in Eq. (15) as the approximate next-to-leading order expression within the CJT-formalism. Strictly speaking we would then need to use the rainbow-ladder Bethe-Salpeter vertex functions in the expression for the self-energy, however to be consistent up to next-toleading order in 1/N c we can also insert next-to leading order expressions. A simpler truncation would be to use the approximation given in Eq. (21).
Using our approximated next-to-leading order contribution of the self energy and neglecting the dependence of the Bethe-Salpeter amplitudes on the propagators then allows to justify the kernel given in Eq. (14). Essentially the kernel guarantees the pion to be massless in the chiral limit. Taking the 1/N c -expansion more serious by performing the ladder summation results in a more complicated expression for the Bethe-Salpeter kernel. We need to evaluate the derivative in Eq. (27). Also including the momentum in the multi-indices we abbreviate Eq.
We then get for the next-to-leading order contribution to the Bethe-Salpeter kernel. The leading order contribution is simply given by the effective one-gluon exchange. Evaluating the derivative by using Eq. (25) Summarizing this analysis we obtain the next-to-leading order contribution to the Bethe-Salpeter kernel to be of the form which corresponds to a resummed ladder in the t-channel and two contributions with two resummed ladders in the s-channel. The contributing diagrams are actually similar to those introduced in [12]. Having discussed these terms we can think of Eq. (26) as an approximation of the first term of Eq. (32), which is identical to the one introduced in subsection II B. A diagrammatic presentation of Eq. (32) is shown in Fig. 8.

D. Gluonic interaction
As discussed in previous subsections we need to specify two different components of the Yang-Mills part of the interaction: the gluon dressing function Z(p 2 ) and the vertex dressing Γ Y M . The gluon dressing function has been calculated numerically from a truncated version of the Schwinger-Dyson equations for the ghost and gluon propagator [13,35]. The resulting dressing functions for the ghost and gluon propagator have been discussed and compared with corresponding lattice results in previous works, see e.g. [36] and references therein. In order to make this paper self-contained we use functional forms for the ghost and gluon dressing functions which represent the numerical solutions to sufficient accuracy [23]. These forms are given by with the running coupling Unquenching effects due to quark-loops in the gluon polarization are included in the expressions (33).
We show (33) together with lattice results for the unquenched (N f = 2 + 1) gluon dressing function [39] and the quenched 1 ghost dressing function [40] in Fig. 9. For the gluon both results are in good agreement in the infrared momentum region (deviations are discussed in [36]). In the intermediate momentum region truncation artifacts in the gluon-SDE lead to an underestimation of the hump in the dressing function. In the ultraviolet the SDE-solutions show the correct logarithmic scaling as expected from resummed perturbation theory, whereas the lattice data suffer from artifacts due to the finite lattice spacing. The ghost dressing function show slight deviations in the infrared, which are discussed in detail in [36]. In general we believe that both solutions are in sufficient agreement such that the expressions (33) serve as a reliable input for the quark SDE.
For the second part of the Yang-Mills interaction, the vertex dressing Γ Y M , we use a procedure first suggested in [25] and further explored in [26]: we devise a functional form for the vertex dressing with parameters fixed by the requirement that lattice data for the quark propagator are reproduced by solutions of the quark-SDE using the complete interaction. The functional form we employ for this procedure is given by with the components where τ = e − 1 acts as a convenient infrared cutoff for the logarithms. As mentioned in section II A other tensor structures in the Yang-Mills part of the quark-gluon vertex might be relevant, an investigation of those is however beyond the scope of this work. We will come back to this point at the end of section II D. In the following we discuss the choice of our ansatz (35) step by step and compare it to the one used in the quenched calculation of ref. [26]. It is well known that the effective interaction in the quark-SDE has to approach the running coupling in the ultraviolet momentum regime [41], i.e.
Up to constants this UV-behavior is represented by Γ 1 with the scheme-dependent scale Λ QCD . The anomalous dimension γ m of the quark propagator is given by γ m = 12/(33 − 2N f ), where we use N f = 2 in our calculations. Furthermore note that Z 2 = Z 3 Z 1F in Landau gauge, where Z 3 is the renormalization factor of the ghost fields. The product Z(k 2 )Γ 2 (k 2 )Γ 3 (k 2 ) goes to a constant for large momenta; the scale Λ Y M is determined from the SDE-solutions for the ghost and gluon propagator (33). The coefficient h is fixed such that the limit (39) is satisfied. The renormalization group invariant G(ζ 2 , µ 2 ) Z 3 (µ 2 ) with the arbitrary but fixed scale ζ is introduced to impose the correct cutoff-and renormalization point dependencies of the effective interaction in the quark-SDE. In our calculations we use ζ = 2.9 GeV, other choices do not affect the result.
The infrared behavior of our interaction is dominated by Γ 2 , which is diverging for p 2 → 0 whereas Γ 1 and Γ 3 go to constants. The degree of divergence is given by This behavior has been derived in ref. [24] from an analytic, selfconsistent analysis of the full SDE for the quark-gluon vertex given in Fig. 1. We wish to emphasize, however, that the precise infrared behavior of Γ Y M is not important for the quark-SDE: this equation is dominated by intermediate loop momenta. One could equally well work with lower degrees of divergence [26] or even with an infrared finite vertex [25].
Finally, from the Slavnov-Taylor identity of the quark-gluon vertex one can see that the vertex is also proportional to the quark wave function Z f . This dependence is taken care of by Γ 3 , which mimics the momentum dependence of Z f . The extra factor Z 2 is vital for ensuring multiplicative renormalizability of the quark-SDE. The dependence of this part of the vertex on the quark mass is expressed in terms of the function  In order to preserve multiplicative renormalizability of the quark-SDE, it is important that the scale ζ is not correlated with the renormalization point. Instead it should be a fixed scale. In our calculations we use the same scale ζ = 2.9 GeV as also present in Γ 2 .
To fit the parameters of the Yang-Mills part of our interaction we use the following procedure: We first determine the values of h, Λ QCD and a(M) at a fixed current quark mass such that (i) the correct ultraviolet behavior (39) of the running coupling is reproduced and (ii) the unquenched lattice quark propagators from ref. [15] are reproduced. The Yang-Mills scale Λ Y M is taken from the SDE-results parameterized in Eqs. (33). This procedure is repeated for the four different current quark masses of ref. [15], which correspond to M(2.9 GeV) = 44, 65, 85, 106 MeV. We then fit Eq. (40) to the results for a(M) and determine a 1 , a 2 and a 3 . The results for the parameters are given in Tab. I.
The numerical results of our fitting procedure are shown in Fig. 10. For the quark mass function M(p 2 ) we find excellent agreement with the lattice data in both, the infrared and ultraviolet momentum regions. For the wave function Z f (p 2 ) the agreement is only slightly less convincing. Particularly interesting is the small spread of the wave functions for different quark masses, which is reproduced by the SDE results. This is in contrast to earlier investigations of the quenched quark propagator, where the spread in the SDE results has been too large [25,26]. We consider the improved result here as an indication that our quark-gluon interaction is more realistic than in previously used models. In particular we attribute this to the fact that the hadronic contributions in the quark-gluon interaction generate all twelve possible Dirac structures of this interaction (only eight of these are independent in Landau gauge), whereas the quenched investigations [25,26] mainly worked with the γ µ -structure only. This is a clear hint that one should include further structures in the modeling of the Yang-Mills part Γ Y M of the interaction and is left for future work. Our solutions for the quark mass function M(p 2 ) and the wave function Z f (p 2 ) for physical up/down quarks are discussed in section III B. Before presenting results we need to detail our numerical procedures. Our choice of the momentum routing in the SDE and BSE is specified in Fig. 11. In the right diagram P denotes the total momentum of the meson and q + = q + P/2, q − = q − P/2 and p + = p + P/2, p − = p − P/2 are the momenta of the internal and external quarks. Since in the rest frame of the meson we have P µ = (0, 0, 0, iM), where M is the mass of the bound state, we need to know the quark propagators in the BSE for complex momenta. In addition, for the pion exchange diagram we need to know the pion Bethe-Salpeter amplitude for complex relative momenta. Whereas solving the quark-SDE in the complex plane is feasible and a standard procedure by now (see e.g. [14,42,43]), the solution of the BSE for complex momenta has not yet been performed to our knowledge. In view of the considerable numerical complexity of the coupled system of SDE and BSE we will not attempt such a solution within the context of this work. Instead we resort to the well explored 'absolute value approximation', which amounts to replacing all momentum arguments in internal quark propagators in the BSE by their absolute value. In addition, we replace the arguments of the pion amplitudes in the meson exchange diagrams by their real parts. This approximation has the merit of making the kernel independent of P µ , which in turn greatly simplifies the Bethe-Salpeter norm (17). We will attempt a complete solution of the problem in the complex plane in future work.
To estimate the quality of this approximation we employed the (quenched) model of [43] and solved it once for complex momenta and once in the absolute value approximation for a physical pion. Using similar current quark masses the resulting difference in the pion mass is less than 1 MeV. The error in f π however is larger: we obtain f π = 94 MeV with complex momenta but only f π = 82 MeV in the absolute value approximation. Roughly 10 MeV of this difference is due to the absolute value approximation in the norm integral (17); the remaining 2 MeV are accounted for by shape distortions in the F , G and H amplitudes of the pion Bethe-Salpeter amplitude (12). We expect similar errors when using the quenched version of our interaction, i.e. only taking into account the gluon exchange part of the interaction given by Eq. (35). In the unquenched calculations the error in the norm integral is backfeeding into the coupled system by the overall normalization of the pion contribution of the interaction. We therefore expect somewhat larger errors in both, the pion mass and the resulting pion decay constant. Nevertheless we consider our results to be accurate and meaningful in their qualitative features.
Further details concerning the numerical techniques needed to solve the quark-SDE are given in [13,14]; corresponding techniques for the BSE are described in [5,14].

B. Quark propagator and pion properties
Our results for the physical up/down quark propagator are shown in Fig. 12. We compare the unquenched (N f = 2) solutions from the SDE/BSE system with partially quenched solutions. Here, 'partial quenching' means that we still include the quark loop effects in the gluon propagator, implicitly present in Eqs. (33), but switch off the hadronic unquenching effects from the pion part of the interaction. In both calculations the renormalized current quark mass is chosen such that the pion mass is M π = 136 MeV. In an MS-scheme with Λ M S = 225 MeV this corresponds to renormalized up/down current quark masses of Note that the unquenched current quark mass is larger than the quenched one. In [14] a small effect in the opposite direction has been found from quark loops in the gluon polarization alone. The stronger increase found here reverses this tendency and indicates that unquenching effects in the quark-gluon vertex are stronger than those in the gluon polarization.
The quark mass function in the upper left diagram of Fig. 12 also shows a sizable reduction of the dynamically generated mass M(p 2 ) in the infrared when the theory is  unquenched. These effect is greatly diminished at p = 2 GeV and vanishes for momenta larger than p = 3 GeV. The corresponding effect in the quark wave function Z f (p 2 ) is much smaller and only visible for momenta lower than p = 0.5 GeV. This finding agrees qualitatively with results from quenched and unquenched (N f = 2 + 1) lattice calculations shown in the lower panel of Fig. 12. Clearly, both approaches underpin the physical intuition: in general, quark loop effects are screening and compensate part of the antiscreening gluonic interaction. This leads to a reduction of the interaction strength in the gluon propagator [13,39]. Here it also reduces the strength of the quark-gluon vertex and consequently entails (further) reduced quark masses. The effects for the quark mass function seen in the SDE/BSE approach are larger than in the lattice calculations. This may or may not be an artifact of the absolute value approximation used in the numerics. In any case, it seems save to conclude that unquenching effects from the pion back reaction (i.e. in the quark-gluon vertex) are considerably stronger than the quark-loop effects in the gluon polarization investigated in refs. [13,14]. The corresponding Bethe-Salpeter amplitudes of the pion are shown in Fig. 13. Similar to the previous findings in the quenched approximation [42] we find a hierarchy in magnitude of these amplitudes: clearly E π is most important in size followed by F π , G π and H π . This hierarchy is not changed by unquenching, although we note a considerable increase in size of G π and H π for N f = 2 compared to the quenched case. Despite these differences both sets of amplitudes are similar enough to justify the approximation made in subsection II B. In the infrared, E π , F π , G π and H π go to constants. In the ultraviolet E π and F π fall like 1/p 2 times a logarithm, whereas G π and H π are proportional to 1/p 4 times a logarithm. This behavior is not changed when the system is unquenched. The relative importance of the Bethe-Salpeter amplitudes can be inferred from table II. We show results for the quark mass function at zero momentum, the pion mass and the contributions to the pion decay constant when part of the amplitudes are omitted. On a quantitative level it is certainly important to include all four amplitudes in the calculation. If, however, one is satisfied with qualitative results, even a calculation involving  only the E π -amplitude would give satisfactory results. Using all four Bethe-Salpeter amplitudes, the resulting values for the pion decay constant calculated from (19) in the absolute value approximation and for a physical value of the current quark masses are Clearly, unquenching reduces the value of f π . However, the magnitude of this reduction is too large. Also both values are too low. One would expect the unquenched value close to the physical point f exp π = 93 MeV and a somewhat larger number in the quenched case. A reduction of the unquenched pion decay constant was also found in recent lattice calculations with quenched [44] and unquenched [45] twisted mass fermions. They obtained for our definition of the pion decay constant and in the chiral limit. In the unquenched simulation the pion decay constant at the physical point was used to determine the physical scale. The uncertainty of this scale, at least in quenched calculations, is usually expected to be on the order of 10%. As discussed in section III A at least part of the deviations of (42) from (43) can be attributed to the absolute value approximation. However, it is by no means clear, whether a calculation in the complex plane will bring f π all the way up to the physical value. An indication that this may well not be the case is a comparison with f π from a modified Pagels-Stokar approximation [46,47] This approximation is not affected by the absolute value problem, but in turn incorporates only the effects of the leading pion Bethe-Salpeter amplitude E π in the chiral limit. Thus the systematic errors are different. By comparison of (44) with the full result in model calculations [42,43] one finds that the approximation (44) should lead to an underestimation of the full f π of the order of 10%. In our case we obtain Compared to (42) we find a considerable increase in the unquenched value and no change for the quenched case. Within our quenched interaction we can also perform a full calculation in the complex plane using the techniques described in [14]. We then find the value which is indeed roughly 10% larger than the result (42) of the absolute value approximation and also 10% larger than the result (45) from the Pagels-Stokar approximation. From this one would expect to find a value of roughly f N f =2 π,complex ≈ 70 − 75 MeV from an unquenched calculation in the complex plane. Clearly both, the calculated result (46) for the quenched case and the conjectured result for the unquenched case fall short of (43) by almost twenty percent.
A possible source of these shortcomings are errors in the overall scale of our calculation. This scale was set by comparison with the lattice results for the gluon [39] and quark propagator [15]. Indeed, already in [25] and again in [26] it was found that the pion decay constant calculated from an interaction fitted to quenched lattice data yields an underestimation of f π by about thirty percent. This agrees with our findings above and supports the argument that there may be problems with the lattice scale. We also note that a rescaling of the lattice units would lead to larger constituent quark masses and in addition reduce the effect of the Goldstone contribution in the selfenergy. This would improve the agreement of the results shown in Fig. (12).

C. Gell-Mann-Oakes-Renner relation and chiral condensate
The Gell-Mann-Oakes-Renner relation can be derived from the inhomogeneous Bethe-Salpeter equations for the axialvector and pseudoscalar currents and the axialvector Ward-Takahashi identity [7]. Since the left-hand side of this equation only contains physical observables, the dependence of the renormalization point of the current quark masses m u (µ) and m d (µ) and the chiral condensate ΨΨ 0 (µ) have to cancel each other. The individual dependence of the pion mass and the pion decay constant on the current quark mass is shown in the upper panel of Fig. 14.
Clearly, one observes the characteristic square-root behavior of the pion mass that becomes massless in the limit m u/d → 0. No qualitative changes can be seen by naked eye between the quenched and unquenched calculation. The situation is somewhat different for the pion decay constant. Here an almost linear behavior for sizable quark masses can be seen, with small deviations towards the chiral limit. Whether the different curvature for m u/d → 0 can be explained by the presence of chiral logs remains to be explored in a more refined calculation. In the lower panel of Fig. 14 we plot M 2 π as a function of the current quark mass. Clearly, the dependence is linear for a wide range of current quark masses. The deviations from the linear fit for small quark masses may well be an artifact of the increasing numerical uncertainties close to the chiral limit. The deviations for larger quark masses show the onset of quantitative corrections to the Gell-Mann-Oakes-Renner relation. It is interesting to note, however, that these corrections play no role up to pion masses of M π = 400 MeV. A similar result has been found in a recent lattice study using Wilson quarks and fine lattices [48].  for the quenched calculation using complex momenta. The error subsumes numerical errors as well as errors in fitting M 2 π (m), which has been performed by a least squares fit in the region from m(2 GeV) = 3 − 40 MeV. By comparison with (46) we note that the GMOR-relation is satisfied at the 4 percent level for the central value in (51). The difference of the corresponding value for the unquenched case and the result Eq. (42) shows once more that the absolute value approximation leads to a large error in the unquenched f π .

D. Analytical properties of the quark propagator
Finally we take a look at the analytic structure of the quark propagator in the complex momentum plane. We wish to stress that the results of this section have to be treated with caution, since the absolute value approximation is likely not to be reliable in this respect. Nevertheless, it may be instructive to have a first glance in the complex plane. To this end we determine the Schwinger-function where σ S (p 2 ) = B(p 2 )/(p 2 A(p 2 ) 2 + B(p 2 ) 2 ) is the scalar part of the dressed quark propagator. (This method has a long history, see e.g. [23] and references therein). According to the Osterwalder-Schrader axioms of Euclidean field theory [49], this function has to be positive to allow for asymptotic quark states in the physical sector of the state space of QCD. Conversely, positivity violations in the Schwinger function show that the corresponding asymptotic states (if present) belong to the unphysical part of the state space. In Fig. 15 we show results for the Schwinger function of the quenched and unquenched quark propagator for the physical up/down quark of Fig. 12. Both functions are qualitatively similar. An excellent fit to the Schwinger function is obtained using the form [50] |σ S (t)| = | b 0 exp(−b 1 t) cos(b 2 t + b 3 ) | , which corresponds to a pair of complex conjugate poles of the propagator in the time-like momentum plane. These poles correspond to a 'quark mass' given by m = b 1 ± ib 2 . In our quenched case this amounts to m = 505(10) ± i 360(10) MeV; the unquenched quark mass is m = 351(10) ± i 254 (10) MeV. The mismatch of the fit with σ(t) at small times shows the presence of additional analytic structure either in form of a cut along the real time-like momentum axis or in form of additional singularities further away from the origin. If this behavior is stable also when the absolute value approximation is overcome it contradicts the Gribov scenario given in [3]. However, as emphasized already above, definite conclusions can only be drawn from a complete calculation including the full complex momentum structure involved in the coupled system of the quark-SDE and the pion-BSE.

IV. SUMMARY AND OUTLOOK
In this work we investigated the pion back reaction on the quark propagator and resulting pion properties. To the best of our knowledge this is the first calculation of this kind in a non-perturbative continuum approach to QCD. We isolated the pion contribution to the quark-gluon vertex and identified an approximation to the quark-SDE which allows for the construction of a Bethe-Salpeter kernel in agreement with the axial Ward-Takahashi identity. This setup is powerful enough to analytically verify Coleman's theorem in 1+1 dimensions: chiral symmetry cannot be broken spontaneously there. Our main interest, however, is QCD 3+1 , where we evaluated the back reaction effects numerically. Here we had to resort to an absolute value approximation, which is reliable as concerns the pion mass but problematic for the pion decay constant.
As a result we have obtained considerable unquenching effects in the quark mass function at small momenta. The screening effect of the pion interaction reduces the quark mass M(0) from M N f =0 (0) = 477 MeV to M N f =0 (0) = 371 MeV. This effect becomes smaller for larger momenta and vanishes above p = 3 GeV in correspondence with the width of the Bethe-Salpeter amplitude of the pion. Almost no corresponding effect in the quark wave function is seen.
The screening effect of the pion interaction also reduces the value of the chiral condensate. The Gell-Mann-Oakes-Renner relation is satisfied with an accuracy less than 4 %. We obtain an almost linear relation between the squared pion mass and the current quark mass for pion masses up to 400 MeV with only slight deviations up to 500 MeV. An unsolved problem is posed by our low value of the pion decay constant. Certainly, part of the problem is the absolute value approximation mainly used in this work. However, we also argued that the overall scale of our calculation, which is based on the lattice results of [15,39], may have to be corrected by about 20 %.
Our findings are certainly reliable on a qualitative basis. To also obtain reliable quantitative results we need to go beyond the absolute value approximation for the arguments in the internal quark propagators and pion amplitudes. Work in this direction is well under way. We then hope to better understand the problem of the very low pion decay constant obtained in this work. Also, working in the complex momentum plane is mandatory to investigate possible unquenching effects in the analytic structure of the quark propagator conjectured by Gribov. Finally, a complex treatment of the system is mandatory to explore decay thresholds in vector, axialvector and the scalar meson sectors. An investigation of this mesons should be simplified by use of the Maximum Entropy