Deuteron Matrix Elements in Chiral Effective Theory at Leading Order

We consider matrix elements of two-nucleon operators that arise in chiral effective theories of the two-nucleon system. Generically, the short-distance piece of these operators scales as 1/r^n, with r the relative separation of the two nucleons. We show that, when evaluated between the leading-order wave functions obtained in this effective theory, these two-nucleon operators are independent of the cutoff used to renormalize the two-body problem for n=1 and 2. However, for n greater than or equal to 3 general arguments about the short-distance behavior of the leading-order deuteron wave function show that the matrix element will diverge.

Introduction: In recent years chiral effective theory (χET) has achieved significant prominence as a technique via which model-independent results can be obtained for few-nucleon systems (see Refs. [1,2] for reviews). χET is based on the realization that the interaction between pions and nucleons is governed by the approximate chiral symmetry of QCD. Therefore the use of heavybaryon chiral perturbation theory, together with Weinberg's proposal that in few-nucleon systems the quantity that has a well-behaved chiral expansion is the nucleonnucleon potential V [3,4], seems to facilitate a systematic calculation with which nucleon-nucleon scattering can be well described-as long as the collision energies are significantly below the chiral-symmetry breaking scale Λ χ [5,6,7]. In this approach V -and also by extension two-nucleon-irreducible operators for other processes-is expanded as a chiral series in the usual chiral perturbation theory expansion parameter where m π is the pion mass, p the momenta of the nucleons involved in the scattering, and M the nucleon mass. Although questions have been raised about the consistency of such a power counting [8,9,10,11] this approach has had considerable success in describing the scattering data in the two-and three-nucleon sector. The theory has also been shown to be consistentin the renormalization sense-in the 3 S 1 -3 D 1 channel, the channel where the two-nucleon bound state deuterium occurs [9,12]. The leading-order wave function of deuterium, |ψ LO can therefore be obtained by solving the Schrödinger equation-in either momentum-space or configuration-space-for two nucleons interacting via the piece of the N N potential which is of chiral order zero. As first realized by Weinberg [3,4] that piece consists of the venerable one-pion-exchange potential, together with a (momentum-independent) four-nucleon contact interaction. This potential is singular: the Hamiltonian it generates is unbounded from below. It therefore requires regularization and renormalization. In practice, the potential is regulated at some momentum scale Λ (or in coordinate space at a distance 1/Λ). The strength of the contact interaction is then adjusted to reproduce some observable, usually the deuteron binding energy. If this can be done over a wide range of Λ, and if other N N scattering observables are independent of Λ up to corrections of higher order in the χET, then we say that the potential has been renormalized. In Ref. [9] Beane and collaborators showed that this could be done for the 3 S 1 -3 D 1 channel-at least if one was only concerned about the results in the chiral limit m π = 0. This conclusion has since also been reached using a different regulator and in a regulator-independent fashion in Refs. [11,12,13].
With these wave functions in hand it is natural to reexamine the many successes that phenomenological potential models have had in describing the deuteron's interaction with external probes, such as electrons, pions, or photons. In such a calculation we use an operatorÔ that is appropriate to the particular external probe under consideration, and has been derived using a chiral expansion in powers of P . That operator is then sandwiched between a wave function obtained from a chiral potential (which presumably should be computed to the same relative order), yielding matrix elements: Many such calculations have been performed (see Ref. [14] for a recent, partial summary), although, with one exception which we will discuss below, none of them have employed the wave functions of Refs. [9,11,12]. In these calculations the bound-state wave function is found by employing some regularization prescription involving a cutoff Λ. Here we will discuss the conditions under which the use of such wave functions in the evaluation of Eq. (2) yields a matrix element that is independent of Λ. The extent of the Λ dependence in the result for M tells us the degree to which the χET prediction for the matrix element is model independent.
Recently, Meißner et al. [15] and Nogga and Hanhart [16] considered this question for pion-deuteron scattering, a reaction originally discussed within this framework by Weinberg in 1992 [17]. There a two-nucleon operator representing the process depicted in Fig. 1(a) gives a large contribution to a πd . The presence of this large double-scattering term can obscure attempts to extract the isoscalar pion-nucleon scattering length from a πd . Pion scattering on the individual nucleons in the deuterium nucleus (see Fig. 1(b)) is O(P 2 ) in the χET, and yields a piece of the πd scattering length [17,18]: with µ = m π /M . Meanwhile the double-scattering diagram of Fig. 1(a) gives [17,18,19]: where here, and throughout what follows, with u and w the 3 S 1 and 3 D 1 radial deuteron wave functions.
The contribution of Eq. (4) to the pion-deuteron scattering length is suppressed by one power of P relative to the nominally-leading contribution (3). However, since corrections to the leading-order N N potential are suppressed by two powers of P [5], a next-to-leading order computation of a πd involves the evaluation of the matrix element: where the subscript LO indicates that the matrix element should be taken with deuteron wave functions obtained with the leading-order χPT N N potential, and the subscript Λ is included in Eq. (6) because a cutoff must be applied to that potential before it is renormalized to give the correct deuteron binding energy. Since there are no short-distance pieces of the π-d operatorÔ at next-to-leading order the matrix element (6) must be cutoff independent over a significant Λ range if Weinberg's approach is to be a consistent way to calculate pion-deuteron scattering. Nogga and Hanhart gave numerical evidence that the matrix element (6) has a Λ → ∞ limit that exists and is finite. (This is in accord with a similar conclusion previously obtained in Ref. [12], where 1 r was evaluated using wave functions derived solely from one-pion exchange. The details of the evaluation of Ref. [12] will be explained further below.) This shows, Nogga and Hanhart claim, that expandingÔ in powers of P is a consistent way to calculate the pion-deuteron scattering length. Consequently, at least in principle, it should be possible to perform a high-accuracy, model-independent extraction of the pion-neutron scattering length from pion-deuteron scattering data.
But other two-nucleon operators that occur in πd scattering computations have higher powers of r in the denominator in their co-ordinate space form. For instance, the next term in the πN multiple-scattering series, depicted in Fig. 1(c), yields a contribution to a πd [20]: This contribution is of order P 5 , and so is N 3 LO in the chiral expansion for a πd , but numerically it is the next two-nucleon effect that must be considered after the inclusion of the matrix element (6) [20]. The failure of the standard chiral expansion to account for the large size of such effects is related to the absence of any suppression of the operator (7) at long distances. It led the authors of Ref. [20] to propose a different power counting for the πd scattering length-one that more accurately captures the relative hierarchy of mechanisms contributing to a πd . Regardless of what the correct counting for the operator (7) is, if the impact of the triple-scattering mechanism of Fig. 1(c) on extractions of the isoscalar pion-nucleon scattering length from a πd is to be assessed we must not only consider the matrix element (6), but also: More generally, all two-nucleon operators in the chiral effective theory will scale as 1/r n at distances ≪ 1/m π and ≪ 1/q (with q the magnitude of any momentum or energy transferred to the nucleus by the probe). The value of n is determined by the process under consideration and the order to which the two-body operator is calculated, with higher values of n being reached as higher orders are computed in the chiral expansion for O. In a perturbative calculation of deuteron matrix elements within χET the highest-order (and most singular) pieces of the operator will contribute to the overall result 3 for M through their expectation value taken with the leading-order wave function. Thus, we are led to a very general question: If one considers the expression then for which values of n does the Λ → ∞ limit exist, and for which is it finite? It is this question, together with associated ones involving operators that connect the S-and D-wave components of the deuteron wave function, that we will answer in this paper.
Theory: Nogga and Hanhart solved the Schrödinger equation in its momentum-space form, i.e. the homogeneous Lippmann-Schwinger equation: where Λ is the scale at which the potential V is regulated, and G 0 (p) = (−B d −p 2 /M ) −1 is the (free, center-of-mass frame) two-nucleon propagator, with B d and M denoting the deuteron binding energy and nucleon mass, respectively. The leading-order potential is given by a one-pion exchange (OPE) contribution and a short-distance piece: with q ≡ p ′ − p the three-momentum of the exchanged pion. In Eq. (11) P t is a projection operator that projects onto the 3 S 1 channel, and C t the strength of the shortdistance potential in that channel. While the OPE contribution is totally determined at leading order through the pion mass m π , the axial coupling constant g A = 1.26, and the pion-decay constant f π = 92.4 MeV, the contact interaction parameter C t has to be determined from N N data and will be a function of the cutoff Λ. From Eqs. (10) and (11) it is straightforward to obtain the coupled one-dimensional differential equations which describe the deuteron wave function: where u and w are, as defined above, the deuteron radial wave functions. The coupled-channel potential is given by with Equations (12)- (14) will be equivalent to Eqs. (10) and (11) provided that r > 1/Λ. The equations (12) are solved subject to the following boundary conditions as r → ∞: with γ = √ M B d the deuteron wave number, A S the normalization constant which guarantees that dr u 2 (r) + w 2 (r) = 1 , and η the asymptotic D/S ratio. By employing Eqs. (12)- (14), and solving them for arbitrarily short distances r, Pavon Valderrama and Ruiz Arriola have calculated u and w using boundary conditions at inter-nucleon distances of order 0.1-0.2 fm [12]. This echoes the much earlier work of Sprung and collaborators [21], as well as the solution of the leading-order potential with a square-well regulating the short-distance behavior [9]. Within such an approach, an analysis of the asymptotic short-distance behavior of the components u and w gives the following result [9,12,21]: C 2A and C 2R are normalization constants which have been determined in [12], R is a new length scale that enters the non-perturbative problem. It is defined by . When Eqs. (12)- (14) are solved in Ref. [12] the phase φ is determined by the boundary condition at r = 0, and so φ is regulator independent, and is a function only of the scales m π , γ, and R. After computing the numerical solution of Eqs. (12) to a sufficiently small radius r, it can be matched to the r → 0 form of the deuteron wave function (17) and an-in principle-regulator-independent wave function can be obtained 1 .
Results: Let us now return to the problem of computing First, we split the integral up as: where R * is sufficiently small that the asymptotic forms (17) apply, but is still ≫ 1/Λ. The piece of the integral from R * to ∞ can be calculated within χET and will depend only on low-energy scales such as γ and m π , and, of course, on R * itself. Meanwhile, as Λ → ∞ the first integral goes to zero, as long as the integrand is integrable. Note that any regulator dependence of the wave functions should be contained in that first piece of our integral, since we have already explained that the wave functions for r > 1/Λ are regulator-independent solutions of Eqs. (12)- (14). Therefore, in order to establish whether or not 1/r n Λ is regulator-independent the pertinent piece of Eq. (19) is If this has a Λ → ∞ limit that is finite then the entire matrix element 1/r n Λ will also be well behaved in that limit. Inside the integral Eq. (9) we may substitute the expressions Eq.(17) for u and w. We then see that shortdistance integrals involving the exponential piece of the wave function will always converge for any n, as the exponential itself regularizes the result. On the other hand, the term including the cosine function does not vanish at the origin and therefore we have to determine the values of n for which the integral 1 Although the final wave function is in principle regulator independent, the boundary conditions employed in the numerical computation of the long-range components can lead to cutoff effects at short distances and this leads to a small numerical uncertainty of the constants C 2A , C 2R quoted in [12]. However, with sufficient care this uncertainty can be made arbitrarily small.
gives a finite result as Λ → ∞. This question is easily answered by simple dimensional analysis. For n = 3 or higher the integral in Eq. (21) will diverge. However, for n = 1, 2 if we ignore the exponential pieces of u and w we can solve the integrals analytically and obtain where: The function f n can be written as a linear combination of incomplete Gamma functions. Its asymptotic expansion is: f n (x) → α n +h n (a)+x 2n−6 sin(x+2φ)+O(x 2n−7 ), (24) as x → ∞ (n = 1, 2), with: The function h n (a) can easily be evaluated, but its value is not relevant for our purposes here, since it cancels in Eq. (22). Therefore for n = 1 and 2 the matrix element has a finite Λ → ∞ limit. However, that limit is approached in the presence of cutoff-dependent oscillations whose amplitude scales as (RΛ) n−3 . Another cycle in this oscillation occurs whenever 1/Λ becomes small enough that we integrate over another node in the deuteron wave function. The presence of oscillations is thus associated with the appearance of spurious bound states in the effective theory-bound states whose binding energies are greater than the theory's breakdown scale. However, for the n = 1 and n = 2 matrix elements, the oscillations vanish as Λ → ∞, so even though such bound states are present, the procedure of evaluating the matrix element with a wave function at a given Λ should yield the correct (i.e. Λ → ∞) answer, as long as a sufficiently high cutoff is chosen. Now adding the Λ → ∞ piece of these matrix elements to the piece that comes from integration between R * and infinity, we find: Cutoff-dependence of the matrix element 1 r Λ, using the radial deuteron wave functions employed in [16]. The cross is the result found using the wave function of Ref. [12], as given in Eq. (26).
where the first result is in agreement with the number given in Ref. [12]. As pointed out there, it is essentially consistent with the range 0.450-0.465 obtained for 1/r using a variety of potential-model wave functions in Ref. [20]. The result (27) is new, and does not fall within the range 0.286-0.345 fm 2 quoted in Ref. [20]. It takes longer to reach asymptopia for this matrix element.
In order to check the results of Eqs. (26) and (27) we have computed the expectation value ψ LO | 1 r n |ψ LO for n = 1 and 2. We used the same wave functions as were used in [16] with cutoff values between 2 and 20 fm −1 . As Fig. 2 shows, we find similar behavior for 1/r as was displayed in that paper, and a limiting value that is consistent with that found by Nogga and Hanhart [22]. But, we can now interpret the Λ-dependent oscillations seen in Fig. 2 as exactly the ones predicted by Eq. (22). Meanwhile, Fig. 3 shows that, as expected based on our analysis of Eq. (22), these oscillations are more pronounced for 1/r 2 , and in consequence a much higher cutoff is needed to achieve a converged result. However, also for this matrix element the numerical results indicate that a finite limiting value exists, which is in accord with our previously presented analytic arguments.
Since Eq. (17) identifies the relevant behavior of the r → 0 piece of the deuteron wave function, we are able to predict that all powers n ≥ 3 will lead to divergent results. This result is supported by numerical calculations of the corresponding matrix elements with the wave functions used in [16].
We also find that, to a very good approximation, the result for agrees (for n = 1 and 2) with the numerical results obtained with the co-ordinate space wave function of Cutoff-dependence of 1 r 2 Λ, using the radial deuteron wave functions employed in [16]. The cross is the result found using the wave function of Ref. [12], as given in Eq. (27).
Ref. [12] that we presented in Eqs. (26) and (27). (The slight disagreement exhibited in Figs. 2 and 3 between the trend of the dots and the cross that represents the result of Eqs. (26) and (27) can presumably be traced to the 3% difference in the πNN coupling constants used in Refs. [12] and [16].) This is an important result, because, together with the even better agreement for matrix elements of positive powers of r (see Table I of Refs. [12,16]) it suggests that the procedure of constructing a limiting sequence of cutoff-dependent wave functions and evaluating the matrix element as a function of Λ is not actually necessary. Instead such matrix elements can be computed using the co-ordinate space results for the deuteron wave function given in Refs. [9,12,21].
Next we consider integrals of the form: For n = 1 this integral occurs in other, numerically less important, corrections to the πd scattering length. For n = 2 contributions of the form (29) occur in the evaluation of sub-leading pieces of two-nucleon effects in the process γd→ π 0 d. Since u and w both display the same short-distance behavior, we arrive at the same conclusion as for the matrix elements considered above: for n = 1 and n = 2 the result is convergent (in spite of oscillations), see Figs. 4 and 5, while for n ≥ 3 it is divergent. For n = 1 and 2 evaluation with the wave functions of Ref. [12] gives:  until Λ ≫ 20 fm −1 , since approximately 20% of the final number accrues in the region between r = 0 and r = 0.1 fm.
Conclusion: In this paper we have shown explicitly that deuteron matrix elements of two-nucleon operators which are proportional to 1/r n at short distances converge for n ≤ 2 and diverge for n ≥ 3 when they are evaluated using the leading-order deuteron wave function. We have given numerical evidence by explicit calculations using a sequence of leading-order deuteron wave functions corresponding to different ultraviolet cutoffs, and the results obtained in this way agree with an analysis based on the short-distance behavior of the deuteron wave function.
As higher orders are calculated in the chiral series for the operatorÔ the divergent case n = 3 will be reached (e.g. γd→ π 0 d at O(P 4 ) in χPT [23]). Suppose this occurs at some chiral order m. If both the wave function |ψ and the operatorÔ are written as a chiral series, then this implies that an mth-order piece of the matrix element (2), specifically the piece ψ LO |Ô (m) |ψ LO , is divergent. This apparently mandates the presence of a contact interaction involving both nucleons and the external probe inÔ (m) , so that this divergence can be absorbed. However, in explicit calculations of a particular process to order m in the chiral expansion this divergence may cancel with contributions to (2) due to corrections to the leading-order wave function. Such corrections will come from higher-order pieces in the chiral expansion of the N N potential, and will generate contributions to (2), e. g., of the type δ ψ (m) |Ô LO |ψ LO . Therefore, the appearance of a matrix element of an operator 1/r n with n ≥ 3 at order m does not immediately indicate that a counterterm is needed at that order. However, if a counterterm is not permitted at that order by chiral symmetry and/or electromagnetic gauge invariance, then the various divergences at order m must cancel each other. This represents a constraint on the sum of all mechanisms that contribute to (2) at that order. (For an example of such a cancellation in the context of pion production, see Refs. [24,25].) Lastly, we note that if an N N potential that is more singular than one-pion exchange is iterated to all orders using the Schrödinger equation, as is done in Refs. [5,6,7], then the behavior of the wave functions u and w will not follow the form (17). Indeed, the more singular the potential, the more convergent u and w will become [26]. Therefore, the conclusion of the previous paragraph is limited to an approach to higherorder calculations in the chiral effective theory where all corrections to the leading-order |ψ are evaluated in perturbation theory, and so the form (17) still represents the dominant short-distance behavior of u and w.