Heavy-quark-mass expansion of vector and tensor currents and intrinsic charm in nucleon form factors

We compute asymptotic expansions for the vector and the tensor current in the limit of large quark masses in perturbative QCD using the framework of the expansion by subgraphs. We apply the results to estimate the influence of heavy quarks on the nucleon’s electromagnetic and tensor form factors.


Introduction
The intrinsic heavy-quark content of light hadrons is a fundamental property of QCD [1][2][3]. As such, it is expected to have significant influence on the nucleon structure or, more explicitly, on those observables that are related to it. Among the most important of these observables are nucleon form factors; their investigation provides a possibility to JHEP01(2016)092 quantify the impact of heavy quarks on the nucleon structure. Form factors parametrize the expectation values of bilinear fermion operators in a single-nucleon state. Thus, the influence of intrinsic heavy quarks on the nucleon structure can be seen in those bilinear fermion operators that contain heavy quark fields. It is natural to employ a large mass expansion to study such operators. One tool to perform large-mass-expansions which is particularly suited to generate operator-product-expansions (OPE's) is the heavy-quarkmass expansion (HQME), which was used in [4] and [5] to expand the vector and the tensor current of heavy quarks in gluon field operators. It was observed in [6] that the HQME-approach is not complete. The problem arises in the estimation of the influence of muons on the electron's anomalous magnetic moment µ e in QED. If the muon mass is treated as large parameter (compared to the electron mass) the HQME can be used to compute the muonic contribution to µ e by expanding the muon's vector current in photon field operators. The first contribution to the expansion is of order 1/m 4 µ [7]. But this result is in contradiction to the outcome of a direct analysis by Feynman diagrams [8], which shows that the contributions of muons to µ e start at order 1/m 2 µ . Due to the similarity of the calculations, the same problems can be expected for the HQME of the vector current of heavy quarks in QCD. In [6] the contradiction is resolved by the suggestion that the 1/m 4 µ -suppressed photon operator has to be renormalized by 1/m 2 µ -suppressed electron operators. Our goal is to generalize this idea to devise a systematic procedure for finding those additional operators that complete the operator product expansions of the HQME.
To formulate a systematic procedure of generating asymptotic expansions of QCDoperators, we use the language of Feynman diagrams instead of that of path integrals, which was used to derive the HQME results in [4].
In this language, rigorous results for the influence of virtual heavy particles are known for a long time. The most prominent is the decoupling theorem [9]: one can find renormalization schemes, in which, to leading order in the mass of the heavy particles, the effects of the heavy degrees of freedom can be neglected altogether. This result was generalized to a systematic expansion that allows to take into account the effects of heavy particles both on the diagrammatic level and on the level of the Green's functions (see for example ref. [10], ref. [11]). The first treatments used momentum subtraction schemes, which allow for a derivation of asymptotic expansions using standard renormalization-techniques. As dimensional regularization became the preferred method for regularization because of its technical simplicity and its preservation of gauge invariance, it was desirable to derive similar results within the minimal subtraction scheme (MS-scheme). In the MS-scheme decoupling is not straightforward: not all effects of heavy particles are suppressed by powers of their masses. Instead, one has to absorb effects of heavy particles in the physical coupling constants of the theory [12]. Nevertheless, a systematic procedure for the asymptotic expansions within the MS-scheme was formulated (see for example [13]). Finally, the technique of expansion by subgraphs [14] as a particular transparent way to find asymptotic expansions of Feynman diagrams was developed. Although this technique works on the diagrammatic level in the first place, the counterterm technique can be used to generalize the expansion to the level of operators as explained in ref. [13].

JHEP01(2016)092
In this work we compute large-mass expansions for the particularly important examples of the vector current and the tensor current using the expansion-by-subgraphs technique and use our results to analyze the validity of the HQME-approach. Section 2 explains the basic idea of asymptotic expansions. In the following sections the expansion is carried out for the vector and the tensor current: in section 3 the relevant Feynman diagrams for the expansion are listed; section 4 outlines the actual computations. Finally, we use the expansion of the vector current to estimate the influence of charm quarks on the electromagnetic and tensor form factors of the nucleon. The results of the operator-productexpansions are summarized in section 5. Applications for the intrinsic charm in nucleon form factors are discussed in sections 6.1 and 6.2.

Asymptotic expansions in heavy quark masses
The theoretical foundations of perturbative asymptotic expansions in renormalizable quantum field theories are discussed in a form that is particularly suitable for our work in refs. [15] and [13]. We shall summarize their results to highlight the main ideas that underlie asymptotic expansions of composite operators in the special case of large-mass expansions.
Consider a renormalizable quantum field theory with one type of particle with large mass M , and several light particles of masses m i ≪ M . The first step to understand the asymptotic behavior of the theory on the perturbative level is the investigation of single Feynman diagrams. We denote a Feynman graph by Γ, its corresponding analytic expression by I Γ ({p}, m, M ) and the superficial degree of divergence of I Γ ({p}, m, M ) by ω Γ .
For an arbitrary Feynman graph Γ we are interested in the kinematic region in which the external momenta, collectively denoted by {p}, are negligible compared to the mass of the heavy particle, that is p ≪ M .
As already mentioned in the introduction, asymptotic expansions can be derived with methods similar to those used for renormalization. To illustrate this remark, we consider a divergent diagram Γ in a theory regularized by a cut-off Λ in a zero-momentum-subtraction scheme. To renormalize the diagram, we state a power-counting theorem: there is a special class of subdiagrams of Γ, called one-particle-irreducible (1PI) graphs with the following property: if all 1PI-subgraphs of Γ have negative mass-dimension, each 1PI-subgraph and Γ are finite. If we succeed in isolating the divergences from the 1PI-subgraphs, we have also isolated the divergences of the complete diagram. To isolate the divergent part of an 1PI-subdiagram without subdivergences, we can perform a Taylor-expansion in its external momenta to reduce the mass-dimension of the diagram until it is convergent. If we perform Taylor expansions on all 1PI-subgraphs in the order as described in the R-operation of the zero-momentum-subtraction scheme, the Taylor expansion is always performed on 1PIsubgraphs without subdivergences; thus, the divergences are contained in the leading terms of the expansion while the remainder of the expansion contains 1PI-subgraphs with negative mass-dimension: it is finite. Expressed differently, the Taylor expansions generate the leading terms in the Λ-expansion of the diagram, as the power counting theorem ensures that the remainder is of higher order in 1/Λ. For large-mass expansions, the large mass M JHEP01(2016)092 plays similar role as Λ does in renormalization. In same way as the divergent terms -the leading terms in the 1/Λ-expansion -are isolated in renormalization, the leading terms in 1/M are isolated in the large mass expansion.
Guided by analogy, we first look for a power counting theorem that distinguishes the class of subdiagrams that rule the 1/M -behavior of the complete graph. Such a theorem can be found in ref. [15].
Let Γ be a diagram with l loops. The asymptotically-irreducible (AI) subdiagrams of Γ are those connected subdiagrams that contain heavy-particlepropagators and cannot be made disconnected by cutting a single light line. The the asymptotic behavior of I Γ ({p}, m, M ) is bounded by where S is that set of disconnected AI graphs containing all heavy lines which has the highest degree of divergence.
Consequently, by introducing an operator that reduces the degree of divergence of the AIdiagrams sufficiently (as the R-operation in momentum subtraction schemes reduces the degree of divergence to render 1PI graphs convergent) one can generate a remainder of an arbitrarily chosen order in the heavy mass. The difference between this remainder and the original (renormalized) diagram will then constitute the correct asymptotic expansion for the diagram. As this procedure is extensively described in [13], we only give the main result. To state it, some notations are required. Let S AI (Γ) be the set of mutually disjoint AI-subdiagrams, S an arbitrary set of disjoint subgraphs. Γ/S denotes the diagram Γ with all elements of S shrunk to a point. I Γ/S • γ∈S V γ is the expression corresponding to the diagram in which the vertex V γ replaces the subgraph γ for each γ ∈ S. R un I Γ/S denotes the R-operation that only acts on those parts of the diagram that do not contain vertices to which the graphs of S were shrunk.
Finally, M aγ γ performs a Taylor expansion of I γ ({p}, m, M ) in the momenta {p} external to the graph γ and the light masses m up to the order a γ .
In [13] it was shown that for a γ = a + ω γ one has, up to O M −a−1 -terms, It was also proven that this result is free of artificial ultraviolet-and infrared-divergences, that is, the divergences of this expression are the same as those for the original diagram. To summarize: in performing an asymptotic expansion in large masses on has to find all sets of mutually disjoint AI-subdiagrams of a diagram Γ. In each subgraph γ of such a set, one has Taylor-expand in external momenta and light masses. The result is a new vertex factor V γ that is reinserted in the original diagram to replace the subdiagram γ. Then the unaffected parts of Γ are renormalized as usual.

JHEP01(2016)092
We now consider the transition from the diagrammatic to the operator level; we want to derive an operator product expansion of an operator O with heavy quarks in terms of light degrees of freedom. We will merely summarize the argument of section 4.4 of ref. [14]. O can be considered as identical to a sum of operators from light degrees of freedom if all their Green's functions coincide. We concentrate on a Green's function with n light external fields Φ and an insertion of O, schematically written as: Peturbatively, this Green's function can be computed in terms of Feynman graphs. For each of the graphs Γ, we can use the relation (2.2). As O contains heavy quarks, it is always part of an AI-subdiagram. We single out this special subdiagram and leave the remaining AI-subdiagrams unexpanded; to complete the elimination of heavy particles in the theory they can be considered after the expansion of O, for example by using an effective lagrangian as obtained in ref. [4]. This procedure yields: We sum over all diagrams Γ contributing to eq. (2.3) and regard Γ/γ as new graph Γ ′ Vγ with a special vertex V γ where γ was shrunk to a point. The sum over Γ can now be divided into a sum over all vertices V that can be generated from AI-subgraphs and a sum over all possible graphs Γ ′ V with vertex V . We can exchange the sum over Γ ′ V and γ. The sum over γ is now over all AI-diagrams γ V with number and type of external lines appropriate for the vertex V : {q}, m, M ) has to be inserted at the location of the special vertex V in Γ ′ V and can be regarded as effective vertex V γ generated by composite operators. We realize that eq. (2.5) is the sum of all diagrams that contribute to eq. (2.3) with an insertion of effective vertices generated by all possible AI-subdiagrams involving O: The operators related to V γ constitute the desired asymptotic expansion. Thus, we have the following recipe: find all AI-diagrams involving O, compute their asymptotic expansion, express the result in terms of effective vertices generated by composite operators and read off the operator product expansion. We now consider the expansion of a heavy operator O in QCD. To expand the operator O of dimension dim(O) to the order a, one has to expand each AI-subgraph γ to the order a γ = ω γ + a on the diagrammatic level. Since in theories with dimensionless coupling constant the dimension of γ is number of external boson-and fermion-fields in γ, it is enough to consider diagrams that satisfy 0 ≤ a + dim(O) − E B − 3 2 E F . For example, at the order 1/M 2 for dim(O) = 3, only diagrams with E B ≤ 5 and E F ≤ 2 have to be taken into account as external fermions have to come in pairs. As we are working in a gauge theory gauge theories, Ward-identities will reduce the actual degree of divergence for diagrams with external vector bosons. This reduces the number of diagrams that have to be taken into account.
In this paper the vector current and the tensor current are considered. An operator product expansion is performed to the order 1/M 2 for the vector current and to the order 1/M for the tensor current. We consider terms up to the order α 3 in the strong coupling constant as this is the lowest order with purely fermionic operators.

The asymptotically-irreducible diagrams for fermion-bilinears
We now use the procedure outlined in section 2 to derive an operator product expansion of the vector and tensor current of a quark of mass M in perturbative QCD for the case of M being large compared to all external momenta, all other mass terms and -to ensure that perturbation theory is applicable -to Λ QCD .
The first step towards an expansion is the identification of the relevant AI-diagrams. First, some the general statements about their structure are discussed. These statements are then used to obtain a list of the AI-diagrams relevant for the expansion of the heavy vector and tensor current to the order α 3 in the strong coupling constant and 1/M 2 or 1/M in the large mass respectively.
The structure of the expansion is constrained by the following argument (see the appendix of ref. [10]). Since the heavy particles should not occur as external particles, every AI-diagram will contain heavy fermions only inside closed loops. Consider the fermion propagators S F (p, M ) of such a loop. If we assume that a regulator that admits the usual definition of γ 5 is used, they satisfy the relations Now, the emerging factors of γ 5 can be absorbed in the vertex factors via This procedure cancels all factors of γ 5 in the heavy-fermion-loop. The diagrams with vector current insertion remain unchanged while the signs of the diagrams with tensor current insertion are inverted. Thus, the expansion of the vector current is in even powers of M , that of the tensor current is in odd powers of M . We can further restrict the types of relevant diagrams. Only gluons can couple directly to the heavy quark loop. This means that each AI-diagram will have a subdiagram that can be considered as off-shell matrix element of the Green's function (3.5)

JHEP01(2016)092
Matrix elements with one external gluon are zero due to color conservation. Let us consider the case of two external gluons. Under charge conjugation, we have Global color conservation requires the matrix element with two gluons to be proportional to the only invariant color-tensor with two indices, δ a b . If one determines the proportionality constant by contracting with another factor of δ a b , one obtains: Applying charge conjugation to both sides of this relation yields: Therefore, for the vector and tensor current one needs to consider only diagrams where the fermion loop is connected to the rest of the diagram by at least three gluon propagators. Similar arguments show that there is only one possible color structure for the matrix elements with three external gluons: due to global color conservation two color structures for this matrix element are possible (8 ⊗ 8 ⊗ 8 contains two different singlets). These can be chosen as the two independent trace structures. In combination with charge conjugation invariance this leads to: For the case of the vector current, QED current conservation requires an additional factor of the momentum entering the diagram at the operator insertion. Dimensional analysis shows that this reduces the power of the heavy mass by one.

JHEP01(2016)092
We will show in section 4.1.1 that if three gluons couple directly to the fermion loop each of the external gluons causes an additional momentum-factor. Taking into account both kinds of momentum factors, the power of the heavy mass is g is the number of external gluons coupled directly to the heavy loop, δ = 1 for the vector current and δ = 0 for the tensor current if no more than three gluons couple directly to the heavy loop. If there are more than three gluons coupled to the loop, we have The general arguments about the large-mass dependence of matrix elements can now be used to give a complete list of relevant AI-diagrams for our expansions as shown in figure 1. Note that for all diagrams given in figure 1 one also has to compute the crossed diagrams. We start with one loop diagrams. Following the arguments given above, the diagrams with three external gluons, shown in figure 1a, give the first contributions for the vector current and the tensor current. For the vector current these diagrams are of order ω γ = −4, for the tensor current of order ω γ = −3. Actually these diagrams can be excluded from the computations at the order we consider. They are computed nevertheless, because we want to compare our results to the HQME-results of ref. [5]. The other diagrams relevant for the expansion have four external legs, as shown in figure 1b. These diagrams have ω γ = −2 for the vector current and of order ω γ = −1 for the tensor current. Our explicit computation and the results of section 4.1.1, however, showed that there are no 1/M 2 or 1/M contributions. This is consistent with the HQME results, which predict that all one-loop contributions are of order 1/M 4 for the vector current and 1/M 3 for the tensor current. Now, consider two-loop diagrams. Since at least three gluon lines must couple to the fermion loop, the relevant AI-diagrams are given by the figures 1c, 1d, 1e and 1f at order g 5 s and by 1g, 1h, 1i and 1j at order g 6 s . The three-loop diagrams are the first which contain external fermions. In the computation we will consider one arbitrary light quark flavor and afterwards take the sum over all light flavors to obtain the complete operator product expansion. If only two of the gluons are coupled to the fermion line, we have at least ω γ = −3 for the vector current and ω γ = −2 for the tensor current. Thus, the only possibility is to couple all three gluon legs to the fermion line, see figure 1k.
We now explain the absence of diagrams with ghost fields. Up to the order α 3 s in the strong coupling constant only external ghost lines can occur. Due to ghost number conservation the number of ghost fields must be even. Charge conjugation invariance shows that diagrams with two external ghost fields vanish: All nonzero diagrams with more than two external ghost fields are of higher order We now turn to those terms in the expansion that contain only gluonic operators. The requirement of local gauge invariance severly restricts their structure: they can contain only field strength tensors and covariant derivatives of them. Since the field strength tensor behaves under charge conjugation as the vector potential, arguments identical to those leading to eq. (3.12) restrict the result for both operators to the schematic structure tr c F 1 , F 2 F 3 . The vector current, which has to be the derivative of some antisymmetric tensor due to current conservation, has the structure ∂ tr c F 1 , F 2 F 3 . Power counting now shows that the gluonic contributions to the expansion should be of order 1/M 4 for the vector current and 1/M 3 for the tensor current to all orders of perturbation theory. Our computations up to O(α 3 s ) confirm these results.

JHEP01(2016)092 4 Computation of the asymptotically-irreducible diagrams
We can now use the relation (2.2) to compute the large mass expansion for the diagrams shown in figure 1. There is one additional simplification. In general, single diagrams in the heavy-quark mass expansion of the vector and tensor current have divergences, but the sum over all diagrams at a given order in the heavy mass and the coupling constant turns out to be finite. Actually even the sums of all diagrams of a given type (as indicated in figure 1) are finite. Thus, it is not necessary to perform the R-operation for them. To obtain the desired expansion, each AI-diagram has to be expanded in small masses and external momenta. For this purpose it is convenient to rescale the integration momenta according to l = Ml (this also makes the integration variables dimensionless). Since the complete Feynman integral is a homogeneous function of dimension ω γ , we have for n-loop integrals I n (M, {m}, {p}): In eq. (4.1) {l i } labels the set of integration momenta, {p} that of external momenta and {m} the set of light masses. Eq. (4.1) clearly shows that an expansion in small masses and momenta is equivalent to an 1/M -expansion. In the discussion of the different Feynman diagrams in the following subsections, we proceed as follows. The algebraic expression for the given diagram is given before rescaling. Prior to the Taylor expansion the momenta are rescaled as described above; the following expressions are given in terms of the rescaled variables.
After the Taylor expansion of the AI-diagrams we obtain several thousand of different tensor integrals and after tensor reduction there are thousands of scalar integrals. In the following we explain how the tensor reduction was done and how each type of diagram was computed. The computation was performed by a Mathematica program which made use of the package 'FeynCalc' [16]. In most cases the expansion of the AI-diagrams given above does not result in scalar integrals. This is either because of uncontracted Lorentz-indices in the diagram or because of scalar products of external and loop momenta occurring in the numerator in the 1/M -expansion of the integrand. We will explain how tensor integrals were reduced to scalar ones. Since the expansion is in all external momenta, the integrals that we obtain effectively correspond to diagrams without external momenta; they are vacuum diagrams. For these, the tensor decomposition is easy: the metric tensor is the only tensor that can occur in the result. All integrals with an odd number of indices vanish.

JHEP01(2016)092
An arbitrary integral F µ 1 ...µ 2n in the expansion can be parametrized as F µ 1 ...µ 2n = A 1 g µ 1 µ 2 · · · · · g µ 2n−1 µ 2n + all other contractions . (4.2) Here A i are Lorentz scalars. In general, there are (2n)!/(n! 2 n ) = (2n − 1)!! different contractions. In many cases, the integral is symmetric under the interchange of some of its indices. Then one can symmetrize both sides of eq. (4.2) such that the number of independent constants A i is further reduced. Now, contracting (4.2) with each of the independent tensor structures in turn, one obtains a system of equations for the A i with scalar integrals on the left-hand site. Solving for the constants and reinserting the solution on the right-hand side of (4.2), one obtains a decomposition of the integral in terms of scalar ones. The computation of the scalar integrals is described for each type of AI integral separately. Note, however, that the Feynman rules for the composite operator O(x) contain an additional factor of e iq·x in comparison to the usual Feynman rules for an identical term in the Lagrangian (this results from the fact that there is no integration over x).
The translation of the results of the Feynman diagrams to the operators occurring in the OPE can be done with different sophistication (for example with the help of the counterterm technique [17] or the method of projectors [18]). Here we simply give operators that yield vertices with the correct external fields and with vertex factors that are those polynomials resulting from the computation of the AI-graphs.
The vertex factor corresponding to an operator O(x) is: It is useful to note that the operator G b αβ corresponds to the projection operator on the transversal component of the momentum that enters the vertex: The generalization to higher powers of gluon fields is straightforward. Additional momentum factors are of course generated by additional derivatives of the gluon field operators. The situation is especially easy for additional factors of the overall momentum flowing into the vertex: . The O(g s )-contribution is of higher order in comparison to the diagrams considered here and can be neglected. We shall denote the OPE as Here '≃' indicates that the relation gives an asymptotic series.  where J α a (x) = Q γ α T a Q (x) denotes the color current. In QCD the color current is not exactly conserved, but it satisfies the relation ∂ α J α a (x) = 0 + O (g s ). The fermion loop that we consider contains only the contributions of LO in g s . Thus, to analyze the loop we can assume that the current is conserved. Taking the divergence of the matrix element, LO terms are only generated by moving the derivative inside the time ordering: The commutators are determined by the transformation properties of the operators under SU(3) c -transformations:

JHEP01(2016)092
This gives the result (4.12) In momentum space, we have (4.13) Similar relations hold for all external momenta p 2 . . . p n . These results allow to draw conclusions for the case of three and four external gluons.
For three external gluons, we have (4.14) The last equality follows from charge-conjugation invariance as discussed above. From this relation (and from the corresponding ones for p 2 , p 3 ) it is clear that (4.15) Now, the power counting arguments given above apply: the matrix element for the tensor current is suppressed as M −3 and that for vector current as M −4 .
In the case of four external gluons: where n = 3 for the tensor current and n = 4 for the vector current. Since the same argument applies for all external momenta, those parts of the matrix elements not suppressed as described above must have at least four projection operators, that is, their Lorentz structure must be such that the contraction of the diagram with any external momentum is zero. Power counting and the general structure of the expansion for tensor and vector current yield a suppression of M −5 and M −6 for contributions with these projectors respectively.

Feynman integrals for the one-loop contribution
The only one-loop integrals to be computed are shown in figure 3. In general, we use the following conventions: momenta of external gluons are incoming, fermion momenta flow JHEP01(2016)092 For the diagrams of the type shown in 3a the analytic expression reads: where we use Γ = γ µ for the vector current and Γ = σ µν for the tensor current. For the diagram 1b with four external lines we obtain where we use l 3 = −q − l 1 − l 2 to simplify the notation (this is the momentum entering the last gluon vertex). The master integrals that occur after the expansion are of the form:

Result for the vector current
First consider the vector current in figure 3a. Up to the order 1/M 4 , the result of expanding the diagram and summing over all permutations of external momenta is: In this form the result clearly shows the expected projector structure: contracting with q µ , (l 1 ) α , (l 2 ) β or (l 3 ) τ we obtain zero. Due to the projector-structure of the results they can be easily rewritten in the form of vertices generated by the insertion of local composite-operators. If we omit the argument x and neglect higher terms in the coupling constant, a short calculation yields: With the application to instanton solutions in mind, we introduce the field strength tensor F µν = g s G µν . Thus, the contribution can be generated by This is exactly the result from the calculation of ref. [5] after continuing to Minkowskispacetime and rescaling the gluon fields. Explicit computation showed that the result for the diagrams of the type shown figure 3b corresponds exactly to the O g 4 s -contribution to eq. (4.22).

Result for the tensor current
For the tensor current a calculation similar to that of subsection 4.1.3 yields the following result: As expected the result vanishes for contractions with (l 1 ) α , (l 2 ) β or (l 3 ) τ . We can again write this in terms of an operator insertion: Thus, the one-loop contribution to the operator product expansion is which coincides with the result of ref. [4] if one takes into account the analytic continuation to Minkowski-spacetime and the rescaling of the gluon fields. As in the case of the vector current, the diagram figure 3b corresponds to the O g 4 s -contribution to equation (4.25).

Two-loop contribution
There are several different kinds of integrals at two-loop level. Their explicit forms are given in the appendices A and B. All relevant integrals can be computed in closed form and after rescaling have the form (4.28)

Three-loop contribution
At the three-loop level, the relevant diagram is shown in figure 4a. The corresponding integrals have the form where m is the mass of any of the light quark flavors. The integrals appearing after expansion have the form (4.30)

JHEP01(2016)092
These integrals are hard to compute directly, but the integration-by-parts (IBP) technique can be used to solve this problem. We discuss this in the next section.

The integration-by-parts technique
The integration-by-parts technique is a well-known tool to reduce large classes of Feynman integrals that are hard or impossible to compute analytically to a small number of so-called master-integrals. The procedure for the case of eq. (4.30) was first formulated in [19]. We will described a slightly modified version of this procedure below. There are several software packages that compute such integrals ('MATAD' [20], 'FIRE' [21]). However, as these integrals occur at an intermediate step of our computation and as one has to compute about 50.000 integrals, it was more convenient to write a program specialized on exactly this type of integral.
We will now show that we can use the IBP-technique to reduce an integral with arbitrary values of the indices α i to integrals which fall in one of the following classes: 1. α 1 ≤ 0, α 2 ≤ 0 or α 3 ≤ 0 2. α 5 ≤ 0 and α 6 ≤ 0.
Afterwards, we show how to compute the integrals for these special cases.
Using the relations one can use the fact that the integral of a total derivative is zero (if the boundary terms vanish), to obtain Each application of eq. (4.37) reduces the sum α 2 + α 4 + α 5 by one. Eq. (4.37) can be used recursively until one of the indices is zero. In the cases α 4 < 0 or α 5 < 0 this relation cannot be used. However, because of the conventional ordering of the indices, this case would imply a negative value of α 6 and it would require α 5 to be the second non-positive index. If finally either α 4 or α 5 is zero, one can stop, if α 6 is already non-positive. If this is not the case, one can interchange the indices to make α 6 the index which is zero. Then one can use the relation: This relation can be used until α 2 or α 5 is 0. We are now in position to discuss special cases. First we consider the case α 1 ≤ 0, α 2 ≤ 0 or α 3 ≤ 0. The cases α 1 ≤ 0, α 2 ≤ 0 can be reduced to the case α 3 ≤ 0. The relevant integral is: If one of the indices α 5 or α 6 is nonnegative, the integral is zero. Expanding the numerator and eliminating the factors of k +l 1 2 ,l 1 ·l 2 andl 2 2 from the numerator by cancellation with terms in the denominator, one is left with integrals of the type:  Further cancellations yield integrals of the form (4.42) If α 1 or α 2 is non-positive, the integral is zero. If all indices are positive, the result is given by (4.27). Even if α 3 is negative, one can use this result as analytic continuation to α 3 < 0. Let us consider the case α 5 ≤ 0 and α 6 ≤ 0.
Renaming and exchanging integration variables in eq. (4.30) yield: If any of the massive lines has non-positive indices, the integral is 0. If α 4 is non-positive, the integral reduces to a product of three integrals like which can be computed using eq. (C.1). If all indices are positive we get: Performing thel 2 -integration and changing integration variables such that p < q, we get

JHEP01(2016)092
Note that in this case all terms in the numerator can be canceled with those in the denominator since all exponents are integers. The result is: At this point we can use eq. (4.27).

Result for the vector current
We now use the techniques of section 4.3.1 to derive the contributions for the vector current.
With α s = g 2 s /(4π) we obtain: where ζ(3) = ∞ n=0 1/n 3 is the zeta-function. If we denote the creation and annihilation operators for the light quark in the three-loop diagram by Ψ and Ψ respectively, it is possible to read off the corresponding operators immediately: (4.50) Note that we can freely interchange derivatives with covariant ones since the difference is of higher order in g s . Using the QCD equations of motion, / ∇ Ψ = −im Ψ and Ψ ← − / ∇ = im Ψ and the fact that the vector current for the light quarks is conserved, ∂ µ Ψγ µ Ψ = 0, we have (4.51) Applying these relations yields: (4.52)

Result for the tensor current
For the tensor the application of the method of section 4.3.1 yields: , (4.53)

JHEP01(2016)092
which, in terms of operators, can be written as: After applying equations of motion and standard relations for the γ-matrices we obtain: (4.55) The result for the tensor current can now be written as: (4.56)

Summary
We can now collect the contributions from the one-, two-and three-loop diagrams to the large-mass expansion of the vector and tensor current of heavy quarks. One has to remember that each light flavor will contribute to the three-loop diagrams. The result is: Here f denotes the sum over all flavors of quarks that are light compared to M . These operator identities can be used to analyze the intrinsic heavy quark content of light hadrons. In the next section we study the intrinsic charm influence on the nucleon's electromagnetic and tensor form factors.

Applications
The expansions of the heavy vector and tensor current in eqs. a single-nucleon state. The sum is over all quark flavors; we concentrate on the contribution of charm quarks, considered as heavy quarks, and regard up, down-and strange-quarks as light. Using eqs. (5.1) and (5.2), those terms involving charm-quarks can be expressed in terms of operators with light degrees of freedom. The expectation values of these operators can be expressed in terms of form factors of the light quarks and gluons. Thus, the influence of charm-quarks on the nucleon structure can be parametrized by the 'light form factors', which is a good starting point to estimate the relative importance of charm-quarks compared to light quarks for the nucleon structure. In section 6.1, this procedure is applied to the electromagnetic current and an estimation of the influence of charm-quarks on the charge radius and anomalous magnetic moment of the nucleon is obtained. Section 6.1 deals with the tensor current and provides estimates for the influence of charm on the tensor form factors of the nucleon.

Heavy-quark contribution to magnetic moment and electromagnetic radii of the nucleon
We can use the results of the previous section to derive a relation for the gluonic contribution to the nucleon's charge radius that is due to heavy quarks. We denote the momentum transfer by q = p ′ − p and the sum of the momenta by P = p ′ + p; the spinors for a nucleon of momentum p and spin σ are denoted by u (p, σ) and u (p, σ) and the nucleon mass by m N . First of all, we need the form factor decomposition for the vector and the tensor current. Using the form factors F 1 q 2 + iσ µν q ν 2m N and F 2 q 2 , the decomposition for the vector current is standard: The tensor current can be parametrized by form factors H T q 2 ,H T q 2 and E T q 2 as [22]: Actually, a parametrization of q ν p ′ , σ ′ | Ψ iσ µν Ψ(0) |p, σ is required: where we introduced the tensor magnetic form factor κ T q 2 = E T q 2 + 2H T q 2 and used the Gordon-identity. Additionally we introduce form factors of the gluon operator

JHEP01(2016)092
Since this operator is an antisymmetric tensor of rank 2, it can be decomposed as the tensor current. Thus, introducing the form factors R q 2 , S q 2 andS q 2 , and Now, we consider the nucleon matrix element of the vector current made of heavy quarks. Using the derived operator identity (5.1) we obtain the heavy quark contribution to the electric form factor G Q E q 2 = F Q 1 q 2 − q 2 /4m 2 N F Q 2 q 2 and the magnetic form We see that if we neglect the chiral corrections of order ∼ m f /m N , the leading ∼ (α s /π) 3 1/M 2 contributions can completely be expressed in terms of singlet electromagnetic form factors. The sub-leading corrections of order ∼ 1/M 4 are expressed in terms of form factors of the gluon operator eq. (6.5). This gives the following result for the heavy quark JHEP01(2016)092 contribution to the electric and magnetic radii r 2 : and (6.10) Finally, the heavy quark contribution to the magnetic moment (in nucleon magnetons) is: We see that the leading ∼ 1/M 2 contribution to µ Q has an additional suppression factor ∼ m f /m N indicating that the subleading correction ∼ 1/M 4 can be numerically more important than the leading one. The result for the heavy quark contribution to the magnetic moments (6.11) coincides with the result obtained in ref. [6]; the results for the electromagnetic radii (6.9), (6.10) are new. From eqs. (6.9), (6.10), neglecting the corrections of order ∼ m f /m N , one obtains the following model independent limits for the electromagnetic radii of the nucleon: We see that the limiting values of the heavy quark contribution to nucleon electromagnetic radii are always negative, so the intrinsic heavy quarks shrink the electromagnetic size of the nucleon. These limiting expression can be useful for lattice QCD simulations. Now we give numerical results for the heavy quark contributions to magnetic moment and electromagnetic radii in the case of the charm quark. We use the quark masses  [23] and the coupling constant α s m 2 c = 0.3. The tensor form factors and the tensor anomalous magnetic moments were computed in the chiral quark soliton model [24,25]. The values for the proton and neutron are H u T,p (0) = H d T,n (0) = 0.92, H d T,p (0) = H u T,n (0) = −0.27, κ u T,p (0) = κ d T,n (0) = 3.03 and κ d T,p (0) = κ u T,n (0) = 1.496 [26,27]. For these calculations isospin symmetry was assumed and the values were evolved from a renormalization scale of µ 2 = 0.36 GeV 2 to µ 2 = m 2 c by the method described in these publications. The numerical result for the charm quark contribution to the electromagnetic radii of the nucleon is: The result for protons and neutrons is the same since the isospin violating effects due to mass difference of u and d quarks are negligible. We see that if the form factors of the gluon operator in eq. (6.5) are numerically of order of one, then the subleading term ∼ 1/M 4 is larger than the leading term for the case of the charm quark. This stresses the importance of the estimates of the gluon form factors in eq. (6.5). One of possibilities to estimate these form factors is the theory of instanton vacuum [28,29]. However, one can easily see that the gluon operator (6.5) is exactly zero on the instanton field. This implies that one has to consider the contribution of the instanton-anti-instanton pair to this operator; hence, the gluon form factors (6.5) are suppressed by the instanton packing fraction in the vacuum. Given such suppression, we can very roughly estimate the size of the gluon form factors as: S ∼ 10 −2 . (6.16) In eq. (6.5) the factor (N c − 2) takes into account that the gluon operator in eq. (6.5) is zero for the case of N c = 2, the factor πρ 2 /R 2 ∼ 1/3 is the instanton packing fraction, and f S ∼ 0.1 is the twist-4 contribution to the nucleon structure functions (see details in refs. [4,30,31]).
Of course, eq. (6.16) provides only a very rough estimation. Nevertheless, it indicates that suppression of the gluon form factors by the instanton packing fraction can be not enough to make the leading ∼ 1/M 2 term in charm quark mass expansion dominant numerically. Therefore, it is important to perform more accurate calculations of the gluon form factors; we shall give detailed a estimate of them elsewhere. Let us note obvious general properties of the gluon form factors (6.5): 1. they are exactly zero for the case N c = 2, 2. they are zero for general self-dual gluon field configurations.
For the heavy quark contribution to the magnetic moment the leading ∼ 1/M 2 term is suppressed by the factor ∼ m f /m N ; therefore we expect that for the case of the charm quark the sub-leading contribution ∼ 1/M 4 is dominant numerically. Indeed, if one takes

Conclusion and outlook
We derived the heavy quark mass expansion of the vector and tensor currents up to three loops using method of 'expansion by subgraphs'. The corresponding operator identities are summarized by eqs. (5.1), (5.2). We applied these results to analyze the influence of the intrinsic charm on the nucleon's electromagnetic and tensor form factors. We showed that the leading orders of the 1/M -expansion are model independent and can be expressed in terms of known light flavor nucleon form factors. For the sub-leading terms one needs the calculation of the nucleon matrix elements of the gluon operators of dimension six. We showed that the corresponding operators in the instanton vacuum are suppressed by the instanton packing fraction. However, the contribution of the gluon operators can still be numerically important for the case of charm quark.
The heavy quark mass expansion eqs. (5.1), (5.2) can used to study the intrinsic heavy quark content of other hadrons, for example of vector mesons. The detailed studies of these topics will be published elsewhere. In addition, the methods developed in this paper can be generalized to the case of other fermionic operators containing heavy quark fields. For example, one can derive the heavy quark mass expansion of the twist-2 heavy-quark operators. This will provide direct access to the intrinsic heavy quark content in the nucleon structure functions. Acknowledgments M.V.P. acknowledges support of RSF grant 14-22-0281.
We thank Arseniy Filin for useful remarks and suggestions concerning the preparation of this paper. there is a set of Feynman diagrams to compute. By changing the order in which the gluons couple to the fermion loop, several topologically different diagrams are generated from each basic type of diagram. Each of these diagrams had to be computed separately. However, all integrals in one class can be computed with the same methods. There are three different integrals of the type shown in figure 5a, three different integrals of that shown in 5b, six integrals of that shown in 5c and ten different integrals of the type shown in figure 5d. Additionally, all permutations of the external gluons have to be considered. This was achieved by symmetrising the results in the sets of external momenta, color-and Lorentz-indices. We now list the analytic expressions corresponding to the types of diagrams in figure 5, denoting the three-and four-gluon-vertices as V µντ abc (l 1 , l 2 , l 3 ) = g f abc g µν (l 1 − l 2 ) τ + g ντ (l 2 − l 3 ) µ + g τ µ (l 3 − l 1 ) ν , (A.1) W µντ ρ abcd = −ig 2 f abe f cde (g µτ g νρ − g µρ g ντ ) + f ace f bde (g µν g τ ρ − g µρ g ντ ) + f ade f bce (g µν g τ ρ − g µτ g νρ ) . (A.2)