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

The framework of the expansion by subgraphs is used to compute asymptotic expansions for the vector and the tensor currents in the limit of large quark masses. We use the results to obtain an estimate for the influence of heavy quarks on the nucleon electromagnetic and tensor form factors.


Introduction
The intrinsic heavy quark content of light hadrons is a fundamental property of QCD [1,2,3]. The influence of intrinsic heavy quark content on the structure of the nucleon can be investigated by several methods. Much information about the nucleon structure like the charge distribution and the anomalous magnetic moment can be obtained from form factors. These form factors parametrize the expectation values of bilinear fermion operators in a single-nucleon state. To describe the influence of intrinsic heavy quarks on the form factors one should therefore consider those bilinear fermion operators that contain heavy quark fields, it is natural to employ a large mass expansion for the investigation of such operators. One way to perform large mass expansions is the heavy quark mass expansion (HQME) which was used in [4] and [5] to obtain an expansion of the vector and the tensor current of heavy quarks in terms of gluon field operators.
One observation that the HQME-approach is not complete was made in [6]. This observation concerns the influence of muons on the anomalous magnetic moment µe of the electron in QED. If the muon mass is treated as large parameter (compared to the electron mass) this influence can be described by using the HQME to expand the muonic vector current in terms of photon field operators. The first photon field operators in the expansion that contribute to µe are of order 1/m 4 µ [7]. However, from the evaluation of the corresponding Feynman diagrams one can show [8] 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] it was argued that the problem in QED arises because the 1/m 4 µ -suppressed photon operator has to be renormalized by 1/m 2 µ -suppressed electron operators. That settles the problem in this particular case, but it would be desirable to have a systematic procedure to find those additional operators that complete the operator product expansions of the HQME.
To give such a systematic procedure, we use the language of Feynman diagrams instead of that of path integrals that was used to derive the HQME results. In this language, rigorous results about the influence of virtual heavy particles are known for a long time. In Ref. [9] the decoupling theorem was proved: 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. Afterwards much work was devoted to the problem of how to obtain 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 functions (see for example Ref. [10], Ref. [11]). The first treatments were done within a momentum subtraction scheme because here the close connection between asymptotic expansions and renormalization was most apparent. However, due to the fact that dimensional regularization is nowadays the preferred method for regularization both because of its technical simplicity and its preservation of gauge invariance, it was desirable to derive similar results within the minimal subtraction scheme. In the MSscheme 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]. Subsequently a systematic procedure for the asymptotic expansions within the MS-scheme was developed (see for example [13]). Finally the technique of 'expansion by subgraphs' [14] was developed, for which it was proved that one obtains an expansion that is free of infrared divergencies and that has the correct 1/M behaviour. This technique works on the diagramatic level in the first place. However, with the help of the counterterm technique these results could be immediately generalized to an operator product expansion.
Here we will use this technique to analyze the validity of the HQME-approach. This is done for the especially important examples of the vector current and the tensor current. The Section explains the basic idea of asymptotic expansions. In the following section the expansion is carried out for the vector current and the tensor current. Finally, we use the expansion of the vector current to estimate the influence of charm quarks on the e.m. and tensor form factors of the nucleon. The reader interested in the final result can directly go to Section 5. Applications for the intrinsic charm in nucleon form factors are discussed in Sections 6 and 7.

Asymptotic expansions in heavy quark masses
Perhaps the simplest way to illustrate the general idea is to cite the power counting theorem for the dependence of a Feynman diagram on the mass of a heavy particle from Ref. [15]: Let Γ be a diagram with l loops. The asymptotically irreducible diagrams (AI-diagrams) of Γ are those connected subdiagrams that contain heavy lines and cannot be made disconnected by cutting a single light line. Then the asymptotic behaviour of Γ is bounded by where S is that set of disconnected AI graphs (spinney) containing all heavy lines which has the highest degree of divergence.
From this it is immediately clear that by introducing an operator that lowers the degree of divergence of the AI-diagrams sufficiently (like the R-operation in momentum subtraction schemes lowers the degree of divergence to make 1PI graphs convergent) one can generate a remainder that is suppressed by an arbitrarily chosen power of the heavy mass. The difference between this remainder and the original (renormalized) diagram will then constitute the correct asymptotic expansion for the diagram. This procedure is systematically described in [13]. We only give the main result. We first have introduce some notation. We denote the large mass in which the expansion is done by M , the external momenta are collectively denoted by {p}. Let Γ be an arbitrary graph and let the diagram I Γ ({p}, m, M ) be the corresponding analytic expression. The degree of divergence of the diagram Iγ is ωγ . Let S AI (Γ) be the set of all AI-spinneys (set of mutually disjoint AI subgraphs) and I Γ/S the diagram I Γ with all elements of the spinney S shrunk to a point. I Γ/S • γ∈S Vγ is the diagram in which the vertex Vγ replaces the subgraph γ for each γ ∈ S. R un I Γ/S is 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 the diagram 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 1 M a+1 -terms, It was also proven that this result is free of artificial UV-and IR-divergences, that is, the divergences of this expression are the same as those for the original diagram. To summarize: To perform an asymptotic expansion in large masses on has to find all AIspinneys of a diagram Γ. In each subgraph γ of this spinney, one performs a Taylor expansion 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.
The transition from the diagrammatic to the operator level is now fairly obvious: If one computes the matrix element of an operator containing heavy particles one will obtain Feynman diagrams containing AI-subgraphs. The expansion of these diagrams in terms of the heavy mass amounts to replacing the AI subgraphs by their Taylor expansions in light masses and external momenta. The terms of this expansion serve as new effective vertices in the complete diagram. The complete dependence on the heavy mass is in these effective vertices. Thus, if we reexpress the effective vertices as vertices due to local operators in the light degrees of freedom, all matrix elements of the heavy operator can be described by matrix elements of light operators. Now, if one wants to expand the heavy operator O of dimension dim(O) to the order a, on the diagrammatic level one has to expand each AI-subgraph to the order aγ = ωγ + a.

Since in theories with dimensionless coupling constant
For example, at the order 1/M 2 for dim(O) = 3 only diagrams with E B ≤ 5 and E F ≤ 2 (remember that external fermion lines have to come in pairs) have to be taken into account. In gauge theories, Ward-identities will lower the actual degree of divergence for diagrams with external vector bosons. This reduces the number of diagrams that have to be considered. In this paper the vector current and the tensor current are considered. For these operators an OPE is done to the order 1/M 2 for the vector current and to the order 1/M for the tensor current and up to the order α 3 in the strong coupling constant. The reason for the order α 3 is the fact that this is the lowest order with purely fermionic operators.

The AI-diagrams for fermion-bilinears
We first give a general constraint for the structure of the expansion (see the appendix of Ref. [10]). For the argument that follows we use a regulator that admits the usual definition of γ 5 . Since the heavy particles should not occur as external particles, every AI-diagram will contain heavy fermions only inside closed loops. If we switch the sign of M , we have Now, the factors of γ 5 can be absorbed in the vertex factors via This cancels all γ 5 in the loop, leaving the diagrams with vector current insertion as they are and changing the sign of the diagrams with tensor current insertion. One can conclude that in the expansion of the vector current only even powers of M can occur while the expansion of the tensor current contains only 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 Matrix elements with one external gluon are zero by color conservation. Let us consider the case of two external gluons. Under charge conjugation, we have Thus, using global color conservation results in This shows that for the vector current and the 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. One can also find the only possible color structure for the matrix elements with three external gluons: From global color conservation we see that there are only two relevant color structures for the matrix element of the current with three external gluons (8 ⊗ 8 ⊗ 8 contains two different singlets). These can be taken as the two independent trace structures. In combination with charge conjugation invariance this leads to: Thus, For the case of the vector current one can conclude from QED current conservation that there is an additional factor of the momentum entering the diagram at the operator insertion. By dimensional analysis, the dimension in the heavy mass is lowered by by 1.
We will show in section 4.1.1 that if three gluons couple directly to the fermion loop each of the external gluons comes with a factor of external momentum. Taking into account both kinds of momentum factors, we have the power in the heavy mass ωγ = ωγ − E h g − δ where E h 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 ωγ = ωγ − δ. Now, we list the AI-diagrams relevant for us as shown in figure 1. We start with one loop diagrams. Following the arguments given above, the diagrams with three external gluons [1a] give the first contributions for the vector current and the tensor current. For the vector current these diagrams have ωγ = −4, while for the tensor current they have ωγ = −3. Actually these diagrams can be excluded from the computations at the order we consider. They are nevertheless, because we want to compare the results to the HQME-results given in [5]. The other group of diagrams that is relevant at the order considered in our expansion is that with four external legs [1b]. By power counting, these diagrams have ωγ = −2 for the vector current and ωγ = −1 for the tensor current. However, our explicit computation and the results of section 4.1.1 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 suppressed up to the 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 AI-diagrams that are relevant for us 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 external fermions. In the computation we will only consider one light quark flavor since all light flavors will, up to Figure 1: AI-diagrams for the HQME of vector current and tensor current: The blob denotes the operator insertion, the double lines denote heavy quarks trivial modifications, yield the same result. 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 Fig. 1k. For external ghost fields charge conjugation invariance shows the absence of contributions with two external ghosts (due to the ghost number conservation the number of ghosts fields must be even): All nonzero diagrams with more external ghost fields are of higher order than 1/M 2 in the expansion. Note that to all diagrams given in Fig. 1 one also has to compute the crossed diagrams.
There is one additional restriction for the purely gluonic operators in the expansion: The expansion should respect local gauge invariance. Thus, it should contain only field strength tensors and covariant derivatives of them. Since the behaviour of the field strength tensor under C is the same as that for the vector potential, the arguments given above show that the result for both operators must have the schematic structure trc F 1 , F 2 F 3 . For the vector current, current conservation requires that the result must be the derivative of some antisymmetric tensor, giving ∂ trc F 1 , F 2 F 3 . From power counting it is now clear 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.

Computation of the AI-diagrams
We can now use the relation (1) to compute the large mass expansion for the diagrams shown in the previous section. There is one additional simplification. For the expansion considered here single diagrams may have divergences, but if one takes the sum over all diagrams that contribute at a given order in the heavy mass and the coupling constant, the divergences cancel. Actually even the sums of all diagrams of a given type (as indicated in Fig. 1) are all 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 its small masses and external momenta. However, it is more convenient to rescale the integration momenta according to l → M l (note that this makes the integration variables dimensionless). This makes it obvious that we can also expand in 1/M .
After the Taylor expansion of the 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 done 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 the tensor integrals were reduced to scalar ones. Since the expansion is in all external momenta, the (tensor) integrals that we obtain effectively correspond to diagrams without external momenta; they are vacuum diagrams. For these, the tensor decomposition is easy: Only the metric tensor is available to represent the tensor structure of a diagram. Thus, an arbitrary integral F µ 1 ...µ 2n (all integrals with an odd number of indices vanish) occurring in the expansion can be written as 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. (2) such that the number of independent constants A I is further reduced. Now, contracting 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 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 from 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(gs)-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.
Now, the commutators are determined by the transformation properties of the operators under SU (3)c-transformations: Here we assume that possible Schwinger-terms can be neglected for vacuum expectation values that we are computing. This gives the result In momentum space we have Obviously, there are similar relations for all external momenta p 2 . . . pn. Now,we can draw conclusions for the case of three and four external gluons. For 3 external gluons, we have The last equality follows from C-invariance as was discussed above. From this relation (and from the corresponding ones for p 2 , p 3 ) it is clear that Here n = 3 for the tensor current and n = 4 for the vector current. Since this applies for all external momenta, those parts of the matrix elements that are not suppressed as described above must have at least four projection operators. From power counting and the general structure of the expansion for tensor current and vector current we see that the contributions with these projectors are suppressed by M −5 and M −6 , respectively.

Feynman integrals for the for the one-loop contribution
The only 1-loop integrals to compute are shown in 1a and 1b. The integrals were computed for a fixed assignment of external momenta. Afterwards the permutations in external momenta were done to generate the crossed diagrams. For the diagram 1a we have the analytic expression: Here we use Γ = γ µ for the vector current and Γ = σ µν for the tensor current. For the diagram 1b with four external lines we have The master integrals that occur after the expansion are: . We use l 3 = −q − l 1 − l 2 to simplify the notation (this is the momentum entering in the last gluon vertex).

Result for the vector current
First consider the vector current in 1a. The result of expanding these diagrams up to the order 1/M 4 reads: . This way to write the result clearly shows the expected projector structure of the result: Contracting with qµ, (l 1 )α, (l 2 ) β or (l 3 )τ we obtain zero.
Due to the projector-structure of the results it is easy to write down the corresponding operators. After a short calculation, we have, omitting the argument x and neglecting higher terms in the coupling constant: With the application to instanton solutions in mind, we introduce the field strength tensor Fµν = gsGµν . Thus, the contribution can be generated by This is exactly the result from the calculation of [5] after continuing to Minkowskispacetime and rescaling the gluon fields. The contribution of figure 1b will not be given here due to its huge size. However, it corresponds exactly to the O g 4 S -contribution to Eq. (3).

Tensor current
For the tensor current a similar computation gives the result: As expected the result vanishes for contractions with (l 1 )α, (l 2 ) β or (l 3 )τ . Again, we can rewrite this in terms of an operator insertion: Thus, we can write this contribution in terms of the operator insertion which coincides with the result of [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 1b corresponds to the O g 4 S -contribution to equation 4.

Two-loop contribution
Since there are several different kinds of integrals at two-loop level, their explicit forms are given in the Appendices A and B. The integrals that have to be computed can be done in closed form and are of the type As was discussed in section 3 there are no gauge invariant purely gluonic operators that can contribute to the order 1/M 2 for the vector current or 1/M 3 for the tensor current. Therefore, the 1/M 2 (1/M)-contributions of the 2-loop diagrams should add up to 0 at each order in perturbation theory. At O g 5 s and O g 6 s the results given in the Appendices A and B confirm this reasoning:

Three-loop contribution
In this case we have to compute the three loop integrals with external fermions (1k). These are given by Here m is the mass of any of the light quark flavors. The integrals that have to be computed are of the form These integrals are hard to compute directly, but there is a way of computing them by the integration by parts (IBP) technique. Here we give a slightly modified version of the approach of [19], where exactly these kinds of integrals where computed. One should also note that there are several software packages that compute such kinds of integrals ('MATAD' [20], 'FIRE' [21]). However, due to the fact that these integrals occur as an intermediate step after expanding in the heavy quark mass and tensor reduction of the result and that one has to do at the order of 50.000 of these integrals it was more convenient to write an own program that is specialized on exactly this type of integral.
These relations will be used such that the series of indices α 1 ,α 2 , α 3 is always decreasing unless one of the indices α 4 , α 5 ,α 6 is non-positive and α 1 , α 2 , α 3 are positive. In these cases the indices α 4 , α 5 ,α 6 will be chosen to be decreasing.
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: Afterwards, we show how to compute the integrals for these special cases.
Using the relations with the standard-notation one can use to obtain This reduces α 2 + α 4 + α 5 . The relation can be used until one of the indices is 0. 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 one finally arrives at α 4 = 0 or α 5 = 0 , 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 0. Then one can use the relation: I (α 1 , α 2 , α 3 , α 4 , α 5 , 0) . This can be used until α 2 or α 5 is 0. We conclude that in any case we end up with one of the special cases given above. Now we show how to compute the integrals for the special values of the coefficients.
The cases α 1 ≤ 0, α 2 ≤ 0 can be reduced to the case α 3 ≤ 0. In this case one has the integral If one of the indices α 5 or α 6 is nonnegative, the integral is 0. Expanding the numerator and doing the obvious cancellations, one is left with integrals of the type: After the some cancellations, one has integrals of the type If α 1 or α 2 are non-positive, the integral is 0. If all indices are positive, this is the standard-integral from the appendix. If α 3 is negative, one can actually use the result (5) as trivial analytic continuation to α 3 < 0.
α 5 ≤ 0 and α 6 ≤ 0: Renaming and exchanging integration variables, we have: If any of the massive lines have non-positive indices, the integral is 0. If α 4 is nonpositive, the integral reduces to a product of three integrals like d d k Doing the l 2 -integration and simplifying one has the easy integrals (take p < q): Here, Note that in this case all terms in the numerator can be canceled since all exponents are integers. This gives: This is again Eq. (5).

Results for the three-loop contributions to the vector current
The techniques described above could be used to compute the contributions for the vector current. We use αs = gs 4π to obtain: Here we use antisymmetrization including the symmetry factors, for example: q [µ γ ν] = 1 2 (q µ γ ν − q ν γ µ ). In this case it is possible to read off the corresponding operators immediately.

Results for the three-loop contributions to the tensor current
Here the computation yields the result: .
In terms of operators, we have: .
One can use the equations of motion and standard relations for the γ-matrices to show: The result for the tensor current can now be written as:

Result
In the final result, we have to remember that each light flavor will contribute to the three-loop diagrams. Collecting the contributions, we have Qσ µν Q = αs π Here f denotes the sum over all flavors of quarks that are light compared to M .
These operator identities can be used to analize intrinsic heavy quark content of light hadrons. One of immediate applications is the study of the intrinsic charm influence on the nucleon e.m. and tensor form factors.

Application: 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 nucleons 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. First of all, we need the form factor decomposition for the vector and the tensor current. The decomposition for the vector current is standard: The tensor current can be decomposed as [22]: We actually need qν p ′ , σ ′ Ψ iσ µν Ψ(0) |p, σ : Here 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 S µν = trc 7 5 F µγ , F γδ F δν + F δγ F δγ F µν . Since this operator is an antisymmetric tensor of rank 2, it can be decomposed as the tensor current. Thus, and Now, we consider the nucleon matrix element of the vector current made of heavy quarks.Using the derived operator identity (6) we obtain the heavy quark contribution to the electric form factor G Q q 2 and the magnetic form factor We see that if we neglect the chiral corrections of order ∼ m f /m N , than the leading This gives the following result for the heavy quark contribution to the electric and magnetic radii, which we define as r 2 : . 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 additional suppression by the light quark mass ∼ 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 (13) coincides with the result obtained in Ref. [6], the results for the e.m. radii (11,12) are new.
From Eqs. (11,12), neglecting the corrections of order ∼ m f /m N , one obtains the following model independent limits for the e.m. radii of the nucleon: We see that the limiting values of the heavy quark contribution to nucleon e.m. radii are always negative, so the intrinsic heavy quarks shrink the e.m. size of the nucleon. These limiting expression can be useful for lattice QCD simulations.  [26,27]. Here isospin symmetry was assumed and the values were evolved from a renormalization scale of µ 2 = 0.36GeV 2 to µ 2 = m 2 c by the method described in these publications. The numerical result for the charm quark contribution to the e.m. radii of the nucleon is: Here we have the same result for protons and neutrons 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 (10) are numerically of order unity, than 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 (10). 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 (10) is exactly zero on the instanton field. This implies that one has to consider the contribution of the instanton-anti-instanon pair to this operator, hence the gluon form factors (10) 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: Here the factor (Nc − 2) takes into account that the gluon operator (10) is zero for the case of Nc = 2, the factor πρ 2 R 2 ∼ 1 3 is the instanton packing fraction, and f (2) S ∼ 0.1 is the twist-4 contribution to the nucleon structure functions (see details in Refs. [4,30,31]).
Surely, the estimate (18) is very rough. Nevertheless, it indicates that suppression of the gluon forma 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 make more accurate calculations of the gluon form factors, we shall give detailed estimate of them elsewhere. Let us note obvious general properties of the gluon form factors (10): 1) they exactly zero for the case Nc = 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 light quark mass ∼ m f /m N , therefore we expect that for the case of the charm quark the subleading contribution ∼ 1/M 4 is dominant numerically. Indeed, if one take the values of the parameters discussed above, one obtains the following result: Clearly the contribution of the gluon form factor is dominant, therefore it is very important to make an estimate of this form factor.

Application: Heavy quark mass contribution to the tensor form factors
The nucleon tensor form factors are defined by Eq. (9). The heavy quark contribution to the tensor form factors can be easily read of the operator identity (7) we derived.
The corresponding equation contains the gluon operator O µν := trc F γν F ρµ Fγρ−F γµ F ρν Fγρ+ FγρF γρ F µν , which nucleon matrix element we parametrize as: packing fraction. However, still the contribution of the gluon operators can be important numerically for the case of charm quark.
The heavy quark mass expansion Eqs. (6,7) 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. Also the methods developed here can be easily generalized to the case of other heavy quark made fermionic operators. For example, one can derive the heavy quark mass expansion of the twist-2 heavy quark quark operators. This gives us direct access to the intrinsic heavy quark content in the nucleon structure functions. All diagrams were computed with a fixed assignment of external momenta. After the results were obtained, symmetrization in external momenta was done. Note that before symmetrization, the diagrams of Fig. 3b come with an additional factor of 1/2 due to the four-gluon-vertex. To simplify the expressions, we denote the three-and four-gluonvertices 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 ) ν 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 νρ .