Applications of Diagrammatic Renormalization Methods in QCD Sum-Rules

In QCD sum-rule methods, the fundamental field-theoretical quantities are correlation functions of composite operators that serve as hadronic interpolating fields. One of the challenges of loop corrections to QCD correlation functions in conventional approaches is the renormalization-induced mixing of composite operators. This involves a multi-step process of first renormalizing the operators, and then calculating the correlation functions in this mixed basis. This process becomes increasingly complicated as the number of operators mixed under renormalization increases, a situation that is exacerbated as the operator mass dimension increases in important physical systems such as tetraquarks, pentaquarks, and hybrids. Diagrammatic renormalization provides an alternative to the conventional operator renormalization approach. Diagrammatic renormalization methods are outlined and applied to a variety of QCD sum-rule examples of increasing complexity. The results are benchmarked, and the diagrammatic method is contrasted with the conventional operator mixing approach. Advantages and conceptual interpretations of the diagrammatic renormalization approach are outlined and technical subtleties are explored.


Introduction
QCD sum-rules are founded on the concept of quark-hadron duality where QCD composite operators, with appropriate quantum numbers and valence content, are used as interpolating fields to relate hadronic properties to QCD correlation functions of composite operators [1,2] (see also QCD sum-rule reviews in Refs. [3][4][5][6]). From this basic formulation, it is immediately evident that renormalization of the underlying composite operators in the correlation function is a necessary aspect of QCD sum-rule calculations. In most cases, the composite operator renormalization effects enter at next-to-leading loop order (NLO), but in cases of cross-correlators between different currents (e.g., glueballs and quark-antiquark mesons [7,8]), composite operator renormalization enters at leading-order (LO). Thus, in the conventional approach of renormalization constants (or counter-terms), renormalization of QCD correlation functions requires two distinct steps: first, replacement of bare quantities (couplings, masses, and fields) with their renormalized counterparts and second, the renormalization of the composite operators themselves. Both steps of this conventional renormalization approach must be implemented correctly to cancel the divergences within QCD correlation functions used to predict hadronic properties.
where D = 4 + 2ǫ represents the spacetime dimension for dimensional regularization and Γ i represents the collective quantum numbers associated with the currents. In many cases, the correlation function is diagonal with Γ 1 = Γ 2 , but non-diagonal correlation functions with Γ 1 = Γ 2 occur when studying mixed interpretations of hadrons (e.g., mixing of glueballs and quark-antiquark mesons). Generally, some projection operator is applied to (1) yielding a scalar function Π(Q 2 ) where Q 2 = −q 2 . (See (9) for an example.) Based on its high-energy behaviour, Π(Q 2 ) satisfies a dispersion relation, where the Π (n) (0) are subtraction constants and t 0 is a field-theoretical threshold. Except for cases where a low-energy theorem exists (see e.g., Ref. [27]), the subtraction constants are unknown and will typically be divergent. However, the divergent subtraction constants are eliminated in QCD sum-rule methods (e.g., by taking a sufficient number of Q 2 derivatives or applying some integral transform) [1][2][3][4][5], and hence, they do not affect physical predictions. From a field-theoretical perspective, the Π (n) (0) are local divergences because of their Q 2 polynomial structure, a property that becomes relevant in diagrammatic renormalization.
In conventional renormalization, a renormalized composite operator (i.e., current) J R (x) is expressed as a sum of renormalization constants and bare operators and hence, in general, there will be renormalization-induced mixing of n composite operators. Calculation of the Z i first involves developing the operator basis from the general principles [9][10][11]15,16] and then calculating Green functions containing the bare operator J B and the replacement of the bare coupling/masses with their renormalized versions. If the entire renormalization procedure is implemented correctly, then all non-local divergences will cancel, and only local (subtraction term) divergences will remain. Thus, in a basis of n mixed bare operators, one anticipates the need to calculate at least n Green functions for renormalization of the composite operator and then n(n + 1)/2 bare-operator correlators, with potentially multiple Feynman diagrams in each case. Thus, the calculation of QCD correlation functions becomes increasingly demanding in conventional renormalization as n increases (e.g., n ∼ 10 for tetraquark systems). Furthermore, there are very few intermediate benchmarks to ensure accuracy in the final result. These challenges of computational efficiency and accuracy in conventional renormalization are the motivation for developing diagrammatic renormalization methods for QCD correlation functions.
The essential idea of diagrammatic renormalization is to renormalize each Feynman diagram G that occurs in the perturbative expansion of bare quantities through a subtraction process. Following Ref. [15], the renormalized diagram R(G) is obtained by first removing all (non-local) subdivergences to construct R(G), and then applying a counterterm C(G) to remove any remaining local divergences In the case of QCD correlation functions, the local divergences correspond to the subtraction constants which are already eliminated when forming QCD sum-rules, so the process is slightly simpler: for a correlation function diagram G, it is only necessary to constructR(G) by recursively subtracting from the bare diagram U (G) subdivergences occurring in subdiagrams γ where T isolates the divergent part of the subdiagram γ (e.g., ǫ-expansion terms that diverge as ǫ → 0 in minimal-subtraction schemes). Thus, in (5), each subdiagram is replaced with a counterterm diagram that subtracts its divergent part. After every diagram has been renormalized via (5), the resulting sum of diagrams is the renormalized correlation function with all mass/coupling parameters understood as their renormalized versions at the renormalization scale ν associated with the subtraction scheme T (typically MS or MS scheme). The entire process of renormalization-induced operator mixing in the conventional renormalization approach is obviated, leading to considerable increases in computational efficiency because it is only necessary to compute and renormalize the diagrams associated with the bare correlator. Because the process (5) must cancel all non-local divergences for each diagram, diagrammatic renormalization has a built-in internal diagnostic to improve calculational accuracy. These features will be illustrated in the applications presented in Section 3.

Vector Mesonic Correlation Function for Light Quarks
The QCD correlation function of light quark vector (mesonic) currents is the most familiar example of a QCD correlation function Because the current is conserved, the correlation function (6) can be expressed via a single form factor It is well known that the vector current (7) does not require any renormalization factors [42] (see also [16]) so that [i.e., Z = 1 in (3)], and hence, the vector current is the simplest example to illustrate diagrammatic renormalization. The Feynman diagrams for Π v are shown in Fig. 1 up to NLO (two-loops), and, in the light-quark chiral limit, Π v can be expressed in the form where the first term represents the LO Diagram a of Fig. 1, ν is the MS renormalization scale, and the NLO Diagrams b-d (i.e., those proportional to α π ) are parameterized by the coefficients {A, B}. Polynomial terms (i.e., non-logarithmic contributions) are omitted in (11) because they represent dispersion-relation subtractions that do not contribute to the QCD spectral function sum-rules. The Fig. 1 Feynman Diagram d represents a gluon-exchange topology, while Diagrams b and c correspond to the self-energy topology. The NLO results for each topology in conventional renormalization (which in this case is trivial) are shown in Table 1 for the gluon propagator in both Landau and Feynman gauges. Although the conventional renormalization process is trivial, in Feynman gauge, the individual topologies contain problematic L/ǫ nonlocal divergences that only cancel in the total result. However, there is considerable simplification in Landau gauge because the quark self-energy is zero in the chiral limit, and therefore, the exchange topology is free of non-local divergences. Taking into account normalization conventions, the total NLO result for Π v in Table 1 agrees with standard calculations (see e.g., Ref. [16]) and are gauge-independent as required for a gauge-invariant current.

Gauge
Feynman Landau Topology Exchange Self-Energy Total Exchange Self-Energy Total  (11) for the Feynman diagrams of Fig. 1 for Feynman gauge (ξ = 1) and Landau gauge (ξ = 0). The self-energy entry represents the sum of the two diagrams, and hence, the total contribution is the sum of the exchange and self-energy entries.
Following the general process outlined in Section 2, diagrammatic renormalization for the light quark vector correlation function first requires isolation of the subdivergences arising from the one-loop subdiagram topologies shown in Fig. 2. The resulting divergent parts of the subdiagrams, referenced to the topology of the original diagram, are given in Table 2. Counterterm diagrams of Fig. 3 are generated from the subdivergences of Fig. 2, and then subtracted to obtain the renormalized diagram given in Table 3 for Feynman gauge and for Landau gauge. As outlined in Section 2, within the renormalized entries, the coupling is interpreted as α(ν) in MS scheme. Note that the counterterm for the LO diagram (Diagram a in Fig. 3) is purely local and, as discussed above in Section 2, is therefore ignored because it corresponds to a dispersion relation subtraction term that does not enter QCD spectral function sum-rules. As in conventional renormalization, Table 3 exhibits considerable simplifications in Landau gauge because the quark self-energy is zero in the chiral limit. The agreement between the total NLO result for Π v in conventional and diagrammatic renormalization (compare Tables 1 and 3) provides a valuable benchmark for the application of diagrammatic renormalization methods for QCD sum-rules.

Gauge
Feynman Landau Topology Exchange Self-Energy Exchange Self-Energy 0 0  In principle, there is an additional subdiagram shown in Fig. 4 along with its counterterm diagram. This counterterm diagram is zero because of the massless tadpole. However, even if the counterterm diagram was non-zero (e.g., the same topology but with a massive line), the external momentum would not enter the counterterm diagram, and hence, the subtraction would correspond to a dispersion-relation subtraction constant which does not contribute to QCD sum-rules.
Another technical subtlety in the diagrammatic method that requires attention is the role of Ddependent factors that occur in the calculation of a subdivergence. Such D-dependent factors could come from contractions of γ matrices, for example. When computing the divergent part of a subdiagram, it is essential to include all D-dependent factors associated with that subdiagram. Failure to do so will result in errors for the finite parts. For example, with computational tools such as FeynCalc [29][30][31], it may be tempting (for coding simplicity) to defer Dirac/Lorentz algebra until computation of the counterterm diagram, but this leads to errors in the final result. Similarly, it is also important to maintain the ordering of the quark propagators and other Dirac matrices when calculating the subdiagrams (e.g., Fig. 2 c and d).
Despite the simplicity of the vector current renormalization (10), some of the advantages of diagrammatic renormalization are already evident. In particular, no knowledge of the renormalization properties of the vector current were needed in the diagrammatic approach, an aspect that will be illustrated more powerfully in the subsequent examples where the renormalization factors are non-trivial. Furthermore, the non-local divergences for every diagram must cancel against the counterterm diagrams generated from the subdivergences, which provides a self-consistency check at the level of each individual diagram. In conventional renormalization, divergences only cancel in the sum of diagrams, making it more difficult to isolate calculation errors at intermediate stages.  Fig. 1, where the square denotes the subdivergence insertion, ⊗ denotes either the vector or scalar composite operator Feynman rule, and the double line represents the external momentum q. Diagram a is a counterterm for the LO diagram in Fig. 1.  (11) for the Feynman diagrams of Fig. 1 in Feynman gauge ξ = 1, then Landau gauge ξ = 0. The bare entries are repeated from Table 1, and the factor of 2 in the exchange counterterm entries is a result of the two identical counterterm diagrams in Fig. 3. The (shaded) renormalized entries are obtained by subtracting the counterterm from the bare result. The self-energy entries also represent the sum of two diagram, and hence, the total contribution is the sum of the renormalized exchange and renormalized self-energy entries (shaded columns).

Scalar Mesonic Correlation Function for Light Quarks
The QCD correlation function of light-quark scalar mesonic currents extends the vector current analysis of Section 3.1 to a situation where the composite operator renormalization is multiplicative but non-trivial (Z = 1). The scalar correlation function is defined as Conventional renormalization of the scalar current is given by [42] (see also [16]) where the notation Z m is a reminder that the scalar current renormalization constant and quark mass renormalization are related. The Feynman diagrams for Π s are shown in Fig. 1 up to NLO (two-loops), and, in the light-quark chiral limit, Π s can be expressed in the form where the first term represents the LO diagram of Fig. 1 and the NLO corrections are parameterized by the coefficients {A, B}. As before, polynomial terms are omitted in (11) because they represent dispersionrelation subtractions that do not contribute to the QCD spectral function sum-rules. The NLO results for each topology in conventional renormalization are shown in Table 4 for the gluon propagator in both Landau and Feynman gauge. For conventional renormalization, the factor Z 2 m (with Z m from (14)) combines with the LO diagram to generate the NLO renormalization-induced corrections in Table 4.
Conventional renormalization is more subtle for the scalar case with cancellations of L/ǫ non-local divergences via the renormalization-induced contributions in both Landau and Feynman gauge. As in the vector case, there are some simplifications in Landau gauge because the quark self-energy is zero in the chiral limit. Taking into account normalization conventions, the total NLO result for Π s in Table 4 agrees with standard calculations (see, e.g., Refs. [43][44][45]) and are gauge-independent as required for a gauge-invariant current.

Gauge
Feynman Landau Topology Renormalization Exchange Self-Energy Exchange Self-Energy Total  (15) for the Feynman diagrams of Fig. 1 for Feynman gauge (ξ = 1) and Landau gauge (ξ = 0). The self-energy entry represents the sum of the two diagrams, and hence, the total contribution is the sum of the renormalization, exchange, and self-energy entries. The gauge-independent renormalization entry represents the LO diagram combined with the renormalization factor Z 2 m expanded to order α. The total result is gauge-independent as expected.
Proceeding in the same way as the vector case, diagrammatic renormalization for the light quark scalar correlation function begins with isolating the subdivergences arising from the one-loop subdiagram topologies shown in Fig. 2. The resulting divergent parts of the subdiagrams, referenced to the topology of the original diagram, are given in Table 5. Note that the self-energy subdivergence does not depend on the current (compare Tables 2 and 5) because the self-energy subdiagram is isolated from the current vertex.
The Fig. 1 renormalized diagrams obtained by subtracting the Fig. 3 counterterm diagrams generated by the Fig. 2 subdivergences (see Table 5) are given in Table 6 for Feynman gauge and for Landau gauge. As before, the counterterm for the LO diagram (diagram a in Fig. 3) is purely local and, therefore, is ignored because it corresponds to a dispersion-relation subtraction term that does not enter QCD spectral function sum-rules. As in conventional renormalization, Table 6 exhibits considerable simplifications in Landau gauge because the quark self-energy is zero in the chiral limit.

Gauge
Feynman Landau Topology Exchange Self-Energy Exchange Self-Energy  Fig. 2 for the scalar current correlation function in both Feynman and Landau gauge. The subdiagrams are classified by the exchange or self-energy topology of the Fig. 1 diagram from which they originate. The results for the two exchange subdiagrams and the two self-energy subdiagrams are equal, and the self-energy entries are identical to Table 2 because the subdiagram is isolated from the current vertex. The momentum p represents the quark loop momentum flowing through the self-energy.
Gauge Topology Exchange Self-Energy  (11) for the bare Feynman diagrams of Fig. 1 and the counterterm diagrams of Fig. 3 in Feynman gauge ξ = 1, then Landau gauge ξ = 0. The bare entries are repeated from Table 4, and the factor of 2 in the exchange counterterm entries is a result of the two identical counterterm diagrams in Fig. 3. The (shaded) renormalized entries are obtained by subtracting the counterterm from the bare result. The self-energy entries also represent the sum of two diagram, and hence, the total contribution is the sum of the renormalized exchange and renormalized self-energy entries (shaded columns).
The agreement between the total NLO result for Π s in conventional and diagrammatic renormalization (compare Tables 4 and 6) is quite remarkable because the diagrammatic method has not used the Eq. (14) renormalization properties of the scalar composite operator. This provides an important benchmark demonstrating that diagrammatic renormalization methods can calculate QCD correlation functions without any knowledge of the underlying renormalization properties of the composite operator.
An important technical subtlety is embedded in the scalar current correlation function. Because the coupling in the renormalized correlation function is the MS scheme coupling α(ν), the non-zero value for B in Π s implies that the renormalization-group (RG) equation for the QCD spectral function sumrule will have an anomalous dimension term. In the conventional renormalization approach, this is most easily implemented with a quark mass prefactor in the current (13) (see, e.g., [45]) because Z m in (14) is the quark mass renormalization factor. Although the diagrammatic method is seemingly oblivious to the renormalization of the scalar operator, its effect implicitly emerges from the diagrammatic calculation, and the anomalous dimension can be extracted from the correlation function. Standard RG methods for QCD sum-rules [46] can then be applied after extracting the anomalous dimension from the correlation function.

Scalar and Vector Mesonic Correlation Function for Heavy Quarks
The analysis of vector and scalar mesonic correlation functions will now be extended to massive quarks with an emphasis on how the diagrammatic renormalization methods are influenced by the quark mass m. Thus, the focus will be on the divergent parts and the diagrammatic renormalization method, and hence, the lengthy expressions for the finite parts will be omitted.
The inclusion of quark mass does not modify the exchange topology subdivergences of Fig. 2, so the exchange topology results of Tables 2 and 5 remain valid. As noted above, the self-energy subdivergences for massive quarks given in Table 7 do not depend on the current.

Gauge
Feynman Landau im α π 1 ǫ Table 7: Divergent parts of the heavy-quark, one-loop self-energy subdiagrams of Fig. 2 (i.e., subdiagrams originating from a self-energy topology in Fig. 1) in both Feynman and Landau gauge. The results for the two self-energy subdiagrams are equal and are identical for the scalar and vector correlation functions because the self-energy is isolated from the current vertex. The momentum p represents the quark loop momentum flowing through the self-energy.
The divergent parts of NLO contributions from Fig. 1, the counterterm diagrams of Fig. 3 obtained via the Fig. 2 subdivergences (Tables 2, 5 and 7), and the resulting renormalized diagrams have the general form where A and B are polynomials in w and Γ ∈ {s, v}. Because we find that all A contributions are local, they correspond to a dispersion-relation subtraction that can be ignored. Similarly the LO counterterm diagram can be ignored as previously discussed. However, theL divergent structure is problematic, and hence, B = 0 must result for the renormalized diagrams. The NLO divergent parts of the heavy-quark vector and scalar correlation functions for the bare diagrams of Fig. 1 and the counter-term diagrams of Fig. 3 generated by the Fig. 2 subdivergences (see Tables 2, 5, and 7) are given in Table 8 for Feynman gauge and for Landau gauge. Similar to the light-quark analyses, Table 8 shows that there are still some simplifications in Landau gauge for the heavy-quark vector and scalar correlation functions.
As is evident from Table 8, the non-local divergences of each diagram are cancelled by its counterterm diagrams, resulting in B = 0 for the renormalized diagrams as required in the diagrammatic renormalization method. Once again, it is remarkable that the diagrammatic renormalization method does not require any knowledge of the conventional renormalization of the underlying vector and scalar composite operators [see Eqs. (10) and (14)] in the correlation function.

Heavy-Light Diquark Correlation Functions
Heavy-light diquark systems are important within constituent diquark models of closed charmccqq tetraquark systems (see e.g., Ref. [47]). These diquark systems have also been studied in QCD sum-rules at NLO for a variety of quantum numbers [48] using conventional renormalization of the diquark composite operators up to two-loop order [49]. Correlation functions of heavy-light diquarks thus provide an interesting QCD system for exploring diagrammatic renormalization methods.   (17) for the bare Feynman diagrams of Fig. 1 and the counterterm diagrams of Fig. 3 in Feynman gauge ξ = 1, then Landau gauge ξ = 0. The vector current is Γ = v and the scalar is Γ = s. The factor of 2 in the exchange counterterm entries is a result of the two identical counterterm diagrams in Fig. 3, and the self-energy entries also represent the sum of two diagrams. Thus, the renormalized diagram is the difference between the bare and counterterm entries, resulting in B = 0 for the renormalized diagrams.
The heavy-light diquark correlation function is defined as where α, ω are colour indices, Γ indicates the quantum numbers, and gauge-invariant information is extracted using the Schwinger string (see Refs. [48,50,51]) where P is the path ordering operator. The heavy-light diquark currents are where C is the charge conjugation operator, T denotes transpose, Q is a heavy quark, and q is a light quark [50,51]. The operators O Γ = γ 5 , I , γ µ , γ µ γ 5 respectively couple to scalar S J P = 0 + , pseudoscalar P (0 − ), axial vector A (1 + ), and vector V (1 − ) heavy-light diquarks. Analogous to (8), the axial vector and vector projections of the diquark correlation functions are given by Up to NLO, the explicit cancellation of the gauge parameter via the Schwinger string (19) has been demonstrated, and it has been shown that S αω [x , 0] = δ αω in Landau gauge [48,50,51]. Thus, up to NLO, the gauge-invariant information content of the diquark correlation function is extracted in Landau gauge where the Schwinger string becomes the trivial identity operator in colour space. Conventional renormalization of the diquark composite operators up to order α is given by [48,49] Thus, in conventional renormalization up to NLO, there will only be renormalization-induced diagrams in the scalar (S) and pseudoscalar (P ) channels. In Ref. [48], the diquark correlation functions have been calculated up to NLO in conventional renormalization, which provides a detailed benchmark for diagrammatic renormalization. The Feynman diagrams for Π (Γ) are shown in Fig. 5 up to NLO (two-loops), their one-loop subdiagrams are shown in Fig. 6, and the counterterm diagrams generated by the subdiagrams are shown in Fig. 7. The diagram topologies are classified by light-quark self energy, heavy-quark self-energy, and gluon exchange. The non-local divergences associated with Figs. 5 and 7 have the general form where m is the heavy quark mass and B is a polynomial in w. The finite parts involve approximately 15 structures including dilogarithms and trilogarithms (see Ref. [48]) and are not presented for brevity, but are discussed below. The subdivergences arising from the Fig. 6 subdiagrams are given in Table 9. The non-local divergent parts of the resulting Fig. 7 counterterm diagrams, the NLO diagrams of Fig. 5, and the renormalized diagrams are given in Table 10. The LO counterterm diagrams are ignored as discussed previously.
As is evident from Table 10, the non-local divergences of each diagram are cancelled by its counterterm diagram, resulting in B = 0 for the renormalized diagrams as required in the diagrammatic renormalization method. It has also been verified that the diagrammatically-renormalized finite parts, where the coupling and mass are the MS-scheme quantities α(ν) and m(ν), agree with the conventional renormalization results of Ref. [48], providing a detailed validation of diagrammatic renormalization methods in QCD spectral function sum-rules. As in all previous examples, it is emphasized that the diagrammatic renormalization J P q Self-Energy Q Self-Energy Exchange Table 9: Divergent parts of the heavy-light diquark one-loop subdiagrams of Fig. 6 in Landau gauge as required for the Schwinger string simplification. The results for the two self-energy subdiagrams are equal to the Landau gauge results in Tables 2, 5, 7 and are identical for all quantum numbers because the self-energy is isolated from the diquark current vertex. The results for the two exchange diagrams are also equal, and O Γ is the appropriate operator for the J P quantum numbers.
method does not require any knowledge of conventional renormalization (22) for the underlying diquark composite operators in the correlation function. At this stage, another calculational efficiency of diagrammatic methods becomes apparent in the diquark analysis. In conventional renormalization, it is somewhat cumbersome to implement mass renormalization in the LO diagram, but no comparable challenges exist in the diagrammatic approach. Furthermore, another diagrammatic self-consistency check becomes evident for the examples developed. In every case where there are no non-local divergences in the bare diagram, all corresponding subdiagrams are finite (i.e., no subdivergence). This property can be used to identify calculation errors within individual diagrams to improve accuracy in loop computations.

Scalar Quark Meson Glueball Mixed Correlation Function
In the diagrammatic renormalization examples presented above, the composite operators were multiplicatively renormalizable, and the correlation functions were diagonal (i.e., contained a single current). The mixed correlator of the scalar light quark meson and scalar glueball operators (relevant to mixing of scalar quark mesons with glueballs [7,8,52]) is defined by J P Topology Exchange Q Self-Energy q Self-Energy  (23) for the bare Feynman diagrams of Fig. 5 and the counterterm diagrams of Fig. 7 in Landau gauge as required for the Schwinger string simplification. The factor of 2 in the exchange counterterm entries is a result of the two identical counterterm diagrams in Fig. 7. The renormalized diagram is thus the difference between the bare and counterterm entries, resulting in B = 0 for the renormalized diagrams.
where J g is the scalar glueball current and J q (x) is identical to the scalar meson current (13). Conventional renormalization of the scalar glueball operator [14,16] is one of the most familiar cases of operator mixing under renormalization, where β 0 = 11 4 − 1 6 n f is the one-loop β function coefficient and, for simplicity, only a single (light) quark flavour of mass m q has been included (extension to additional light flavours is straightforward). At first order in α and to leading order in the chiral expansion of the light-quark mass m q , the general form of Π gq is and, as before, polynomial terms are omitted in (27) because they do not contribute to the QCD spectral function sum-rules.
In conventional renormalization to first-order in α, Π gq is given by the (two-loop) Feynman diagram of Fig. 8 combined with the renormalization-induced Diagram a of Fig. 1 and the −4m q α π 1 ǫ prefactor of the scalar quark operatorqq in (26). (The β 0 ǫ α π prefactor of the glueball operator in (26) leads to a higher-order α contribution to Π gq .) The (gauge-independent) conventional renormalization results given in Table 11 show that the non-local L ǫ divergence from the bare loop is cancelled by the renormalization-induced diagram. In the diagrammatic renormalization method for Π gq , the one-loop subdiagrams for Fig. 8 are shown in Fig. 9. The subdiagrams are classified as either a gluon-vertex or quark-vertex topology, and their resulting divergent parts are given in Table 12. The results for the bare diagram of Fig. 8, the corresponding counterterm diagrams of Fig. 10, and the renormalized diagram are given in Table 13.
The agreement between the total result for Π gq in conventional and diagrammatic renormalization (compare Tables 11 and 13) is an impressive illustration of diagrammatic renormalization methods in situations where there is composite operator mixing in conventional renormalization. The diagrammatic method required no prior knowledge of the underlying composite operator renormalization in (26). Table 11: Conventional renormalization contributions to Π gq in (27). The bare entry represents the Feynman diagram of Fig. 8. The renormalization entry represents renormalization-induced Diagram a of Fig. 1 (with the scalar quark meson ⊗ operator) combined with the renormalization factor −4m q α π 1 ǫ . The total contribution is the sum of the (gauge-independent) bare and renormalization entries. This benchmark example for the quark/glueball mixed correlation function can also be used to build conceptual understanding of the relationship between conventional and diagrammatic renormalization. It is easily seen that subtraction of the Table 12 subdivergence for the gluon vertex topology is identical to the renormalization-induced diagram generated by the coefficient ofqq in (26), and hence, the diagrammatic counterterm diagram and the conventional renormalization-induced term are operationally identical. The analogy can be taken even further by recognizing that the gluon-vertex subdiagram would be used in conventional renormalization to project out theqq operator's contribution in (26). Furthermore, the absence of a quark-vertex divergence in Table 12 illustrates that, in conventional renormalization, theqq operator does not mix with gluonic operators as illustrated in the multiplicative renormalization result (14) for thē qq operator.

Discussion
The Section 3.5 analysis of the mixed correlation function of scalar glueball and quark meson operators provided a conceptual interpretation of the relationship between diagrammatic and conventional renormalization methods. The conceptual connections between diagrammatic and conventional renormalization are still present in the other benchmark examples presented, but these connections are a bit more difficult to discern. For example, in conventional renormalization for Feynman gauge, the Fig. 2 gluon exchange subdiagram for the light-quark vector correlator is used to find the conventional renormalization of the vector composite operator (10). The Table 2 divergence in the gluon exchange subdiagram corresponds to the quark field renormalization for the external quark lines, which establishes that no additional renormalization factor is needed for the vector composite operator (i.e., Eq. (10) corresponds to Z = 1). The  Table 12: Gauge-independent divergent parts of the subdiagrams of Fig. 9 (i.e., subdiagrams originating from Fig. 8).  Figure 10: Counterterm diagrams generated by the subdiagrams of Fig. 9 and associated with the underlying diagram in Fig. 8, where the square denotes the subdivergence insertion, ⊗ denotes the G 2 composite operator Feynman rule, ⊕ denotes theqq composite operator Feynman rule, and the double line represents the external momentum q.

Diagram Bare Gluon Vertex Counterterm
Quark Vertex Counterterm Renormalized   (27) for the bare Feynman diagram of Fig. 8 and counterterm diagrams of Fig. 10. The bare entries are repeated from Table 11. The renormalized entries are obtained by subtracting the counterterms from the bare result. All contributions are gauge-independent.
1 3ǫ prefactors in Table 2. When the self-energy counterterm diagram of Fig. 3 is calculated, the / p factor conspires to cancel one of the propagators so that the self-energy and gluon-exchange counterterm diagrams have the same loop-integration structure and cancel as seen from Table 3.
A similar (but more subtle) conceptual interpretation applies to the scalar current case in Feynman gauge, except that the Fig. 2 gluon exchange subdiagram divergence corresponds to a combination of the quark field renormalization for the external quark lines and the renormalization factor in (14). However, the Fig. 2 self-energy subdiagram is still related to the quark field renormalizations, and hence, the prefactors of the subdivergences in Table 5 are now different. The self-energy subdivergence / p factor again conspires to cancel one of the propagators so that the Fig. 3 self-energy and gluon-exchange counterterm diagrams have the same structure as the Fig. 1 LO bare diagram. However, because Z = 1 for the scalar operator renormalization, the two counterterms do not cancel, as seen in Table 6, but instead, their residual contribution is from the scalar operator renormalization factor (14). From Table 6, the residual contribution from the counterterms is 4 − 2 ǫ , which when subtracted, is identical to the conventional renormalization contribution (which is also proportional to the Fig. 1 LO bare diagram) in Table 4.
Thus, the diagrammatic renormalization results (e.g., Tables 3 and 6) can be parsed in two ways. The bare and counterterm contributions for each diagram can be combined to cancel the divergence in each diagram leading to the total renormalized result emerging from the sum of the renormalized diagrams. Alternatively, the counterterm contributions can be combined to reconstruct the conventional renormalization-induced contribution (e.g., Table 4).
With the conceptual connections between diagrammatic and conventional renormalization now established, the advantages of the diagrammatic method become more apparent. In both the conventional and diagrammatic approaches, similar subdiagrams need to be computed. In the diagrammatic case, these subdiagrams naturally emerge from the underlying diagrams of the loop expansion, but, in the conventional method, one is calculating Green functions containing the composite operator and fundamental fields with only the general guiding principles of composite operator renormalization. In the conventional approach, it may be necessary to carefully choose and analyze the Green functions to project out and disentangle the mixings in the underlying composite operator's renormalization structure. After the entire conventional operator renormalization process is complete, it is then necessary to go back and compute the new diagrams resulting from composite operator renormalization, including any operator mixing. By contrast, the dia-grammatic method does not need to disentangle the structure of the subdiagrams and it does not matter if the underlying structure comes from multiple operators; all that is necessary is to include the subdiagram's divergence into the associated counterterm diagram. Thus, diagrammatic renormalization methods provide a considerable increase in computational efficiency for QCD sum-rule correlation functions.
A consistency check exists at each stage of the diagrammatic approach because the non-local divergences of an individual diagram must be cancelled by the counterterms generated by the subdiagrams. By contrast, in the conventional approach, the consistency check comes only at the final stage of the calculation where the non-local divergences must cancel in the total result summing all bare diagrams and all renormalizationinduced operator mixing diagrams. In the situation where conventional renormalization may involve mixing of many operators (e.g., dimension-six four-quark operators of Refs. [17]), it may be difficult to isolate calculation errors that could exist within any of the renormalization constants, individual renormalizationinduced diagrams, or individual bare diagrams.
In summary, the motivation for this paper is the need to develop efficient renormalization techniques for NLO QCD sum-rule analyses, what with numerous experimental discoveries of exotic hadrons such as tetraquarks, pentaquarks and hybrids (see e.g., Refs. [19][20][21] for reviews). For these NLO QCD sum-rule analyses, it is necessary to renormalize correlation functions of QCD composite operators of high mass dimension (e.g., 15/2 for pentaquarks).
The conventional renormalization of composite operators becomes increasingly complicated as mass dimension increases, possibly requiring a large basis of operators mixed under renormalization. In conventional renormalization, the calculation of a QCD correlation function is a two-step process of first renormalizing the composite operator, and second, calculating renormalization-induced diagrams arising from the operator mixing under renormalization. Thus, there are strong motivations to develop more efficient composite operator renormalization methodologies for QCD sum-rule applications.
Diagrammatic renormalization methods provide a compelling alternative to conventional renormalization by obviating the need to disentangle operator mixing. QCD sum-rule examples for vector and scalar mesons (Sections 3.1, 3.2, 3.3), diquarks (Section 3.4), and scalar quark/glueball mixing (Section 3.5) have been presented in both conventional and diagrammatic approaches to illustrate the validity and advantages of diagrammatic methods in QCD sum-rules, and to highlight subtleties in the diagrammatic method. Key steps of the diagrammatic renormalization methodology in QCD sum-rules are summarized in Section 5 along with results for extracting divergent parts of subdiagrams. It is hoped that the efficiency and internal self-consistency checks of diagrammatic methods will provide the necessary tools for QCD sum-rule practitioners to tackle the very challenging prospect of NLO QCD sum-rule analyses for exotic hadrons.
Acknowledgments TGS and DH are grateful for research funding from the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix: Methodological Summary and Key Results
In this Appendix, a methodological summary is provided to guide the application of diagrammatic renormalization methods to QCD sum-rules for two-point correlation functions, with a particular view to applications to NLO contributions in multi-quark systems such as tetraquarks and pentaquarks. Key results for the divergent parts of one-loop integrals are also provided to support the outlined methodology.
The diagrammatic renormalization methodology begins by constructing all (NLO) loop diagrams for the desired two-point QCD sum-rule correlation function. Each individual diagram is then renormalized using the following steps. Each step in the process includes an action, possible consistency check(s), and items to note.
After all individual diagrams have been renormalized, the renormalized diagrams are combined to obtain the final result for the renormalized QCD sum-rule correlation function, with coupling and mass parameters interpreted as α(ν), m(ν) for renormalization scale ν in the chosen scheme. The final result for the correlattion function should be gauge independent, providing an additional consistency check supplementing those outlined in the above summary. Anomalous dimensions for the correlation function can be extracted from the final result.
Divergent parts of one-loop Feynman integrals used in our examples and that are anticipated to occur for one-loop subdiagram topologies in future NLO applications to QCD sum-rules (including multiquark systems) are now outlined for our dimensional regularization convention D = 4 + 2ǫ. For two-point subdiagram topologies containing one external momentum scale in loop integrals (e.g., diagrams a and b in Fig. 6), useful results for divergent parts are as follows: where the divergent parts are independent of the propagator masses m 1 , m 2 and Eqs. (31) and (32) emerge from gauge parameter dependence and are therefore written for m 1 = 0 to emphasize this point. Similarly, for three-point subdiagram topologies containing two external momentum scales in loop integrals (e.g., Fig. 9), useful results for divergent parts are as follows: g µν 2p λ p ρ + 2q λ q ρ + p λ q ρ + q λ p ρ + g µν g λρ p 2 + q 2 − p · q + permutations + O(ǫ) where the permutations in (37) indicate those necessary to form a completely symmetric tensor in {µ, ν, λ, ρ} and zero in Eqs. (33) and (34) indicate that the integral is finite.