Deuteron disintegration in three dimensions

We compare results from the traditional partial wave treatment of deuteron electro-disintegration with a new approach that uses three dimensional formalism. The new framework for the two-nucleon (2N) system using a complete set of isospin - spin states made it possible to construct simple implementations that employ a very general operator form of the current operator and 2N states.


Introduction
The theoretical description of electromagnetic processes is based, like most nuclear physics formalisms, on a partial wave decomposition of relevant operators. This restricts theoretical methods to systems where a relatively small number of partial waves is important. Recently, several formalisms for the three dimensional description of few-body systems and processes have been developed. In this paper we fill a gap in the development of the three dimensional framework and present an approach that allows for a simple implemen-tation of three dimensional electromagnetic currents. Our final expressions can be translated to a numerical implementation via direct substitutions of 16 dimensional square matrices representing operators in the 2N isospinspin space. Finding the matrix representation of relevant operators is greatly simplified by using symbolic programming in the Mathematica R [1] software package. Our approach allows us to use a very general operator form of current operators and can therefore be used for a wide class of processes.
In this paper we apply it to the case of electron induced deuteron disintegration and compare the results with traditional partial wave calculations. It is worth noting that the methods presented in the following sections can be applied to the description of other processes, where electroweak probes interact with the 2N system. Electromagnetic currents can be replaced by any operators acting on the same degrees of freedom; this makes our implementations also useful for calculations involving for example muon capture or neutrino induced deuteron disintegration, performed recently still with the use partial wave expansion [2,3].

Formalism and notation
We adopt a notation in which capital letters describe the total momentum of a two particle system (P = p 1 + p 2 ), lower-case letters describe the relative momentum (p = 1 2 (p 1 − p 2 )). Subscripts denote individual particles and superscripts assign a momentum to a particular quantum eigenstate. The two particle momentum eigenstates are normalized such that: and the transition from the total and relative to the individual momenta can be achieved using: where in (3) and in the following we assume that the difference between the proton and neutron mass is negligible. We examine the case of deuteron disintegration (e + 2 H → e + p + n) where the 2N system is treated in the non-relativistic approximation. In the initial state the deuteron is at rest (P = 0) and the electron has a momentum magnitude of q e . We assume that the rest mass of the electron is negligible in comparison to its kinetic energy therefore the initial electron energy E e ≈ q e . The final electron momentum magnitude is q ′ e , the final energy E ′ e ≈ q ′ e and the electron scattering angle is θ e The magnitude of the three momentum transferred to the 2N system is Q = q 2 e + q ′2 e − 2q e q ′ e cos θ e ≈ E 2 e + E ′2 e − 2E e E ′ e cos θ e .
We work in a reference frame, where the momentum transfer is parallel toẑ. In this frame momentum conservation leads to the expression for the total momentum of the proton and the neutron in the final state: The magnitude of the final relative momentum can be calculated from the energy conservation: where E d is the (negative) deuteron binding energy and the direction of p f can be arbitrary. The crucial nuclear matrix element M µ between the initial deuteron state (where the total momentum P = 0 and the two particle total angular momentum has aẑ projection m d ) and the final 2N scattering state can be expressed in terms of the full 2N current operator (j µ 2N ) and the t operator: In this equation the final state is anti-symmetrized: In the first term on the right hand side of (8) particle 1 has a spin (isospin) z projection m 1 (ν 1 ), particle 2 has a spin (isospin)ẑ projection m 2 (ν 2 ) and in the second term on the right hand side of (8) the particles are exchanged. The j(1) (j(2)) operator is a single nucleon current acting on the degrees of freedom of the first (second) nucleon. The j(1, 2) operator accounts for processes where two nucleons are involved, t is the 2N transition operator and G 0 is the free 2N propagator. The energy argument of the transition operator and the propagator is E = (p f ) 2 /m. Finally, the µ index denotes the component of the current operator. In particular µ = 0 stands for charge density operators while µ = 1, 2, 3 stand for spatial components. In the following sections the way we calculate the individual parts of M µ in (7) will be discussed separately for a specific choice of the coordinate system and value of µ; we will drop µ for brevity. 4

Deuteron bound state
The structure of the deuteron wave function can be written in the operator form, following [4,5,6]: where In equation (9) | χ(m d ) is a state in which the spins of the two spin 1 2 particles are coupled to a total spin 1 with aẑ projection m d . The isospins of the two nucleons are coupled to the total isospin 0 state | 0 0 . φ l are scalar functions of the relative momentum and σ(1), σ(2) are doubled spin operators in the spin space of one nucleon and identity operators in the spin space of the other nucleon, respectively: Vector components in (10) and (11) are operators written in terms of the tensor product ⊗. They act in the 4 dimensional spin space of the two nucleon system and can be represented by 4 × 4 matrices -tensor products of identity operators and Pauli matrices. Note that in equation (9) states and operators in the 2N isospin -spin space are placed inside the [. . . ] brackets; in Sec. 7 we will introduce a simple way to implement these expressions as 16 dimensional vectors and 16 × 16 matrices for use in our numerical treatment.
Scalar functions φ l in expansion (9) can be calculated using three dimensional formalism, see for example [6]. Nowadays deuteron bound state calculations can use any 2N potentials given in the operator form and do not require substantial computational resources.

Single nucleon currents in three dimensions.
Single nucleon (1N) currents act on the degrees of freedom of one particle. Their matrix elements in the momentum space depend only on the initial and final momenta and are operators in the isospin -spin space. For example the matrix element for the second nucleon, j(2), reads: where in view of the standard nonrelativistic current, the dependence on the difference and sum of the initial and final momenta is used. Implementing the transition from the individual particle momenta to the relative momenta leads to Again, the expressions inside the square brackets [. . . ] can be easily represented using the notion of the Kronecker product, see Sec. 7. The action of j(2) on the deuteron state at rest can be worked out: where we used P = 0 and the normalization of momentum eigenstates (1)- (2) so O 1N is the resulting single particle operator. Equation (14) gives the full isospin -spin state for the final p ′ , P ′ momenta.

2N currents in three dimensions.
For a wide class of 2N current operators, their matrix elements in the momentum space (operators in isospin -spin space) are given in the form: see for example [7,8,9]. The right hand side of (15) is as a linear combination of scalar functions (f jS i , f j i ) and products of spin space operators (O jS , O j ) and isospin space operators (T i ): where the subscript S distinguishes density and current components. In Ref.
[10] a general operator basis for the local 2N current operator was introduced. In Sec. 7 some examples of operators from (16), (17) will be used to demonstrate our matrix representation of expressions inside [. . . ]. 6 Again, using (3), current matrix elements become: In the following we restrict ourselves to this class of momentum dependences. Our approach can, however, be generalized to include any type of momentum dependence. The action of j(1, 2) on the deuteron state can be worked out using (1),(2), (3) and P = 0. In the laboratory frame it yields: where O 2N is the resulting two-particle operator. Equation (19) gives the full isospin -spin state for the final p ′ , P ′ momenta.

t operator in three dimensions
The t operator satisfies the Lippmann -Schwinger equation: where G 0 (E) is the free propagator depending on the energy E and V is a 2N potential. It follows that, as shown in [11] also t can be written as a linear combination of scalar functions (t γ,i ) and operators (W γ,i ) in the isospinspin space: Here are again operators in the isospin -spin space (matrix elements between momentum states), with w i (i = 1, 2, ..., 6) acting in the 4 dimensional spin space of the 2N system. The decomposition (22) is not unique; our choice of the six w i operators is consistent with [6] . Scalar functions arising in the decomposition of t can be calculated in the three dimensional formalism. 7 Calculations can be performed for any type of the NN potential satisfying a similar expansion (22). For details see Ref. [6].
The C γ =| γ γ | isospin operators project onto one of the four 2N isospin states: The | γ states are chosen in this way, because t(E) conserves the total 2N isospin. The rescattering part of the matrix element M in (7) can be written as: where O is either O 1N from Eq. (14) or O 2N from Eq. (19), E = p ′2 /m is the relative energy of the final 2N state and The integral (30) with the cut-off valuep can be easily calculated numerically. In the next section we show how to prepare its component [O(p, P ′ )].

Representation of spin-isospin operators
Now that the form of expressions in (7) has been established it remains to find a way to represent operators and states inside [. . . ]. Once an appropriate matrix representation is found, numerical calculations can be constructed 8 Table 1 Reference quantum numbers for our KP states.
using equations (14), (19) and (30) where m spin(isospin) j are the spin (isospin) projections of nucleon j and the corresponding quantum numbers are given in Tab. 1 Any operator or state in this space can be constructed using the notion of the Kronecker product (KP). The Mathematica R [1] symbolic programming software contains definitions for the KP, which makes translating expressions for isospin -spin operators a simple task. However, it is important to remember that the order of operators in the KP must be preserved. Tab. 1 can serve as a reference to keep consistence with this paper.
The deuteron operators [B 1 (p)] and [B 2 (p)] with p = (p x , p y , p z ) from equation (9) where B is a 4 × 4 matrix: Operators appearing in (16) and (17) have a more complicated isospin operator form. A simple example of their matrix representation is given below in Fig. 1.
States can be constructed in a similar manner using the built in Mathematica R definitions for Clebsch-Gordan coefficients, for example [| 0 0 ⊗ χ(1)] results in:

Results
In the following we will present results obtained using a chiral NNLO potential [12] with Λ = 550 MeV/c andΛ = 600 MeV/c. The operator form of such a potential was briefly described in Appendix C of Ref. [6], where also an example set of necessary parameters was given for its neutron-proton version. The same parameters will also be used in the present paper.
There are three basic ingredients in our calculations: the deuteron wave function, the 2N t-matrix and the 2N current operator. Before we show selected observables for the 2 H(e, e ′ p)n reaction, we will describe our numerical performance and the way we verify the quality of our calculations.
As described in Ref. [6] and equation (9), the deuteron in the operator form is represented by two functions φ 1 (p) and φ 2 (p). The corresponding Schrödinger equation for φ 1 (p) and φ 2 (p) can be rewritten as an eigenvalue problem, which is of the same type and dimension as the one solved for the deuteron wave function in the standard partial wave representation, where one deals with the s-and d-components, ψ 0 (p) and ψ 2 (p). The connection between the solutions is very simple [5] and can be used to check the numerical performance.
In Fig. 2 we show directly the φ 1 (p) and φ 2 (p) functions required for the operator expansion of the deuteron. In Fig. 3 the resulting s-and d-wave components in momentum space are compared to the results obtained by firstly decomposing the NN potential into partial waves and then solving the resulting Schrödinger equation in its standard form. The agreement for the two wave function components is perfect for all their significant values. In Ref. [6] we solved the Lippmann-Schwinger equation (LSE) for the 2N t-matrix directly in three dimensions. At that time we focused mainly on the on-shell behavior of the expansion coefficients t i (p ′ , p, x; E 2N ), that is we were interested in t i p 0 , p 0 , x; E 2N = p 2 0 m , which are sufficient to calculate the Wolfenstein parameters and the nucleon-nucleon scattering observables. Furthermore, we solved the LSE in such a form (Eq. (2.6) from Ref. [6]) that the magnitude of the initial p momentum could be fixed. Clearly, for the 2 H(e, e ′ p)n reaction we need a "left" version of Eq. (36), which allows us to find the half-shell t-matrix for a fixed final relative momentum, p f , given now by the reaction kinematics. The starting point for this new version is equation (20). Repeating the algebra outlined in Ref. [6], we prepared a numerical realization of this "left" version of the LSE, leading to the scalar expansion coefficients t i p f , p; where x ≡p f ·p. Our numerical scheme was again based on matrix inversion and used the standard LU decomposition of Numerical Recipes [13]. In order to achieve a unique and smooth solution also for p = p f , it was sufficient to calculate the average with δ p f ≈ 0.01 fm −1 . Actually, this effort turned out to be unnecessary and provided merely an additional check of numerics, since for the most general rotational, parity and time reversal invariant form of the NN force. That means that the left coefficients, t i p f , p, x; E 2N = (p f ) 2 m , can be obtained directly from the "right" version of LSE.
In order to further check our t-matrix coefficients, we used them to calculate the partial wave representation of the t-matrix: where l (l ′ ) is the initial (final) angular momentum of the 2N system, s is the 2N (conserved) spin and j is the 2N (conserved) total angular momentum. These matrix elements can be compared with direct solutions of LSE obtained in the standard partial wave representation. We performed the projection of the three-dimensional t-matrix on partial waves, employing the simple method proposed in Ref. [14] for NN forces. In Figs. 4-7 we show examples for the uncoupled and coupled channels, with the 2N isospin t = 0 and t = 1. We chose p f ≈ 1.9 fm −1 , which corresponds to a relatively high NN center of mass energy, E 2N = 150 MeV. For such an energy many partial waves contribute to the NN scattering observables and the question arises if the partial contributions are consistent with the full three-dimensional calculations. From Figs. 4-7 (and many other cases which are not shown here) we infer that this is really the case. The agreement between results based on the two quite different approaches is very good. The final ingredient in our framework is the 2N current operator. It consists of the single-nucleon and two-nucleon operators. For the purpose of this paper we assume that its single-nucleon part comprises the standard nonrelativistic charge density as well as the convection and spin current operators. In the 2N part we take for simplicity only the leading one-pion-exchange current operator in the chiral effective field theory representation. In our threedimensional treatment of the 2 H(e, e ′ p)n reaction, we calculate the spin and isospin matrix elements of the current operator directly, using simple matrix representations of the spin and isospin operators and the concept of the Kronecker product to deal with the two-nucleon spin and isospin spaces.
In the traditional calculations, a partial wave decomposition of the current operator is required. It is a rather easy task for the single-nucleon part of the current. For the one-pion-exchange current operator it is known analytically (see for instance Ref. [15]). It can also be obtained using the method proposed in Refs. [16,17], where even more complicated two-pion-exchange current operators were considered.
To give examples of our results on deuteron electro-disintegration we chose several electron kinematics given in Table 2. They allowed us to study the reaction for three different internal nucleon-nucleon energies (corresponding to the three values of p f ) and for five values of the three-momentum transfer Q. The first parameter is the input for the t-matrix calculations and the second one specifies the properties of the current matrix elements. Table 2 The six electron kinematics considered in the paper for the exclusive 2 H(e, e ′ p)n process. The initial electron energy (Ee), the electron scattering angle (θe), the final electron energy (E ′ e ), the final relative nucleon-nucleon momentum (|p f |), the energy transfer (ω), and the magnitude of the three-momentum transfer (Q) are given. In addition to the information given in Table 2, we need to label the exclusive kinematics. For the fixed "electron arm", we deal in fact with a two-body kinematics in the final proton-neutron system. We restrict ourselves here to the case where the proton is ejected in the electron plane, where θ p would be the angle between the three-momentum transfer Q and the final proton momentum p p . Since we have to distinguish between the φ p =0 deg and φ p =180 deg cases, we ascribe the negative sign to θ p for φ p =0 deg. This is shown in Fig. 8. Note that the six electron kinematics provide a unique solution for any θ p value and that θ p changes from 0 do 180 degrees.
We are now ready to show our results for several selected observables. We chose first of all the unpolarized cross section, d 5 σ/(dE e ′ dΩ e ′ dΩ p ). We take also into account one example of the spin-dependent helicity asymmetry, where h is the initial electron helicity and the projection of the initial deuteron total angular momentum (J d ) on Q, J dz , is equal 1. In addition we show our predictions for the deuteron tensor analyzing powers T kq . Note that they are calculated in the system, where Q ẑ. Our primary goal was to compare results based on the partial wave decomposition for the t-matrix and the nuclear current operator with new pre-dictions resulting from the three-dimensional scheme. We observed a perfect agreement for all the electron kinematics and for all the considered observables, if a sufficient number of partial waves in the first type of calculations is included. The six kinematics can be divided into two groups: (K1 , K3 , K5) and (K2 , K4 , K6). In each group a similar type of convergence of the observables with respect to the number of partial waves is observed. That is why in Figs. 9-14 we show predictions for two representative (K3 and K6) kinematics only. In the first case we see a rapid convergence and partial wave based results with j ≤ 4 are already very close to the full three-dimensional prediction. In the second case all partial waves with j ≤ 9 are necessary to achieve convergence.
It is interesting to see that slow convergence for the K2, K4 and K6 kinematics does not result from the higher p f values (that is from the t-matrix) but is related to the Q values and thus to the partial wave decomposition of the current operator. It is well known (see for example Ref. [16]) that especially the partial wave decomposition of the single nucleon current requires many partial waves. However, even if the initial bound state is given in the partial wave representation, the single nucleon current can be applied directly in the case of the plane wave amplitudes. This holds not only for the two-but also for the three-nucleon system [18]. In order to demonstrate this behaviour, we showed in Fig. 15 observables for the K6 kinematics. In this case the single-nucleon current contribution to the plane wave amplitude is calculated without partial wave decomposition. We clearly see that the convergence is significantly improved, even if the two-nucleon current contribution to the plane wave part of the nuclear matrix element and the whole rescattering part of the nuclear matrix element is calculated with the partial wave decomposition. Re ( Fig. 8 The kinematics for the exclusive 2 H(e, e ′ p)n process. k (k ′ ) is the initial (final) electron momentum. We neglect the electron mass, so | k |= E and | k ′ |= E ′ .  Fig. 9 The unpolarized cross section d 5 σ/(dE e ′ dΩ e ′ dΩp ) as a function of the proton scattering angle θp for the K3 electron kinematics from  Fig. 10 The same as in Fig. 9 for the spin-dependent helicity asymmetry A .  Fig. 11 The same as in Fig. 9 for the deuteron analyzing powers T kq .  Fig. 12 The same as in Fig. 9 for the K6 electron kinematics from Table 2. In the right panel partial wave based results with j ≤ 4 (dash-dotted line), j ≤ 7 (dotted line) and j ≤ 9 (dashed line) are compared with the full three-dimensional prediction (solid line).   Table 2. The single-nucleon current contribution to the plane wave part of the nuclear matrix element is now calculated without partial wave decomposition. The remaining parts of the nuclear matrix element are still calculated using partial waves. Results obtained with j ≤ 1 (dash-dotted line), j ≤ 2 (dotted line) and j ≤ 3 (dashed line) are compared with the full three-dimensional prediction (solid line).

Conclusions and outlook
The presented method to treat several electroweak processes involving 1N and 2N current operators in three dimensions can successfully replace standard partial wave treatment. We showed, for the case of electron induced deuteron disintegration, that results obtained using the new approach agree very well with those obtained using PWD. For all observables considered in this paper, the traditional results converge to the three dimensional predictions when the number of partial waves is sufficiently high (Figs. 9-14). Our formalism employs the two nucleon bound state, the 2N t matrix and the current operators in the joined isospin -spin space of the 2N system using the three dimensional formalism in 2N momentum space. Each element of this framework has been separately tested and compared with the standard PWD approach. Our method seems to be more flexible and can deal with the rich structures of the 2N force and the current operator, especially derived within the higher orders of the chiral effective field theory [12,19,10]. We plan to use our framework for other processes such as muon capture or neutrino induced deuteron disintegration. In [2] muon capture on 2 H and 3 He was considered using the PWD approach. It would be interesting to compare those results with three dimensional calculations. A similar convergence to three dimensional results as in Figs. 9-14 is expected.