Short distance repulsion in 3 nucleon forces from perturbative QCD

We investigate the short distance behavior of 3 nucleon forces (3NF) defined through Nambu--Bethe--Salpeter wave functions, using the operator product expansion(OPE) and calculating anomalous dimensions of 9--quark operators in perturbative QCD. As is the case of NN forces previously considered, we show that 3NF have repulsions at short distance at 1--loop, which becomes exact in the short distance limit thanks to the asymptotic freedom of QCD. Moreover these behaviors are universal in the sense that they do not depend on the energy of the NBS wave function for 3 nucleons.


Introduction
Realistic nuclear potentials between two nucleons (2N), determined precisely from 2N scattering data together with the deuteron binding energy, have often been used to study nuclear many-body problems. These two-nucleon forces (2NF), however, generally underestimate the experimental binding energies of light nuclei [1,2] and this fact indicates the necessity of taking into account three-nucleon forces (3NF). In addition, a clear indication of 3NF is observed in high precision deuteron-proton elastic scattering data at intermediate energies [3].
The 3NF may also play an important role for various phenomena in nuclear physics and astrophysics, which include (i) the backward scattering cross sections in nucleus-nucleus elastic scattering [4], (ii) the anomaly in the oxygen isotopes near the neutron drip-line [5], and (iii) the nuclear equation of state at high density relevant to the physics of neutron stars [6]. Universal short distance repulsion for three baryons (nucleons and hyperons) is also suggested to explain the observed maximum mass of neutron stars [7].
Despite of its phenomenological importance, a microscopic understanding of 3NF is still limited, due to difficulties in studying 3NF experimentally. Pioneered by Fujita and Miyazawa [8], the long range part of 3NF has been modeled by two-pion exchange [9], which is known to be attractive at long distance. In addition a repulsive component of 3NF at short distance is introduced in a purely phenomenological way [10].
To go beyond phenomenology, it is most desirable to determine 3NF directly from the fundamental degrees of freedom , the quarks and the gluons, on the basis of QCD. Recently the first investigation of this kind has been attempted using lattice QCD simulations, where 3NF have been extracted from the Nambu-Bethe-Salpeter (NBS) wave function for a specific alignment of 3 nucleons [11,12,13]. The method used there had been previously employed to extract nucleon-nucleon potentials (i.e. 2NF) [14,15,16,17] as follows. The NBS wave function for 2 nucleons is defined by ϕ E ( r) = 0|N ( x + r, t)N ( x, t)|2N, W , (1.1) where |2N, W is a QCD eigenstate with energy W = 2 m 2 N + k 2 with m N being the nucleon mass, E = k 2 /m N represents the non-relativistic kinetic energy, 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). The non-local but energy independent potential (more precisely the half off-shell T -matrix) is extracted from this NBS wave function as where H 0 = −∇ 2 /m N . The non-local potential can be expanded in terms of the velocity (derivative) with local function as at the lowest few orders, where r = | r|, σ i represents the Pauli-matrices acting on the spin index of the i-th nucleon, S = ( σ 1 + σ 2 )/2 is the total spin, L = r × p is the angular momentum, and is the tensor operator. This method has been shown to work well. The central potentials at the leading order in the expansion have qualitatively reproduced common features of phenomenological 2N potentials: the force is attractive at medium to long distance while it has a characteristic repulsive core at short distance. See also refs. [18,19] for a summary of results and recent developments. The present authors have investigated short distance behaviors of the 2NF defined in the framework mentioned above, using the operator product expansion(OPE) and perturbation theory thanks to asymptotic freedom of QCD [20,21,22] . (See also a similar attempt for the solvable models in 2 dimensions [23].) The behavior of the NBS wave function ϕ E ( r) at short distance (r → 0) is encoded in the operator product expansion (OPE) of the two nucleon operators: where {O k } is a set of local color singlet 6-quark operators with two-nucleon quantum numbers. Asymptotically the x-dependence and energy dependence of the wave function is factorized into Standard renormalization group (RG) analysis gives [21] the leading short distance behavior of the OPE coefficient function as where ν k is related to the 1-loop coefficient of the anomalous dimension of the operator O k , d k is the tree-level contribution of D k ( 0) and finally r 0 is some typical non-perturbative QCD scale. Assuming its matrix element does not vanish, the operator with largest RG power ν k dominates the wave function (1.7) at short distances. We denote the largest power by ν 1 and the second largest one by ν 2 . If ν 1 is non-zero, this leads to the leading asymptotics of the s-wave potential of the form which is attractive for ν 1 > 0 and repulsive for ν 1 < 0. If ν 1 = 0, the situation is more complicated. The relative sign of the ratio R between the leading and the subleading contributions becomes important and we find: If R is positive, the potential is repulsive, while it is attractive for negative R. A system of two nucleons corresponds to this degenerate case. Unfortunately in this case R depends on the energy E. In [21] it is argued that in the relevant energy range the relative coefficient R is positive, so that the short distance limit of the nucleon potential is repulsive.
In this paper, we extend the above OPE analysis to the 3NF. The corresponding equal time NBS wave function for 3 nucleons is given by where E 3N and |E 3N denote the energy and the 3N state. We introduce Jacobi coordinates From this wave function, the three nucleon potential is defined by where V 2N ( r ij ) with r ij = x i − x j denotes 2NF between (i, j)-pair, V 3N F ( r, ρ) the 3NF, µ r = µ ρ = m N /2 the reduced masses. In sect. 2, we start with renormalization group considerations and OPE, which are relevant for 3NF. The anomalous dimensions of 9-quark operators are computed in sect. 3.
Finally we discuss the short distance behavior of 3NF in sect. 4. For the convenience of the reader we give a brief summary of our results here. The OPE analysis shows that the 3N central potential at short distance behaves as 14) where N f is the number of dynamical quarks. Unlike the 2NF where the situation was not completely determined by PT alone, it is shown that the 3N potential always has a repulsive core. Furthermore it is universal in the sense that it does not depend on details of the 3N state used to define the NBS wave function such as the energy of the state.

Renormalization group equation for composite operators
In QCD, using dimensional regularization in Summation of repeated indices is assumed throughout this paper unless indicated otherwise. The meaning of the above formula is that we obtain finite results if we insert the right hand side into any correlation function of the fundamental gluon and quark fields, 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 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 and λ is the covariant gauge parameter. 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 Here the RG beta function is where β D (g, ǫ) is the beta function in D dimensions and the RG gamma function (for quark fields) is It is useful to introduce the RG invariant lambda-parameter Λ by taking the ansatz 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 (2.14) Finally the RG equations for n-quark correlation functions with operator insertion are of the form (2.16)

OPE and RG equations
Let us recall the operator product expansion for ρ, r → 0: We will need it in the special case where the operators O 1 , O 2 , O 3 on the left hand side are nucleon operators and the set of operators O B on the right hand side are local 9quark operators of engineering dimension 27/2 and higher. All operators in (2.17) are renormalized ones, but from now on we suppress the labels (ren) . As we know, the nucleon operators are renormalized diagonally as and we can define the corresponding RG gamma functions by For the nucleon operator, (2.20) Next we write down the bare version of (2.17) (in terms of bare operators and bare coefficient functions): Comparing (2.17) to (2.21), 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 (2.23) 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 has the perturbative expansion where s 2 = r 2 + ρ 2 and tan ω = ρ/r with r = | r|, ρ = | ρ|, and Ω r , Ω ρ are solid angles of the vectors r and ρ, respectively.
is the dimension of the coefficient function. Note that 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 ρ, 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:γ (2.27) In such a basis the solution of (2.25) in perturbation theory takes the form can be rewritten as Since, because of asymptotic freedom (AF), for s → 0 alsoḡ → 0 as can be calculated perturbatively using (2.26) and (2.28). Putting everything together, we find that the right hand side of the operator product expansion (2.17) can be rewritten: There is a factorization of the operator product into perturbative and non-perturbative quantities: F α 1 α 2 B (Λs, ω, Ω r , Ω ρ ) is perturbative and calculable (for s → 0) thanks to AF, whereas the matrix elements ofÕ B are non-perturbative (but s-independent). Note that d C = d B .

3NF from OPE
Using results in the previous subsection, the NBS wave function for 3N can be written at short distance as Since a r α 1 ρ α 2 term produces angular momenta l 1 ≤ α 1 and l 2 ≤ α 2 , we can write up to less singular terms at short distances. Then the sum of 2N and 3NF potentials V 2N +3N F is extracted as we have Since terms withd A = 0 dominate in Ψ 3N at short distance, contributions fromd A = 0 terms to 2N+3NF potentials are suppressed by an r α 1 ρ α 2 factor, so that they do not contribute at short distance. Therefore we consider terms withd A = 0 (α 1 = α 2 = 0) hereafter and do not write α i dependence in coefficients. We then have so that the ω dependent part gives a 1/s 2 contribution unless ω = 0, π/2, π, 3π/2, where either r = 0 or ρ = 0. We assume r = 0 and ρ = 0 in our calculation. Since an ω dependence appears only at 1-loop or higher orders in D A , we can neglect it unless an operator O A which appears at ℓ A loop has large anomalous dimension such that β A − ℓ A is larger than other β B corresponding to operators O B appearing at tree level. As we will see later such operators are absent; it is then enough to consider the tree level contribution in D A , so that ω dependent terms in D A can be neglected. The largest eigenvalue among operators appearing at tree level is thus denoted by β A , which corresponds to ν 1 discussed in the introduction for 2N forces.
We then obtain The NBS wave function is dominated by the term with largest β A . If we assume β A is non-zero, we finally obtain . (2.44)

OPE for 3N 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 c = 3, the spinor index from 1 to 4, and the flavor index from 1 to N f . 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 (3.1).
The usual nucleon operator which is employed in lattice simulations 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 an 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 3N operators given above at tree level becomes where + · · · denote higher dimensional operators, which do not contribute at short distance. The leading operators couple only to the states with L = 0 (the S-state).
If we construct local 3N operators at L = 0 from B f α (x) for nucleons (which only involve 2 different flavors), there is only one such operator, which has I = 1/2 and S = 1/2, due to the Pauli statistics of nucleons. Explicitly it is given by where f = g, f, g = u, d and α = β, α, β = 1, 2,1 = 3,2 = 4 for the explicit form of the γ matrices. Above no summation is taken for f and α.

General formula for the divergent part at 1-loop
As shown in ref. [21], the gauge invariant part of the divergence from diagrams involving exchange of a gluon between any pair of quark fields is given by where for either α 1 , β 1 ∈ {1, 2}(right-handed) or α 1 , β 1 ∈ {3, 4}(left-handed), and it vanishes for other combinations. The operator in eq. (3.5) can be written as a linear combination of simple operators According to this 1-loop formula for divergences, such a simple operator mixes only with operators [BBB] F A F B F C Γ A Γ B Γ C which preserve the set of flavor and Dirac indices in the chiral basis as (3.8) Note however that such operators are not all linearly independent. In the case of a 2N operator, we have the following constraint which comes from the general identity Here [i, j] means a simultaneous exchange between the i-th indices of F 1 , Γ 1 and the j-th indices of F 2 , Γ 2 . This identity can be generalized to from which we have where the i-th index of ABC, the j-th index of DEF and the k-th index of GHI are cyclically interchanged in (ABC, DEF, GHI)[i, j, k]. For example, Note that the cyclic interchange of indices occurs simultaneously for both Γ i and F i in the above formula. Both 2N and 3N identities are incorporated in our calculation.
The results in the tables show some notable patterns. Firstly the eigenvalues γ j /2d are all even integers; this is non-trivial since there appears to be considerable mixing in our initial operator bases. Secondly there is a tendency for the operators with larger isospin to have smaller (more negative) eigenvalues. Thirdly there are relations between the entries in the tables for different Dirac indices; e.g. the isospin degeneracies for the indices 111113344 and 111122223 in table 1 are identical and furthermore all the corresponding eigenvalues related by a common shift of −32. These observations suggest that there is an underlying algebraic structure that we have unfortunately not yet been able to reveal.
Note that a combination obtained from an other one by i) the interchange (1, 2) ↔ (3,4) or ii) the simultaneous interchange of 1 ↔ 2 and 3 ↔ 4 or iii) the interchange 1 ↔ 2 has obviously the same spectrum of anomalous dimensions and for this reason not listed separately.
The star symbol * next to the eigenvalues means that there is a corresponding operator which overlaps with the tree operator in eq. (3.5). Since the tree operator has I = 1/2, this can only happen if the corresponding n 1 is different from zero. The tree operator is invariant under the symmetry i) above whereas its two spin components are exchanged under the symmetry ii). This means that whenever the tree operator overlaps with a particular operator, it also overlaps with these transformed ones. The meaning of the symbol ♮ is an overlap between the tree operator and the transform of the operator under symmetry iii). It is intriguing that the tree operator generally overlaps with the largest eigenvalue for given Dirac indices; again a fact for which we do not yet have a simple explanation.
As a simple example, let us consider the Dirac index distribution 111112222 (third entry in table 1). The part of the tree operator relevant for this case is There are altogether 53 local 9-quark operators with this distribution of Dirac indices, but the number of independent ones is reduced to 2 after imposing all the gauge identities (3.9) and (3.12). A possible choice is (3.14) Using the gauge identities we have 15) which is proportional to the operator (of I = 1/2) corresponding to anomalous dimension -60 in our units. The other combination has anomalous dimension -84 and I = 3/2 and has no overlap with T 1 . As explained above, among the operators with Dirac index distribution 333334444 there is one with anomalous dimension -60 which also overlaps with the spin=1 component of the tree operator and among the ones corresponding to 222221111 one which overlaps with its spin=2 component. The space of operators corresponding to other Dirac index distributions are considerably larger than in this example. For example, the case 111223344 (last entry in table 4) has 1369 operators before, and 117 operators after imposing the gauge identities.
As can be seen from tables 1-4, the largest eigenvalue of γ AB is 16d (occurring already at the tree level). Therefore, the largest eigenvalue ofγ (1) AB becomes 2d × (8 − 36), which is negative, so that β tree A = −14/(33 − 2N f ). Therefore, in conclusion, the operators at the tree level in OPE dominate at short distance in the 3N NBS wave function.

Short distance repulsion of 3NF
As discussed before, the 3N potential at short distance is given by where β tree A is given in eq. (1.14). Since this result dominates over the one appearing in the 2N potential, which is of the form (1.10), the above behavior of V 2N +3N F ( r, ρ) at short distance must come solely from V 3N F ( r, ρ). Unlike for the 2NF no additional nonperturbative considerations are required in this case and therefore we can conclude that the 3NF is always repulsive at short distance, which is universal in the sense that it does not depend on the details of the 3N state, used to define the NBS wave function, such as its energy E. Note however that this conclusion is valid for the 3NF in our definition (i.e. defined from the NBS wave function) since potentials are not observable in general and therefore scheme-dependent. Unless one fixes the scheme for the definition of the potential , it is meaningless to ask whether the 3NF has a repulsive core or not. By using the potential scheme considered in this paper, we can show that the 3NF universally have repulsive cores.
Finally we would like to note that as listed in Tables 1-4, the total number of 9-quark local operators is several hundred and the spectrum of anomalous dimensions is rather dense. We singled out the ones with largest anomalous dimensions which dominate at short distance but it is not evident at what distance scale this leading behavior sets in. Even if it turns out that these individual operators really dominate at extremely short distances only, our main conclusion may remain valid due to the fact that all the eigenvalues of the effective gamma matrix are negative (corresponding to short distance repulsion). We think a simple explanation of this fact should exist (maybe related to the Pauli principle).
An interesting and important extension of the present analysis is to investigate the short distance behavior of the three baryon force (3BF) by the same method. Its results can tell us whether there is a universal short distance repulsion also in the 3BF, which has been suggested to explain the observed maximum mass of neutron stars [7].