Application of the operator product expansion to the short distance behavior of nuclear potentials

We investigate the short distance behavior of nucleon-nucleon (NN) potentials defined through Bethe-Salpeter wave functions, by perturbatively calculating anomalous dimensions of 6-quark operators in QCD. Thanks to the asymptotic freedom of QCD, 1-loop computations give certain exact results for the potentials in the zero distance limit. In particular the functional form of the S-state central NN potential at short distance $r$ is predicted to be a little weaker than $r^{-2}$. On the other hand, due to the intriguing character of the anomalous dimension spectrum, perturbative considerations alone can not determine whether this potential is repulsive or attractive at short distances. A crude estimation suggests that the force at short distance is repulsive, as found numerically in lattice QCD. A similar behavior is found for the tensor potential.


Introduction
In a recent paper [1] a proposal has been made to study nucleon-nucleon (NN) potentials from a first principle QCD approach. In this field theoretic framework, potentials are obtained through the Schrödinger operator applied to Bethe-Salpeter (BS) wave functions defined by where |2N, E is a QCD eigenstate with energy E (suppressing here other quantum numbers), and N is a nucleon interpolating operator made of 3 quarks such as N (x) = ǫ abc q a (x)q b (x)q c (x). Such wave functions have been measured through numerical simulations of the lattice regularized theory [1,2,3,4]. Although many conceptual questions remain to be resolved, the corresponding potentials indeed qualitatively resemble phenomenological NN potentials which are widely used in nuclear physics. The force at medium to long distance (r ≥ 2 fm) is shown to be attractive. This feature has long well been understood in terms of pion and other heavier meson exchanges. At short distance, a characteristic repulsive core is reproduced by the lattice QCD simulation [1]. No simple theoretical explanation, however, exists so far for the origin of the repulsive core. For an approach based on string theories, see ref. [5]. By writing 2N, nπ, E|f n ( r, E) , (1.2) where |2N, nπ, E is a state with the energy E containing two nucleons and n pions (and/or nucleon-antinucleon pairs), we see that ϕ E ( r) = f 0 ( r, E). (Our normalization is 2N, nπ, E|2N, n ′ π, E ′ = 2Eδ nn ′ δ(E − E ′ ).) We may thus interpret the wave function ϕ E ( r) as an amplitude to find the QCD eigenstate |2N, E in N ( x + r, t)N ( x, t)|0 . The behavior of the wave functions ϕ E ( r) at short distances (r = | r|) are encoded in the operator product expansion (OPE) of N ( x + r, t)N ( x, t). An OPE analysis [6] of BS wave functions in the case of a toy model, the Ising field theory in 2-dimensions, successfully described the analytically known behavior. In this case the limiting short distance behavior of the potential does not depend on the energy (rapidity) of the state, and further it only mildly depends on energy (for low energies) at distances of the order of the Compton wave length of the particles.
In this report we perform an operator product expansion (OPE) analysis of NN BS wave functions in QCD, with the aim to theoretically better understand the repulsive core of the NN potential, (at least that of the measured BS potential). Thanks to the property of asymptotic freedom of QCD the form of leading short distance behavior of the coefficent functions can be computed using perturbation theory. A short summary of our results has been published in ref. [7] 1 .
In sect. 2 we start with some general considerations on BS potentials, and sect. 3 presents some standard renormalization group considerations. The anomalous dimensions of 3-and 6-quark operators are computed in sect. 4. Finally in sects. 5, 6 we discuss the application of the results to NN potentials. In appendix C we make a similar analysis for the I = 2 two pion system. For the convenience of the reader we give a brief summary of our results here. The OPE analysis shows that the NN central potential at short distance behaves as for the S-state (L = 0) , where S and I are total spin and isospin of the NN system, respectively, and β SI is negative and explicitly obtained as 4) where N f is the number of quark flavors, and the overall coefficient C E depends on the energy E. Unfortunately the OPE analysis is not as conclusive as that in the toy model referred to above, in particular the sign of C E is not determined by perturbative cosiderations alone. The latter requires additional non-perturbative knowledge of matrix elements of composite operators. A crude estimation using the non-relativistic quark model indicates that C E is positive, which implies a repulsive core with a potential diverging a little weaker than the generically expected r −2 at short distances.

Operator Product Expansion and potentials at short distance in 3 dimensions
In this section we discuss the application of the operator product expansion (OPE) to the determination of the short distance behavior of the BS potential. We consider the equal time Bethe-Salpeter (BS) wave function defined by where |E is an eigen-state Note that r dependence appears solely in D C AB ( r) while the E dependence exists only in 0|O C ( 0, 0)|E . As we will see, in the r = | r| → 0 limit, the coefficient function behaves as where θ, φ are angles in the polar coordinates of r, so that We now assume that C has the largest contribution at small r: for ∀ C ′ = C. The potential can be calculated from this wave function. As will be seen later, α C = α C ′ = 0 for the NN case in QCD. Furthermore states with zero orbital angular momentum (L = 0) dominates in the OPE, so that the wave function at short distance is given by we obtain the following classification of the short distance behavior of the potential.
2. β C = 0: In this case we have The sign of the potential at short distance depends on −β C ′ D C ′ (E)/D C (E).
If there are two or more operators which have the largest contribution at short distance, we have In this case, the above analysis can be applied just by replacing On the lattice, we do not expect divergences at r = 0 due to lattice artifacts at short distance. The above classification holds at a ≪ r ≪ 1 fm, while the potential becomes finite even at r = 0 on the lattice.

Renormalization group equation for composite operators
In QCD, using dimensional regularization in D = 4 − 2ǫ dimensions, bare local composite operators O (3.1) (Summation of repeated indices is assumed throughout this paper.) The meaning of the above formula is that we obtain finite results if we insert the right hand side into any correlation function, provided we also renormalize the bare QCD coupling g 0 and the quark and gluon fields. For example, in the case of an n-quark correlation function with operator insertion, which we denote by G (0) n;A (g 0 , ǫ) (suppressing the dependence on the quark momenta and other quantum numbers) we have We recall from renormalization theory that for the analogous n-quark correlation function (without any operator insertion) we have where the coupling renormalization is given by The renormalization constant Z 1 in the minimal subtraction (MS) scheme we are using has pure pole terms only: where Similarly for the fermion field renormalization constant, we have where γ F 0 is given by (4.10). The gluon field renormalization constant is also similar, but we do not need it here. Finally the matrix of operator renormalization constants is of the form The renormalization group (RG) expresses the simple fact that bare quantities are independent of the renormalization scale µ. Introducing the RG differential operator the RG equation for n-quark correlation functions can be written as where the RG beta function is where β D (g, ǫ) is the beta function in D dimenions and the RG gamma function (for quark fields) is It is useful to introduce the RG invariant lambda-parameter Λ by taking the Ansatz Λ = µ e f (g) (3.13) and requiring DΛ = 0. The solution is the lambda-parameter in the MS scheme (Λ MS ) if the arbitrary integration constant is fixed by requiring that for small coupling (3.14) Finally the RG equations for n-quark correlation functions with operator insertion are of the form (3.16)

OPE and RG equations
Let us recall the operator product expansion (2.2) We will need it in the special case where the operators O 1 , O 2 on the left hand side are nucleon operators and the set of operators O B on the right hand side are local 6-quark operators of engineering dimension 9 and higher. All operators in (3.17) are renormalized ones, but from now on we suppress the labels (ren) . As we will see, the nucleon operators are renormalized diagonally as and we can define the corresponding RG gamma functions by Next we write down the bare version of (3.17) (in terms of bare operators and bare coefficient functions): (3.20) Comparing (3.17) to (3.20), we can read off the renormalization of the coefficient functions: and using the µ-independence of the bare coefficient functions we can derive the RG equations satisfied by the renormalized ones: where the effective gamma function matrix is defined as

Perturbative solution of the RG equation and factorization of OPE
We want to solve the vector partial differential equation (3.22) and for this purpose it is useful to introduceÛ AB (g), the solution of the matrix ordinary differential equation and its matrix inverse U AB (g). We will assume that the coefficient functions are dimensionless and have the perturbative expansion where r = |y|. For the case of operators with higher engineering dimension 9 + α the coefficients are of the form r α times dimensionless functions and the analysis is completely analogous and can be done independently, since in the massless theory operators of different dimension do not mix. (In the full theory quark mass terms are also present, but they correspond to higher powers in r and therefore can be neglected.) We will also assume that the basis of operators has been chosen such that the 1-loop mixing matrix is diagonal:γ (3.26) In such a basis the solution of (3.24) in perturbation theory takes the form where R AB (g) = O(g 2 ), with possible multiplicative log g 2 factors, depending on the details of the spectrum of 1-loop eigenvalues β A . Having solved (3.24) we can write down the most general solution of (3.22): Here the vector F A is RG-invariant. Introducing the running couplingḡ as the solution of the equation F B can be rewritten as Since, because of asymptotic freedom (AF), for r → 0 alsoḡ → 0 as Putting everything together, we find that the right hand side of the operator product expansion (3.17) can be rewritten: There is a factorization of the operator product into perturbative and non-perturbative quantities: F B (Λr) is perturbative and calculable (for r → 0) thanks to AF, whereas the matrix elements ofÕ B are non-perturbative (but r-independent). An operator O B first occurring at ℓ B -loop order on the right hand side of (3.17) and corresponding to normalized 1-loop eigenvalue β B has coefficient F B (Λr) with leading short distance behavior (3.34) In principle, an operator with very large β B , even if it is not present in the expansion at tree level yet, might be important at short distances. This is why it is necessary to calculate the full 1-loop spectrum of β B eigenvalues. As we shall see, no such operators exist in our cases, and therefore operators with non-vanishing tree level coefficients are dominating at short distances. The corresponding coefficient functions have leading short distance behavior given by

OPE of two nucleon operators at tree level
The general form of a gauge invariant local 3-quark operator is given by where α, β, γ are spinor, f, g, h are flavor, a, b, c are color indices of the (renormalized) quark field q. The color index runs from 1 to N = 3, the spinor index from 1 to 4, and the flavor index from 1 to N f . In this paper a summation over a repeated index is assumed, unless otherwise stated. Note that B f gh αβγ is symmetric under any interchange of pairs of indices (e.g. B f gh αβγ = B gf h βαγ ) because the quark fields anticommute. For simplicity we sometimes use the notation such as F = f gh and Γ = αβγ as indicated in (4.1).
The nucleon operator is constructed from the above operators as where P +4 = (1 + γ 4 )/2 is the projection to the large spinor component, C = γ 2 γ 4 is the charge conjugation matrix, and τ 2 is the Pauli matrix in the flavor space (for N f = 2) given by (iτ 2 ) f g = ε f g . Both Cγ 5 and iτ 2 are anti-symmetric under the interchange of two indices, so that the nucleon operator has spin 1/2 and isospin 1/2. Although the explicit form of the γ matrices is unnecessary in principle, we find it convenient to use a (chiral) convention given by As discussed in the previous section, the OPE at the tree level (generically) dominates at short distance. The OPE of two nucleon operators given above at tree level becomes Knowing the anomalous dimensions of the 6-quark operators appearing in the OPE, which will be calculated later in this section, the OPE at short distance (r = | y| ≪ 1, where the coefficient functions behave as and β A,B,C are related to the anomalous dimensions of the 6-quark operators O f g,A αβ , of those with two derivatives O f g,B αβ,kl and of those with one derivative O f g,C αβ,k . The wave function defined through the eigenstate |E is given by for the anti-symmetric states, while for the symmetric states.
In this paper, we consider only 6-quark operators without derivatives and calculate the corresponding anomalous dimensions.

General formula for the divergent part at 1-loop
Following the previous section, we define the renormalization factor Z X of a k-quark operator X = [q k ] through the relation where q 0 (q) is the bare (renormalized) quark field. The wave function renormalization factor for the quark field is given at 1-loop by where λ is the gauge parameter and C F = N 2 −1 2N . At 1-loop the renormalization of simple k-quark operators (those involving no gauge fields) is given by the divergent parts of diagrams involving exchange of a gluon between any pair of quark fields. The 1-loop correction to the insertion of an operator q a,f α (x)q b,g β (x) in any correlation function involving external quarks is expressed as the contraction of where tr T A T B = δ AB /2 in our normalization. Since two identical contributions cancel the 2! in the denominator, the contraction at 1-loop is given by where the free quark and gauge propagators are given in momentum space as The above contribution can be written as where whose divergent part is independent of the momenta p, q and is given by We then obtain the divergent part of the 1-loop contribution as where (bold-faced symbols represent matrices in flavor and spinor space) Here we use the notation (4.21) Using the following Fierz identities for spinor indices where P R , P L are the chiral projectors we can simplify T 0 as

Renormalization of local 3-quark operators at 1-loop
In this subsection we calculate the anomalous dimensions of general 3-quark operators at 1-loop. In terms of the renormalization factor defined as λ ) is the λ-independent (dependent) part at 1-loop, the divergent part of the insertion of the 3-quark operator B F Γ = B f gh αβγ defined in (4.1) at 1-loop is given by a linear combination of insertion of baryon operators, and (with a slight abuse of notation) we express this as The λ-dependent contribution from T 1 in (4.19) is diagonal and given by so that the λ-dependent part of ζ vanishes: Therefore ζ is λ-independent, as expected from the gauge invariance. We remark that we leave N explicit in some formulae to help keep track of the origin of the various terms, but in our case we should always set N = 3 at the end. The λ-independent part of Γ (1) from T 0 in (4.25) leads to (N = 3): where α, β, γ ∈ {1, 2} (right-handed), whileγ ∈ {1 = 3,2 = 4} (left-handed). Note that the same results hold with hatted and unhatted indices exchanged. These relations can be easily diagonalized and the combinations which do not mix are given by where d is given by (4.36) The square bracket denotes antisymmetrization [αβ] = αβ − βα, and curly bracket means {αβ} = αβ + βα, {ααβ} = ααβ + αβα + βαα. The totally symmetric case corresponds to the decuplet representation (for N f = 3) and contains the N f = 2 , I = 3/2 representation. The antisymmetric case corresponds to the octet representation (for N f = 3) and contains the N f = 2 , I = 1/2 representation. The anomalous dimension at 1-loop is easily obtained from (4.37) Therefore we have

Anomalous dimensions of 6-quark operators at 1-loop
In this subsection we consider the renormalization of arbitrary local gauge invariant 6quark operator of (lowest) dimension 9. Any such operator can be written as a linear combination of operators with A = (Γ 1 , F 1 ) and B = (Γ 2 , F 2 ). Note O A (x) and/or O B (x) may not be operators with proton or nucleon quantum numbers and separately may not be diagonally renormalizable at one loop. The reason for considering the renormalization in more generality is that in principle there may be operators in this class which occur in the OPE of two nucleon operators at higher order in PT, but are relevant in the analysis because of their potentially large anomalous dimensions.

Linear relations between 6-quark operators
According to the considerations in subsect. 4.2 the operators in eq. (4.42) mix only with Note however that such operators are not all linearly independent. Relations between them follow from a general identity satisfied by the totally antisymmetric epsilon symbol which for N labels reads (4.43) For our special case, N = 3, this identity implies the following identities among the 6-quark operators Note that the interchange of indices occurs simultaneously for both Γ 1 , Γ 2 and F 1 , F 2 in the above formula. The plus sign in (4.44) appears because the quark fields are Grassmann. An immediate consequence of the identity is that the divergent part of the λ-dependent contributions, calculated from T 1 in (4.17), must vanish, after the summation over the 9 different contributions from quark pairs on the different baryonic parts A, B is taken. The λ-dependent part of the contribution of quark contractions on the same baryonic parts is compensated by the quark field renormalization. Thus the renormalization of the bare 6-quark operator is λ-independent as expected from gauge invariance.
As an example of identities, we consider the case that Γ 1 , Γ 2 = ααβ, αββ (α = β and where minus signs in the first line come from the property that There are no further relations among 6-quark operators beyond (4.44).

Divergent parts at 1-loop
We thus need only to calculate the contributions from T 0 , which can be classified into the following 4 different combinations for a pair of two indices: The computation can be made according to the following steps: i.) Select the total flavor content e.g. 3f + 3g or 4f + 2g (f = g). These are the only cases we will consider since in this paper we are mainly restricting attention to baryon operators with N f = 2, but the approach is also applicable to more general cases (N f > 2).
ii.) Given a flavor content classify all the possible sets of Dirac labels in the chiral basis e.g. 111223, 112234, ... It is obvious from the rules above that some have equivalent renormalization at 1-loop e.g. 111223 and 112223 with 1 ↔ 2, and also those with hatted and unhatted indices exchanged e.g. 111223 and 133344.
iii.) For given flavor and Dirac sets generate all possible operators 3 . Then generate all gauge identities between them and determine a maximally independent (basis) set {S i }.
iv.) Compute the divergent parts of the members of the independent basis: v.) Finally compute the eigenvalues and corresponding eigenvectors of γ T to determine the operators which renormalize diagonally at 1-loop.
An example of the procedure is given in Appendix A. Some of the steps are quite tedious if carried out by hand. e.g. in the case 3f + 3g and Dirac indices 112234 there are initially 68 operators in step iii. with 38 independent gauge identities, and hence an independent basis of 30 operators. However all the steps above can be easily implemented in an algebraic computer program using MATHEMATICA or MAPLE.
If the quarks f, g belong to an iso-doublet e.g. we identify f with u having I 3 = 1/2 and g with d (I 3 = −1/3), then if an eigenvalue is non-degenerate the corresponding eigenvector belongs to a certain representation of the isospin group. If the eigenvalue is degenerate then linear combinations of them belong to definite representations. For the 3f + 3g case they can have I = 0, 1, 2, 3. Eigenvectors with I = 0, 2 are odd under the interchange f ↔ g and those with I = 1, 3 are even. The operators in the case 4f + 2g have I 3 = 1 and hence have I = 1, 2, 3. The eigenvectors in this case can be obtained from those of the 3f + 3g case by applying the isospin raising operator.
The complete list of eigenvalues and possible isospins are given in Tables 2-4 in Appendix A. Here we summarize the most important results.
1) For the 3f + 3g (and 4f + 2g) cases all eigenvalues γ j ≤ 48d = 2γ N , where γ N is the 1-loop anomalous dimension of the nucleon (3-quark) operator. We have not found an elegant way of proving this other than computing all cases explicitly.
2) It is easy to construct eigenvectors with eigenvalue 2γ N e.g. operators of the form since there is no contribution from diagrams where the gluon line joins quarks in the different baryonic parts.
3) Operators with higher isospin generally have smaller eigenvalues.

Results for anomalous dimensions
It is very important to note here that operators B SI V I for both cases (SI = 01 and 10) have the maximal anomalous dimension at 1-loop, since as noted in point 2) above, no 1-loop correction from T 0 joining quarks from the two baryonic components exists for B F 1 ,F 2 αβγ,α ′β′γ′ type of operators. Therefore we always have some operators with β A = 0 which dominate in the OPE at short distance.
The 1-loop corrections Γ (1) to 6-quark operators B SI are computed in appendix A and are summarized as: for SI = 01. The last two results can be written as where Similarly we have for SI = 10 where Denoting the eigenvalues of the anomalous dimension matrix by γ C , we give the values of γ SI defined by

Short distance behavior of the nucleon potential
We consider the following structure of the potential.
where r = | y|, and is the tensor operator. Here σ i acts on the spin labels of the i th nucleon.
Since, as shown in the previous section, 6-quark operators appeared at tree level in the OPE of NN which have the largest and the second largest anomalous dimensions, we mainly consider 6-quark operators at tree level in the OPE, which is written as where c V I and c II are some constants, and · · · represents other contributions, which are less singular than the first two at short distance. (We here write the spinor and flavor indices α, β and f, g explicitly for later use.) The anomalous dimensions β SI 0 are given in (4.73) and (4.74).

Potential for S = 1 and I = 0
We here consider the spin-triplet and isospin-singlet state(S = 1 and I = 0). Since I = 0 and I z = 0 in this case, we drop indices I and I z unless necessary. In this case, the leading contributions in eq.(5.3) couple only to the J = 1 state, which is given by where | 3 S 1 , J z = 1 = |L z = 0, S z = 1 L=0,S=1 , (5.13) , (5.14) x is the mixing coefficient, which is determined by QCD dynamics, and 2S+1 L J specifies quantum numbers of the state. Relevant matrix elements are given by for i = II and V I, where B 0 i are non-perturbative constants, and Using the above results, we have By applying ∇ 2 we obtain where we use the formula in the appendix B, V 10 c (r) = V 0 0 (r) + V 0 σ (r), and with {αβ} −1 = δ α2 δ β2 and {αβ} 0 = (δ α1 δ β2 + δ β1 δ α2 )/ √ 2. By comparing eq.(5.18) with eq.(5.19), we have This shows that the central potential V 10 c (r) diverges as F 10 (r) in the r → 0 limit, which is a little weaker than 1/r 2 , while the tensor potential V T (r) becomes zero in this limit at the tree level in the OPE.

Higher order in the OPE and the tensor potential
While the short distance behavior of the central potential is determined by the OPE at the tree level, the determination of the tensor potential at short distance requires the OPE at higher order, whose relevant contribution is given by where [y k y l ] = y k y l − r 2 δ kl /3, and we assume that the third term with the tensor-type operator B f g,kl T,αβ appears first at ℓ SI T (> 0) loop of the perturbative expansion in the OPE. Therefore, with this assumption, we have where ∆ SI T = γ T /(2d) with γ T being the anomalous dimension of the operator B f g,kl T,αβ . The calculation of anomalous dimensions for all 6-quark operators in the previous section shows that ∆ T ≤ 24, so that β SI T < β SI 0 < 0. An extra matrix element we need is given as where B T is a further non-perturbative constant.
Using the above results, we have for (S, I) = (1, 0) By applying ∇ 2 we obtain This shows that the central potential V c (r) diverges as F 10 (r) in the r → 0 limit, which is a little weaker than 1/r 2 , while the tensor potential V T (r) diverges as F 10 T (r) in this limit, which is not stronger than F 10 (r) since β 10 0 − 1 ≥ β 10 T .

Evaluation of matrix elements
We rewrite 3-quark operators in terms of left-and right-handed component: for X, Y = R or L. In terms of these we have Note that we take ( x, t) = ( 0, 0) in the above operators. We need to know 0|(B i ) f g αβ |2N, E (6.10) for i = II, V I. For f = g, Lorentz covariance leads to where s = E 2 = 4 p 2 + m 2 N with the total energy E in the center of mass frame, σ i (i = 1, 2) is the spin of the i-th nucleon, and C AB XY is an unknown function of s. Using invariance of QCD under the parity transformation P B X P −1 = γ 4 BX whereR = L and L = R, we rewrite eq.(6.11) as where γ 4 u(− p, σ 1 ) = u( p, σ 1 ) is used. The above relation implies CĀB XȲ = C AB XY . Using this property for the unknown functions C AB XY , we have Taking p = (0, 0, p z > 0) and Dirac representation for γ matrices [9], we have where E N = p 2 + m 2 N . For I = 1 ( f g + gf ) and S = 0 (σ 1 = + and σ 2 = − ) the above explicit form for the spinors gives while, for I = 0 ( f g − gf ) and S = 1 (σ 1 = + and σ 2 = + ), we have We finally obtain for f g + gf and (σ 1 , σ 2 ) = (+, −) ( 1 S 0 ), and for f g − gf and (σ 1 , σ 2 ) = (+, +) ( 3 S 1 ), where s = 4E 2 N . Unfortunately, we can not determine the sign of the ratio for these matrix elements. As a very crude estimation, we consider the non-relativistic expansion for constituent quarks whose mass m Q is given by m Q = m N /3. In the large m Q limit, γ 4 q 0 = q 0 and γ 4 u 0 = u 0 , where a subscript 0 for q and u means the 0-th order in the non-relativistic expansion. In this limit, it is easy to show C AB XY = C for all X, Y, A, B, so that C RR RL+LR = 2C and C RL RL = C. Furthermore the first order correction to C AB XY = C vanishes in the expansion. Therefore in the leading order of the non-relativistic expansion, we have for (σ 1 , σ 2 ) = (+, +) ( 3 S 1 ). For both cases, we have positive sign for the ratio, which gives repulsion at short distance, the repulsive core.

Conclusions and discussion
The OPE analysis leads to conclusion that the S-state potential at short distance behaves as in (1.3) with (1.4). However perturbative considerations alone can not tell the crucial sign of the overall coefficient C E . Moreover we found that the latter was also not directly predicted by chiral PT . A crude estimation using non-relativistic quarks suggests that C E is positive, hence predicting a repulsive core, which diverges a little weaker than r −2 at small r. The leading corrections involve small powers of logs and hence it could happen that the dominant asymptotic behavior appears only at extremely short distances. Our analysis suggests that the repulsion of the NN potential at short distance is related to the difference of anomalous dimensions between a 6-quark operator and two 3-quark operators at 1-loop, and to the structure of the composite operators which probe the NN states. The explicit 1-loop calculation indicates that a combination of fermi statistics for quarks and the particular structure of the one gluon exchange interaction determines the sign and size of the β's. The appearance of zero effective gamma eigenvalues is simply explained by chiral symmetry, however we were unable to find a simple proof of the absence of positive eigenvalues established by explicit calculation. One speculation is that the latter intriguing pattern can be explained by the relation between the 1-loop QCD anomalous dimensions and those of super YM theory [10].
At higher order in the perturbative expansion, tensor operators appear in the OPE. Using this fact, we also found that the tensor potential also diverges a little weaker than r −2 as r → 0.
There are several interesting extensions of the analysis using the OPE. An application to the 3-flavor case may reveal the nature of the repulsive core in the baryon-baryon potentials. Since quark masses can be neglected in our OPE analysis, the calculation can be done in the exact SU(3) symmetric limit. It is also interesting to investigate the existence or the absence of the repulsive core in the 3-body nucleon potential. Such an investigation would require the calculation of anomalous dimensions of 9-quark operators at 2-loop level. Certainly more precise evaluations (also involving numerical simulations) of matrix elements 0|O X |E will also be needed to theoretically predict the nature of the core of the NN potential.

A. Explicit calculations of the divergent part for 6-quark operators at 1-loop
In determining the divergent parts for 6-quark operators we will consider the various cases in turn adopting the following mechanical (if rather inelegant) procedure. We first list the various operators which can appear and determine their linear relations. Then we compute the divergent parts using (4.25), initially keeping N explicit in the formulae to indicate from which part of (4.25) the terms originate. Finally we set N = 3 and use the constraint equations to express the result in terms of linearly independent operators.
A.1 S = 0 and I = 1 case We prepare following 11 operators which are all the possible operators corresponding to the Dirac labels 111222 in Table 2 (in the 4f 2g case). In terms of these we write Taking B 1,2,3,4 as independent operators, the constraints from (4.44) are: The 1-loop corrections Γ (1) to B 1,2,3,4 can be calculated from T 0 as where an overall factor d/ǫ is dropped on the rhs for simplicity. i.e. using the gauge identities eq. (A.3) we obtain with the matrix γ given by  Table 2.
The other eigenvalues of γ are obtained similarly.

B 01
We prepare the following 6 operators, , (A.9) in terms of which we have The 1-loop corrections to B 1,2,3,4 can be calculated as (A.11) Therefore we have , we introduceB i = B i (α ↔α), so that we have We prepare the following 15 operators: in terms of which we write (A.14) Taking B 1,2,3,4,5 as independent operators, others can be expressed as The 1-loop corrections to B 1,2,3,4 can be calculated as We therefore obtain We have 3 more structures: and B f f g . We introduce for the first one,B i = B i (α ↔α) for the second, and B i = B ′ i (α ↔α) for the third. We then have Γ 01 We prepare the following 12 operators: in terms of which we write The constraint from gauge invariance leads to The 1-loop corrections can be calculated as (A.21) Therefore we have , so that we have and B 01 In this case we have to prepare 29 operators as follows: in terms of which we have The constraint from gauge invariance gives where B 1,2,3,4,5,6,7,8 , X 1 , Y 1 , Z 1,2 and V 1 are taken to be independent. The 1-loop corrections to B 1,2,3,4 can be calculated as On the other hand, constraints give We prepare the following 11 operators: in terms of which Constraints are given by The 1-loop corrections Γ (1) can be calculated as and therefore we have , the result can be obtained from the above by the interchange of α and α. In terms ofB i = B i (α ↔α), we have We prepare the following 6 operators: B ggf ααβ , The 1-loop corrections can be calculated as and therefore we have , we introduceB i = B i (α ↔α), so that we have We prepare the following 15 operators: in terms of which we write Constraints are given by The 1-loop corrections can be calculated as . We introduce for the first class,B i = B i (α ↔α) for the second, and B i = B ′ (α ↔α) for the 3rd. We then have We prepare the following 16 operators: in terms of which we write Constraints give The 1-loop corrections can be calculated as We prepare the following 30 operators: in terms of which we have Constraints give The 1-loop corrections can be calculated as We introduce operatorsB i = B i (α ↔α) for B 5,6 , so that we have

B.1 Eigenstates
At given J, there are 2 distinct states, the spin-singlet (S = 0) state and the spin-triplet (S = 1) state.
The singlet state is denoted as 1 J J , since it has S = 0 and J = L. The fact that I + L + S must be odd to satisfy fermion anti-symmetry gives I = 0 for odd J and I = 1 for even J. The eigenstate with J z can be easily obtained as where |J z , S z J,S = |J z J ⊗ |S z S . The spin-triplet state is classified into 3 types: 3 J J , 3 (J ± 1) J . For the first one, I = 0 (even J) or I = 1 (odd J), and vice versa for the other two types. By the Wigner-Eckart theorem, the matrix elements of the five operators do not depend on J z . Therefore it is enough to know eigenstates with J z = J only. Explicitly we have

B.2 Evaluation of each operator
Using these eigenstates, it is easy to see for 1 J J , 3 J J , 3 (J − 1) J and 3 (J + 1) J , respectively. For S 12 defined in (5.2) the results are more complicated due to the mixing between 3 (J − 1) J and 3 (J + 1) J . After a little algebra we obtain, C. The I = 2 2-pion system Here we consider the operator product expansion of two iso-vector pseudoscalar densities in QCD.

C.1 Anomalous dimensions
The local composite operators with π + π + quantum numbers in QCD with lowest dimension are 4-quark operators with dimension 6. There are 5 independent such (bare) scalar operators O 1 =dγ µ u ·dγ µ u +dγ µ γ 5 u ·dγ µ γ 5 u, where in this appendix we use the notation u = q u , d = q d and suppress explicit color and Dirac indices of the quark fields. There are also 3 independent such (bare) traceless tensors operators T µν 1 =dγ µ u ·dγ ν u +dγ µ γ 5 u ·dγ ν γ 5 u − and similarly for tensor fields. The results for the one-loop anomalous dimensions of the scalar fields can be found in [11]. The non-vanishing entries of the mixing matrix for the scalar case are We have extended the analysis of ref. [11] to the tensor case. Here we find The OPE for two π fields can be written in QCD as π(x)π(0) = α γ α (x)B α + · · · (C.25) Here π(x) is the field annihilating π + , B α are the renormalized doubly charged dimension 6 operators discussed in the previous section, x is spacelike (and for simplicity we assume its time component vanishes) and γ α (x) are c-number coefficient functions. [Here we "pretend" all operators are scalar, although in fact three of them are symmetric traceless tensors. Taking into account their tensor structure however does not change any of our conclusions here.] The dots stand for higher dimensional operators with less singular coefficient functions. The short distance asymptotics of the I = 2 wave function is given by Ψ(x) = 0|π(x)π(0)|2 ∼ α γ α (x)B α + · · · (C.26) where B α = 0|B α |2 (C.27) are the (energy dependent) matrix elements of the local operators. In QCD we can write the divergence of the axial current as is a (bare) quark bilinear field with π + quantum numbers and the quark mass parameter m 0 is the sum of the u and d quark masses. [Here we used the fact that the axial current, being partially conserved, has renormalization constant Z = 1. There is a subtlety in dimensional regularization where because of the presence of γ 5 the renormalization constant is not equal to unity. It is finite nevertheless and this does not alter our conclusions at the 1-loop level.] The (canonically normalized) pion field is defined as where Ω is a constant: (C.31) From this we see that the field Φ R renormalizes with the inverse of the mass renormalization constant.
The RG analysis of the pion-pion wave function goes along the same lines as in the main text for the nucleon-nucleon case. By inspecting the spectrum of anomalous dimensions we see that, again, only operators already present in the tree level expansion contribute to the leading short distance part of the wave function and even from this set we need only the operators with the largest anomalous dimensions. Since the coefficient of such an operator X

(R)
A is asymptotically proportional to (− ln r) Note that the ratio ψ 1 /ψ 0 may be energy dependent. We need to calculate this ratio (or at least its sign) nonperturbatively, to be able to determine whether the potential in this channel is attractive or repulsive. ChPT is not applicable for this problem since there are too many extra low energy constants characterizing the matrix elements of 4-quark operators and in the end the sign of this ratio is left undetermined. In the absence of a reliable non-perturbative method to calculate the above matrix elements we try to estimate them by inserting a complete set of states in the middle of the operator and truncating the sum after the 1-particle contribution. This is very similar in spirit to the vacuum insertion method [12], (oft rightly criticized) however surprisingly successfully applied to ∆S = 2 weak matrix elements in the past. In this approximation 0|(dΓ 1 u) · (dΓ 2 u)|2 ≈ 0|(dΓ 1 u)|1 1|(dΓ 2 u)|2 Thus the ratio ψ 1 /ψ 0 is positive in this naive approximation and the potential is repulsive, as indicated by the (quenched) lattice measurements [13].