Bethe-Salpeter studies of mesons beyond rainbow-ladder

We investigate the masses of light mesons from a coupled system of Dyson-Schwinger and Bethe-Salpeter equations. The dominant non-Abelian and sub-leading Abelian contributions to the dressed quark-gluon vertex are explicitly taken into account. We also include unquenching effects in the form of hadronic resonance contributions via the back-reaction of pions. We construct the corresponding Bethe--Salpeter kernel that satisfies the axial-vector Ward-Takahashi identity. Our numerical treatment fully includes all momentum dependencies with all equations solved completely in the complex plane. This approach goes well beyond the rainbow-ladder approximation and permits us to investigate the relative impact of different corrections beyond rainbow-ladder on the properties of mesons. We find that sub-leading Abelian corrections are further suppressed dynamically, and that our results supersede early qualitative predictions with significantly simpler truncation schemes.


Introduction
In hard scattering experiments we can probe the structure of mesons and baryons, thus revealing their partonic building blocks; the quarks and gluons. This picture is natural at high-energies where we have asymptotic freedom and mirrors the apparent degrees of freedom in the QCD Lagrangian. However, at low to intermediate energies where the effects of confinement become more important, a description in terms of hadronic degrees of freedom is more appropriate. In the description of bound-states, in order to study their properties and interactions in terms of their constituent particles, we require non-perturbative tools that are capable of relating both the small and large scales together. This is necessary to provide, for example, a connection between the symmetries inherent in the QCD Lagrangian and their manifestation in the hadronic spectrum.
Suitable tools include chiral perturbation theory, quark models and of course lattice QCD. Though heavy quark systems can be adequately described with non-relativistic models, and lattice calculations are quickly approaching the physical point, we are also interested in how the light meson spectrum arises from dynamical chiral symmetry breaking. This suggests that we use an approach in which we can make a direct connection to the limit of vanishing quark mass; moreover, we still wish to employ a description in terms of quarks and gluons rather than effective degrees of freedom. One such approach that is intrinsically non-perturbative, formulated in the continuum and capable of exploring both chiral and heavy quarks are the Dyson-Schwinger equations (DSEs) and Bethe-Salpeter equations (BSEs), which describe one-particle irreducible Green's functions and two-body bound states.
The DSE-BSE framework has been explored extensively in the literature with regards to the light-meson a e-mail: richard.williams@physik.tu-darmstadt.de sector over the last few decades, for reviews see Refs. [1,2]. However, the vast majority of studies have been undertaken under the umbrella of the rainbow-ladder (RL) approximation, in combination with a phenomenological approach to the quark-gluon interaction. The main reason for the focus on this truncation is clear: firstly, it satisfies the vector and flavour non-singlet axial-vector Ward-Takahashi identities (axWTI) so necessary for the identification of the pion as a (pseudo)Goldstone boson; and secondly because of its sheer simplicity. The last few years have seen these same studies extended to the description of baryons via the quark-diquark approximation and more recently through the full covariant three-body Faddeev equations [3,4,5,6].
With the successes of RL in the meson sector [7,8,9,10,11], and also for the baryons, it is prudent to attempt to improve upon our approximation schemes. The most important object of interest here is the quark-gluon vertex since its specification ultimately determines the interaction between our fundamental degrees of freedom. An extension of the truncation scheme can give rise to Yang-Mills corrections, unquenching effects and eventually to coupled channel processes and decays; all of these are desirable properties that we wish to include in a more complete description of mesonic bound-states. Since a well-known prescription exists for constructing a DSE/BSE truncation based on diagrammatic expansions of the vertex [12,13], this has been the main focus over the last few years. Because of the sheer complexity of such studies which require higher loop integrals to be computed, most studies [13,14,15,16,17] employ the simplistic Munczek-Nemirovsky delta-function ansatz [18] for at least one of the exchange gluons. This kills a loop integral making studies tractable, and in some cases analytical, at the expense of barring the dominant vertex corrections from consideration. In some cases this also restricts the study to just the pseudoscalars and vectors with other mesons becoming unbound. We will later see that the use of such a non-propagating ansatz gives rise to truncation artefacts that are hard if not impossible to disentangle from the results, invalidating some of the qualitative predictions such investigations give rise too. Regardless of these concerns, past studies are pioneering in themselves when one considers the resources available at the time with much being learned as a result.
While the investigation of Bethe-Salpter equations is important, we must also explore and improve our understanding and truncations of the Dyson-Schwinger equations that provide their input. For example, the Yang-Mills (YM) sector of QCD has also been much explored, resulting in solutions of the ghost and gluon propagators that are in qualitative agreement with lattice calculations, see Refs. [1,19]. These achievements encourage us to steer away from phenomenological studies and attempt to construct higher n-point Green's functions using all of the information available [20,21,22]. Recently, Dyson-Schwinger studies of the quark-gluon vertex and quark propagator have been made in which the dominant non-Abelian corrections and sub-leading Abelian contributions were considered for the first time [23]. In lieu of these developments, and following our programme of investigating Bethe-Salpeter equations beyond rainbow-ladder [24,25], a significant extension of the Bethe-Salpeter equation beyond the rainbow-ladder truncation has been proposed [26]. There, the dominant non-Abelian corrections to the quark-gluon vertex are considered in a symmetry preservation truncation of the Bethe-Salpeter equation, with the corresponding light-meson spectrum calculated.
Because of the technical challenges in such a state-ofthe-art calculation the study was exploratory, relying upon phenomenological input for the gluon propagator. Though this is undesirable in light of what we know about the basic Green's functions of QCD, this does not trivialize the calculation itself. The integrals in Euclidean space remain two-loop in nature, and require our input quantities to be evaluated for complex momenta. The aim of the present paper is not to take the next step and introduce a quantitatively reliable truncation scheme. Rather, our intention is provide further corrections to the quark-gluon vertex in addition to those considered thus far [26]. Hence, we will include both the sub-leading Abelian corrections and the pion back-reaction that represents hadronic unquenching effects. As a result, we obtain a model that is orders of magnitude more complicated than the monumental step already made. Thus, for the time-being, we satisfy ourselves with a qualitative study of corrections beyond rainbow-ladder and their relative impact on mesonic bound states. The extension to a consistent and qualified truncation scheme will be reported in detail elsewhere. This paper is organised as follows: in section 2 we introduce the Dyson-Schwinger equation for the quark-propagator and the quark-gluon vertex. The truncation scheme is introduced, being sub-divided into a non-resonant Yang-Mills part and a resonant hadronic contribution. In section 3 we introduce the Bethe-Salpeter equation for a quark-antiquark bound-state, detailing the relativistic decomposition of the bound-state amplitude, general solution methods, their normalisation and the interaction kernel. In section 4 we give the results of our calculation and compare to existing literature. Finally, we conclude and provide an outlook in section 5.

Quark propagator
The quark propagator and the gap equation it satisfies is one of the most important quantities in the covariant description of mesons. It encodes such non-perturbative properties as dynamical mass generation and the realisation of a non-zero condensate through the dynamical breaking of chiral symmetry. Moreover, as we will see in section 3.3 a symmetry preserving truncation of the Bethe-Salpeter kernel can be constructed directly from the self-energy part of the quark Dyson-Schwinger equation.
In Euclidean momentum space, the renormalised dressed gluon and quark propagators in the Landau gauge are given by where Z(p 2 ; µ 2 ) is the gluon dressing function, Z f (p 2 ; µ 2 ) is the quark wave-function and M(p 2 ) is the renormalisation point independent quark mass function. The dependence of such functions on the renormalisation point µ 2 will be implicitly assumed from here on. The quark dressing functions A(p 2 ) and B(p 2 ) can be recombined into the quark mass and wave-function by M(p 2 ) = B(p 2 )/A(p 2 ) and Z f (p 2 ) = 1/A(p 2 ). These propagators may be obtained by solving their respective Dyson-Schwinger equations. The DSE for the quark propagator, shown diagrammatically in Fig. 1, is written where Σ(p) is the quark self-energy, k − p − q and the Casimir C F = 4/3 stems from the colour trace. We introduced the reduced quark-gluon vertex Γ ν (q, p) defined by Γ a ν (q, p) = ig λ a 2 Γ ν (q, p). The bare inverse quark propagator is S −1 0 (p) = ip + m. The renormalisation factors are Z 1F = Z 2 / Z 3 for the quark-gluon vertex, Z 2 for the quark propagator and Z 3 for the ghost dressing function.
The scalar dressing functions of the quark DSE are solved for by appropriate projections of Eq. (3). This is a coupled non-linear integral equation that is solvable provided we know the gluon dressing function and the structure of the quark-gluon vertex. In the rainbow approximation both are specified by Ansätze. Here, the quark-gluon vertex will be provided by solving its respective DSE in a truncation scheme to be introduced in the next section. Following the lead of past studies and our own investigations we employ the momentum dependent ansatz [27] which is a Gaussian distribution whose scale and strength of the effective gluon interaction are provided by D = 16 and ω = 0.5, respectively. Since it cuts out the dynamics at large momenta all of the equations are finite and do not require renormalisation. Such factors may then be set to one, though for reference we still show them explicitly. The use of such an interaction with this simple (but non-trivial) form simplifies some aspects of the calculation which will remain somewhat involved.

Quark-Gluon vertex
Though the quark propagator is one of the most important inputs to our Bethe-Salpeter equation, the most important quantity with regards to the dynamics of the theory is the quark-gluon vertex. This describes the non-perturbative interplay between our dressed quarks and gluons. We are left with two choices as to how to specify the details of this interaction: either we provide some physically motivated Ansatz for the quark-gluon vertex; or we attempt to solve its corresponding Dyson-Schwinger equation, shown in full in Fig. 2(a). Since a priori we don't know how the quark-gluon behaves non-perturbatively, the provision of an ansatz is challenging. The most well-studied QFT in the literature with regards to vertex functions and the construction of ansätze is QED; they usually employ restrictions from the Ward-Takahashi identities and multiplicative renormalisability [28,29,30]. However, since QCD is by no means an Abelian theory it is not clear how much information can be directly transferred between gauge theories [31,32,33]. One thing we do know for sure is how to decompose our vertex into an appropriate basis. In general, the quarkgluon vertex can be written as twelve independent tensors built-up of two independent four-momenta p µ 1 , p µ 2 and the Dirac matrices ½, γ µ . An oft used basis has been provided for us by Ball and Chiu [28]. The twelve components are split into longitudinal and transverse parts where the longitudinal and transverse basis components L µ , T µ in Euclidean space are given in Table 1. The scalar Table 1. Longitudinal and transverse components of the quarkgluon vertex in the Ball-Chiu basis. Here, the incoming gluon momentum is p 3 = p 2 − p 1 with p 1 , p 2 the incoming and outgoing quark momentum. We use σ µν = 1 2 (γ µ γ ν − γ ν γ µ ).
coefficients that parametrise this vertex, λ i and τ i , are in general functions of p 2 1 , p 2 2 , and p 2 3 . We note that this basis is only favoured for comparative purposes. Since we work in Landau gauge where the gluon propagator is explicitly transverse (and always accompanies the quark-gluon vertex), a reduced basis of just eight components is necessary. However, the choice of basis is irrelevant to the final solution though a different basis can provide significant simplifications to the calculation.
Rather than attempting to provide an ansatz for the quark-gluon vertex, we choose a more systematic approach and attempt to solve its Dyson-Schwinger equation. Since this is an intricate and highly non-trivial coupled integral equation, containing basically unknown four-and five-point functions, we must impose a truncation scheme. Thus, we approximate the DSE following the detailed investigation of Refs. [34,23]. This results in the approximate DSE portrayed in Fig. 2(b) in which two-loop diagrams have also been neglected. Here, the first non-Abelian diagram subsumes the first two diagrams in the full DSE to first order in a dressed skeleton expansion of the four-point functions. The quark-antiquark kernel in the third loop diagram of Fig. 2(a) is expanded in terms of resonance contributions to the kernel and one-particle irreducible Green's functions [34]. This gives rise to the (a) second Abelian diagram of Fig. 2(b) and one that involves pion exchange. In Fig. 2(b) we have shown explicitly the colour factors of the non-Abelian and Abelian diagrams. While the non-Abelian diagram is associated with colour factors N c /2, the Abelian diagram is N 2 c suppressed with a colour factor of −2/N c . Similarly, the resonant contribution is also proportional to 1/N c due to the implicit 1/ √ N c dependence of the two pion amplitudes. This gives us an indication as to the relative strength of the contributions from each diagram. It is clear that in the approximation scheme considered here all corrections are sub-dominant compared to the diagram containing the three-gluon vertex.
In Ref. [26] we discarded the sub-leading corrections, keeping only the non-Abelian diagram. Here we provide an extension by including the sub-dominant diagrams and investigate their contributions in turn. For distinction we will classify our diagrams as follows: the non-resonant non-Abelian and Abelian diagrams will be referred to collectively as the Yang-Mills part; the remaining resonant diagram corresponding to unquenching effects will be addressed as the pion back-coupling. In other words, we separate the hadronic resonance contributions from the other corrections to the quark-gluon vertex. We introduce each of these in turn in the following.

Yang-Mills contribution
The corrections beyond rainbow-ladder that correspond to non-resonant corrections will now be considered. In order to keep our calculation tractable we employ the wellestablished strategy of absorbing all internal vertex dressings of the diagram into effective dressing functions for the two internal gluon propagators. The resulting Dyson-Schwinger equation for the quark-gluon vertex Γ µ (p 1 , p 2 ) with quark momenta p 1 and p 2 and gluon momentum p 3 reads 4 , the renormalisation factors Z 1F ,Z 1 and Γ 3g the bare three-gluon vertex.
For the gluon, we employ the momentum dependent ansatz Eq. (4). Naturally, such an ansatz provides only a first step towards a full calculation of the Abelian and non-Abelian diagrams, including input from the DSEs for the three-gluon vertex and the gluon propagator. Nevertheless we believe that the ansatz of Eq. (4) is sufficient to provide for reliable qualitative results as concerns the relative effects of vertex corrections onto meson properties. In particular it is not sensitive to the question of scaling vs. decoupling [19] in the deep infrared, p < 50 MeV: both scaling and decoupling lead to a combination of three-gluon vertex and gluon propagator dressings that is vanishing in the infrared in qualitative agreement with the ansatz Eq. (4). In addition, quantitative effects in the interaction below p < 50 MeV are not expected to affect observables in the flavour non-singlet sector since the dynamical mass of the quark, M ≈ 350 MeV, suppresses all physics on scales p ≪ M (see however [35]). Finally, the ansatz Eq. (4) is not sensitive to details of the Slavnov-Taylor identity (STI) for the three-gluon vertex: exactly those longitudinal parts of the vertex that are constrained from the STI are projected out in any Landau gauge calculation by the attached transverse gluon propagators.
19th International IUPAP Conference on Few-Body Problems in Physics

Unquenching Effects
In the previous section we introduced the leading and subleading non-resonant corrections to the quark-gluon vertex. Here we consider the dominant hadronic resonant contributions corresponding to pion exchange. The contribution is further approximated according to [34,24] such that a symmetry preserving Bethe-Salpeter kernel is constructable. This results in the one-loop pion exchange contribution to the quark propagator, as shown in Fig. 3. The DSE for the quark, Eq. (3) is modified accordingly where D π (p) = 1/ p 2 + M 2 π is the pion propagator. The factor of three in the pion contribution comes from the flavour trace and represents contributions from π 0,± which we treat equally in the isospin limit. To further reduce the complexity of such a contribution, we employ the simple off-shell prescription for the pion used in [25] where we take only the leading component of the pion amplitude in the exchange diagram. Moreover, we replace it by its exact chiral limit value where B χ (p 2 ) is the scalar quark dressing function in the chiral limit to ensure the correct asymptotic behaviour, and f π is the leptonic decay constant. This ensures that the Bethe-Salpeter kernel has a non-trivial momentum dependence are produces results that are qualitatively in accord with phenomenology.

Solution
The system of Dyson-Schwinger equations for the quarkgluon vertex and quark propagator are solved by applying suitable projectors and taking the Dirac trace. This gives rise to fourteen coupled integral equations for the (scalar) basis coefficients shown in Eqs. (2) and (5) which we solve for on a discretised momentum and angular grid via gaussian quadrature. This can be accomplished very efficiently and without requiring large computing resources. For the Bethe-Salpeter equation it becomes necessary to evaluate the quark propagator, and hence the quark-gluon vertices at complex values of the Euclidean momenta. By choosing the momentum routing in an appropriate fashion, we can avoid the unconstrained analytic continuation of the gluon propagator and three-gluon vertex. To obtain solutions, we employ the 'shell-method' described in the appendix of Ref. [34].
In Fig. 4 we show the quark wave-function and massfunction for different truncations of the quark-gluon vertex: RL is rainbow-ladder; NA is the non-Abelian diagram; AB is the Abelian diagram; and PI represents the pion back-reaction. As can be seen, inclusion of the non-Abelian diagram leads to a much larger change in the mass function than the Abelian diagram alone. Taking both diagrams together it is obvious that the effect of the Abelian diagram is largely suppressed and is hidden underneath the NA curve (the difference in the mass-function is approximately 1 MeV at zero momenta). The only significant contributions arises when the pion back-reaction is also included, giving rise to the customary reduction of the quark mass-function. However, although the Abelian contribution appears suppressed here its impact on the Bethe-Salpeter equation must still be investigated since the dynamics associated with crossedladder kernels are often expected to be important.

Bethe-Salpeter equation
The Bethe-Salpeter equation describing a relativistic bound-state of mass M is calculated through Here Γ(p; P) ≡ Γ (µ) (p; P) is the Bethe-Salpeter vertex function of a quark-antiquark bound state, specified below. It is a homogeneous equation, with a discrete spectrum of solutions at momenta P 2 = −M 2 i corresponding to λ P 2 = 1. The lightest of these M i pertains to the ground state solution. The momenta k + = k + ηP and k − = k − (1 − η)P are such that the total momentum P of the meson is given by P = k + − k − and the relative momentum k = (k + + k − ) /2. The momentum partitioning parameter is η, of which physical observables are independent. For convenience we choose it to be 1/2 without loss of generalisation. The object K rs tu (p, k; P) is the Bethe-Salpeter kernel, whose Latin indices refer to colour, flavour and Dirac structure.
The Bethe-Salpeter vertex function Γ (µ) (p; P) can be decomposed into eight Lorentz and Dirac structures. The structure is constrained by the transformation properties under CPT of the meson we wish to describe [36]. In particular, our pseudoscalar, scalar and vector have quantum numbers J P of 0 − , 0 + and 1 − , respectively. The axial-vector can be parametrised in two different ways, dependent on its transformation under charge conjugation, J PC = 1 ++ and 1 +− . Taking this into consideration, we can write down a general basis for each desired meson amplitude, where the components T (µ) i are given in Table 2. For the axial-vector we point out that for the J PC = 1 +− (1 ++ ) state the first (last) four components significantly dominate and the remainder can be neglected. Also note that we can impose different charge conjugation properties on these covariants by changing the odd/evenness of the F i (p, P) functions i.e. by swapping filled for empty circles in Table 2. In this way, states with exotic quantum numbers may be constructed, see [37].
To give a definite example, the relativistic covariants required to describe a pion are as follows Note, as stipulated in the text, the extra factor of (p · P) with F 3 that imposes the even Chebyshev expansion for equal mass constituents. We solve the BSE for the meson amplitude Γ (µ) (p; P) via matrix methods using the following procedure. First, we project our BSE onto the scalar amplitudes of our meson decomposition. This gives rise to either four or eight coupled integral equations for the F i (p; P). To make manifest the angular dependence of the amplitude functions, we treat the total momentum P 2 as a parameter and expand the function as a series of Chebyshev polynomials in the angle p · P = p · P/|pP| The functions F k i (p 2 ; P 2 ) are projected out through use of the orthonormal properties of the Chebyshev polynomials. With the angular dependence made explicit, we can evaluate numerically the two non-trivial angles appearing in the integration measure. We cast the remaining radial integral in the form of a matrix equation by matching the external momenta to the radial loop momenta, p 2 j = k 2 j at the abscissae of our integration nodes. Thus the amplitude of our BS equation is projected onto the decomposition F k i p 2 j ; P 2 . Schematically we are solving for the column vector Γ as a parametric equation in P 2 = −M 2 with M the mass of the meson. A bound-state corresponds to solutions with λ = 1, which as an eigenvalue problem is equivalent to satisfying det (½ − K) = 0. However, since this equation is homogeneous the overall normalisation of the amplitude is not determined. This is obtained by imposing an auxiliary condition on the solution, discussed in the next section.

Decay constants and normalisation
The calculation of observables from the Bethe-Salpter amplitudes, such as leptonic decay constants for the pseudoscalar and vector mesons [7,8,9,10,11], requires that they are properly normalised. This normalisation condition is derived by demanding that the residue of the pole in the four-point quark-antiquark Green's function (from which the BSE is derived) is unity [38,39].
Using conventions such that f π = 93 MeV, the standard Leon-Cutkosky condition reads 3 where Q 2 = −M 2 is fixed to the on-shell meson mass, the trace is over Dirac matrices and the Bethe-Salpeter wavefunction χ is defined by χ j π (k; P) = S (k + P/2)Γ j π (k, P)S (k − P/2) .
The conjugate vertex functionΓ is given bȳ with the charge conjugation matrix C = −γ 2 γ 4 . For the rainbow-ladder truncation scheme, which we introduce below, the kernel is independent of the total momentum P and so vanishes under the derivative. For more sophisticated truncations this will not be the case; in Fig. 5 we show the diagrams that would arise in the truncation scheme that arises from our discussion of the quark-gluon vertex. These involve the evaluation of three-loop integrals over non-perturbative quantities. This can be satisfactorily tackled with the aid of standard Monte-Carlo techniques, with one caveat; the integration kernel for the crossed-ladder diagram is prohibitively large after naively taking the Dirac trace. Either we can spend some time massaging the expression to make it manageable -which still leaves a horrific calculation to be done -or we seek an alternative normalisation procedure. In fact, such an equivalent normalisation condition exists in the literature, courtesy of Nakanishi [40], and was recently explored in [26]. Using the momentum dependent eigenvalue λ(P 2 ) found from solving Eq. (9) we have the formula For comparative purposes, we tested this truncation against the non-trivial kernel's present in Refs. [16,25]. This is a one-loop expression and consequently requires significantly less numerical effort to evaluate. We find it simple to apply to all truncations solved via the homogeneous Bethe-Salpeter equation, such as the one presented here. The leptonic decay constant characterising the pion coupling to the point axial field is subsequently given by [39] where again the trace is over Dirac matrices, and k + = k + P/2, k − = k − P/2. There exist analogous expressions for the vector mesons [9], where an additional factor of 1/3 must be included due to summation over the polarisation tensor.
One may write a similar equation to (18), which corresponds to the residue of the pseudoscalar vertex: The axWTI imposes a relationship between these two residues, known as the generalised Gell-Mann-Oakes-Renner relation [7,8,41], which must hold at and beyond the chiral limit: where m u , m d are the masses of the up and down quarks at the renormalisation point µ 2 , and in this work considered to be degenerate. Confirming this relation serves as a check of both our numerics and indicates how well our truncation satisfies the axWTI. For the Yang-Mills part of our truncation this is exactly satisfied, whilst for our pion back-reaction contribution it is satisfied to better than 1% for pion masses of the order of 600 MeV.

Constructing a symmetry preserving truncation
In constructing a truncation scheme for our quark DSE and meson BSE, we must be careful not to break the symmetries of the theory. One of the most important of these is chiral symmetry. Satisfaction of this ensures the Goldstone nature of the flavour non-singlet psuedoscalars irrespective of the details of the model. This is an essential feature that any dynamical description of light mesons must reproduce, where explicit breaking due to the introduction of quark masses gives rise to the physical spectrum of mesons. Explicit and implicit mass dependences of our interaction should allow for an interpolation between light and heavy systems of quarks. The consequences of chiral symmetry can be expressed through the axial-vector Ward-Takahashi identity for flavour non-singlet mesons with p and P the relative and total momentum of the meson respectively, p ± is the combination p ± P/2, Γ 5 µ the axialvector vertex and Γ 5 the pseudoscalar vertex. This equation can be recast in a form that makes explicit the connection between the quark-gluon vertex and the Bethe-Salpeter kernel, Fig. 6. However, it is not an easy task to construct a Bethe-Salpeter kernel given a particular truncation of the quark DSE or arbitrarily specified quark-gluon vertex. Regardless of this difficulty, progress has been made and constructive schemes do exist in the literature as we will discuss below. Fig. 6. The flavour non-singlet axWTI written such that the connection between the quark self-energy and Bethe-Salpeter kernel is explicit. All propagators are dressed, with wiggly and straight lines showing gluons and quarks respectively. The crossed circle indicates the insertion of a γ 5 .
Suppose we choose the simplest possible truncation of the quark-gluon vertex and replace the (reduced) quarkgluon vertex by a minimally dressed tree-level vertex Here, Γ YM represents a dressing function, dependent only upon the square momentum of the gluon. With regards to the quark DSE, this is called the rainbow approximation due to the rainbow-like appearance of diagrams that it effectively resums. The dressing function is chosen to provide a non-perturbative enhancement of the quark-gluon vertex and compensates for lost strength from the missing basis structures. From the axWTI portrayed in Fig. 6 one finds, by inspection, that the corresponding symmetry preserving interaction kernel is This represents a single gluon exchange, forming the socalled ladder part of the rainbow-ladder truncation. Before we consider the constructive scheme for the BS kernel, we make mention of a different approach followed in Ref. [42]. There, rather than constructing a BS kernel for a given truncation of the quark-gluon vertex, they instead constrain the BSE for the meson amplitude itself giving rise to a chiral symmetry preserving truncation scheme. One caveat, however, is that the quark-gluon vertex is forced to satisfy the Abelian Ward-Takahashi identity which fixes the relative strengths of the longitudinal vector and scalar components. The relative signs and strengths of these terms are important in generating the right spin-orbit splitting in different meson channels; any deviations from QCD must be compensated by introducing corrections in different vertex components. Regardless of these limitations, the relative impact of the individual transverse vertex components on meson masses can be successfully tested in this approach and provides for important physical insight.
Since our truncation has a diagrammatic representation, we wish to employ the constructive scheme to obtain our chiral symmetry preserving kernel. It is well-known that one may relate the Bethe-Salpeter kernel K to the quark self-energy Σ by means of the functional derivative [12] a process which is understood to take place in the presence of bilocal external sources. In momentum space, this procedure can be thought of as the cutting of internal quark propagators and their replacement by the Bethe-Salpeter wavefunction, together with the appropriate injection of the mesons total momentum, P.
Since we include the pion back-reaction according to Refs. [34,24,25], we separate the Bethe-Salpeter kernel into non-resonant Yang-Mills corrections (YM) and resonant contributions (pion) The kernel K pion tu;rs (q, p; P), consistent with Eq. (7) and satisfying the axWTI in the chiral limit is K pion tu;rs (q, p; P) = where l ± = (p + q ± P)/2 and D π (p) = 1/ p 2 + M 2 π . The remaining kernel K YM tu;rs (q, p; P) is obtained from Eqs. (3, 6) via the cutting procedure that follows from Eq. (24). This yields the Bethe-Salpeter equation portrayed in Fig. 7, where the pion exchange contributions are indicative of the kernel given above. Note, that to preserve the appropriate charge conjugation properties we have symmetrised the kernel. Now we discuss the different diagrams appearing on the right-hand side of Fig. 7. The first term is that which arises from the rainbow-ladder truncation, and is here associated with the inhomogeneous term present in the DSE of the quark-gluon vertex. The next two loop diagrams stem from the dominant non-Abelian corrections to the vertex that were considered in Ref. [26]. These appear as selfenergy contributions to the quark-gluon vertex and may be included without difficulty in the Bethe-Salpter kernel. The only problem is that these vertices must be calculated for complex incoming and outgoing quark momenta, but is just an exercise in using appropriate numerical methods. Introduction of the sub-leading Abelian vertex correction provides the next three new terms. Two of these are again in the form of self-energy contributions to the vertex and can trivially be added. However, the crossed-ladder kernel deserves further comment being as it is a significantly more complicated non-planar two-loop integral. On application of standard projection methods in order to access the scalar dressing functions of our Bethe-Salpeter amplitude, we obtain a two-loop integration kernel that most compilers will baulk at [43,44]. There are several different paths one may take to reduce this complexity. We chose to construct our Bethe-Salpeter kernel in spinor space using an explicit matrix representation of the Dirac algebra. This has the merit of reducing a million line kernel to a few lines of matrix multiplications. It also permits us to, in essence, 'uncross' the gluon-ladder which simplifies the overlapping integration boundary. The remaining two diagrams of Fig. 7 correspond to the pion back-reaction kernel, which is given explicitly in Eq. (26). This is the full truncation of the Bethe-Salpeter equation that we consider here. Only the first diagram comes from the rainbow-ladder truncation, with everything else classified as 'beyond the rainbow'. Because we do not trivialize the momentum dependence of our propagators, the calculation is genuinely two-loop and we are not restricted as to which meson channels we can investigate. Moreover, such a truncation allows us to consider the dominant non-Abelian corrections [26], together with the sub-leading Abelian corrections [16]. We can also include hadronic unquenching effects according to Refs. [34,24,25]. We are now in a position to perform a systematic comparison of the impact of each vertex correction on the spectrum of light mesons.

Results
Using the truncation outlined in the previous sections, we consider the effect of systematically including different corrections beyond the rainbow-ladder. Because of simplified form of the gluon propagator we employ such a comparison is qualitative; the consistently truncated quantitative study employing Green's functions as calculated from other DSEs will be reported elsewhere. However, because the approximations we make do not trivialise the dynamics and we evaluate the numerics in full, the qualitative statements we make should be valid in general. An further advantage of using the gluon propagator of Eq. (4) is that it allows us to make a direct comparison with other truncations in the literature that have employed the same form of the interaction [27,16,26,42].
We present results in Table 3. The benchmark for comparison, for which the parameters of the interaction are fixed to meson observables, are those of the rainbowladder (RL) approximation. By including the dominant non-Abelian vertex correction (NA) we find only a small change in the pion mass, as expected since it is a pseudo-Goldstone boson. For the remaining bound states we typically see an enhancement of the mass of the order Table 3. Masses for a variety of mesons calculated using rainbow-ladder (RL), additional corrections from the Abelian (AB) and non-Abelian (NA) diagrams, and with the pion back-reaction (PI). Masses are given in MeV.
† For Ref. [16] we quote our own results as obtained from their truncation scheme without use of the real-axis approximation.
of 100-200 MeV, depending upon the meson channel.
That is, inclusion of the non-Abelian diagram is repulsive in the meson channels considered here. In Fig. 8 we show the vector meson mass as a function of the pseudoscalar mass with the non-Abelian correction included. Note that the slope will be changed with respect to rainbow-ladder due to the implicit quark-mass dependence now present in the quark-gluon interaction. We compare these results with the study contained in Refs. [15] (see also [17]). Therein the authors employ a resummation of the Abelian vertex and use the Munczek-Nemirovsky (MN) delta-function to model the gluon. They recognise that this correction is sub-leading with respect to the non-Abelian diagram (considered here) and compensate by altering the colour factor accordingly. However, this gives rise to substantial attraction in the accessible pseudoscalar and vector mesons, as seen in Table 3. Since the interaction and parameters are not similar to ours, we quote their rainbow-ladder results in parenthesis. The precise magnitude of the effect is not important, but we see attraction of the order of 5% in the pseudoscalar and of the order of 100 MeV in the vector. This is in contradistinction to the repulsive corrections that we observe and is seen to be an artefact of the MN interaction or a result of the different kinematical dependences of the vertices.
If we turn off the dominant non-Abelian correction and instead consider the sub-leading Abelian diagram (AB) we obtain the results labelled AB in Table 3. Based on the relative strength of the correction due to its colour factors we expect the results to be N 2 c suppressed, and indeed this is the case. The mass of the pion is protected by chiral symmetry and so it receives negligible contributions from such corrections beyond rainbow-ladder. For the scalar, vector and a 1 axial-vector we see mass reductions of the order of 20 − 40 MeV, while for the b 1 there is a slight repulsion of 3 MeV. This gives a strong indication that such Abelian corrections to rainbow-ladder are generally small and attractive.
The result of our Abelian correction can be directly compared with the study of [16]. In Table 3 we present the results of their truncation using our updated numerical methods; that is, we calculated the quark propagator for complex values of the momenta without recourse to the real-value approximation. For the pseudoscalar, scalar and vector mesons we have qualitatively similar behaviour; they also see attraction in these meson channels. However, in this hybrid truncation with a Munczek-Nemirovsky gluon in combination with Eq. (4) we are not able to find bound-states of the axial-vectors. Instead, our quark propagator develops non-analyticities in the complex plane which restrict our bound-state calculation to masses below 800 MeV. If we consider the original results of Ref. [16] for the axial-vectors, a 300 MeV repulsion was observed in the b 1 channel, compared with our 3 MeV repulsion in the full two-loop calculation. This gives a strong indication that use of the Munczek-Nemirovsky delta-function gluon can give rise to large model artefacts. This gives rise to a degree of uncertainty in any qualitative predictions made. Now, if we consider both the Abelian and non-Abelian corrections to the quark-gluon vertex together we might expect the effect on the meson bound-state to stack. That is, for the vector meson we would expect to see a ∼ 120 repulsion from the NA correction, with a ∼ 25 attraction from the AB correction, resulting in a bound-state mass of ∼ 860. Instead, we see that the Abelian corrections are heavily suppressed by the non-Abelian ones for all meson channels, giving results that are almost identical to the those from the NA diagram alone. The conclusion is that naively adding and subtracting the results of different independent studies without performing the combined dynamical calculation can be misleading. It is thus reasonable, with the interaction model presented here, to ignore the Abelian diagram completely whenever we take the dominant non-Abelian diagram into account.
Finally, having demonstrated that the Abelian diagram is significantly suppressed in the meson spectrum, we consider the dominant non-Abelian corrections beyond rainbow-ladder with unquenching effects in the form of pion exchange (PI). This gives rise to the fifth row of Table 3. As expected, the inclusion of a pion exchange kernel is generally attractive, with an 80 MeV reduction of the vector meson mass with respect to the RL+NA result. This demonstrates, as suspected in [13,26] the near cancellation of beyond-the-rainbow corrections with unquenching effects in this channel. Here, this is not an 19th International IUPAP Conference on Few-Body Problems in Physics exact mechanism but the result of dynamical combinations of attractive and repulsive components of the quark-gluon vertex; a more accurate picture will be revealed when improved approximation schemes are employed that permit quantitative study. Finally, we make a comparison with the truncation of Ref. [42]. Rather than extend the rainbow-ladder truncation by diagrammatically representable corrections to the quark-gluon vertex, they constrain the longitudinal part of the quark-gluon vertex with the aid of the Abelian Wardidentity. This is equivalent to writing the quark-gluon vertex as Γ µ (k, p) = Γ NA q 2 Γ µ BC (k, p) , which is the Ball-Chiu Ansatz as obtained in QED, combined with a dressing function that provides the non-Abelian character of the vertex. Note that such a separable ansatz has been extensively used in studies of the quark-propagator [47,48]. Their results, taken from the figure contained in Ref. [42], are shown in Table 3. This ansatz has the virtue of providing a significant scalar contribution to the quark-gluon interaction that will prove useful in the study of mesons containing significant orbital angular momentum.

Discussion and Outlook
We improved upon previous beyond-the-rainbow investigations, which focused upon pion unquenching corrections [34,24,25] and the dominant non-Abelian corrections to the quark-gluon vertex [26] by including the sub-leading Abelian corrections to the quark-gluon vertex. This enables us to make a direct comparison with existing works in the literature that employ trivial [13,14,15,17] and finite width [16,49] ansätze for the gluon propagator. We investigated the impact of these corrections beyond-rainbow ladder independently and in combination, finding that the Abelian-like corrections are further suppressed relative to the non-Abelian ones in dynamical calculations. This suggests that naively adding and subtracting the relative mass shifts from separate investigates prove to be misleading. We also found that the MN delta-function ansatz gives can give rise to spurious results and model artefacts, thus rendering many conclusions with these models qualitatively unreliable. By including the pion back-reaction in combination with the leading non-Abelian vertex correction we see that there is a degree of cancellation between the two corrections in the vector meson channel.
The study we presented involves many technical developments in the calculation of the quark propagator and quark-gluon vertex at complex momenta. For highly non-trivial two-loop diagrams such as the crossed-ladder kernel, new and improved solution methods were required to render the calculation tractable. Despite these developments, the investigation remains qualitative and serves only as an exploratory study of the relevance of the different contributions beyond rainbow-ladder considered thus far.
Future investigations will focus upon including realistic input from the ghost and gluon sector of the theory [19], appropriate dressings the internal vertices contained within the Abelian and non-Abelian corrections [23] and with input provided from the DSE of the three-gluon vertex [21]. An extension to strange and charm quarks is desirable and possible within our scheme. These are currently under investigation and the results will be reported elsewhere.
In addition to choosing a more theoretically motivated truncation scheme, it would be interesting to introduce further unquenching effects to our Bethe-Salpeter kernel and also our gluon propagator. Such sophisticated studies such as these are left for the future.