Relativistic effects in two-particle emission for electron and neutrino reactions

Relativistic effects in two-particle emission for electron and neutrino reactions I. Ruiz Simo,1C. Albertus,1J. E. Amaro,1M. B. Barbaro,2J. A. Caballero,3and T. W. Donnelly4 1Departamento de Física Atómica, Molecular y Nuclear, and Instituto de Física Teórica y Computacional Carlos I, Universidad de Granada, Granada 18071, Spain 2Dipartimento di Fisica, Università di Torino and INFN, Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy 3Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, Apartado 1065, 41080 Sevilla, Spain 4Center for Theoretical Physics, Laboratory for Nuclear Science and Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA (Received 17 May 2014; published 18 August 2014) Two-particle two-hole contributions to electroweak response functions are computed in a fully relativistic Fermi gas, assuming that the electroweak current matrix elements are independent of the kinematics. We analyze the genuine kinematical and relativistic effects before including a realistic mesonexchange current operator. This allows one to study the mathematical properties of the nontrivial sevendimensional integrals appearing in the calculation and to design an optimal numerical procedure to reduce the computation time. This is required for practical applications to charged-current neutrino scattering experiments, in which an additional integral over the neutrino flux is performed. Finally, we examine the viability of this model to compute the electroweak two-particle–two-hole response functions. DOI: 10.1103/PhysRevD.90.033012 PACS numbers: 13.15.+g, 13.60.Hb, 25.30.Fj, 25.30.Pt I. INTRODUCTION The understanding of intermediate-energy (0.5–10 GeV) neutrino-nucleus scattering cross sections is an important ingredient to atmospheric and accelerator-based neutrino oscillation experiments [1–4]. The analysis of these experiments requires having good control of nuclear effects. The simple description based on a relativistic Fermi gas (RFG) model does not accurately describe the recent measurements of quasielastic neutrino and antineutrino scattering [5–8]. Mechanisms such as nuclear correlations, final-state interactions, and meson-exchange currents (MECs) may have an impact on the inclusive neutrino charged-current (CC) cross section. In particular, explicit calculations support the theoretical evidence [9–11] for a significant contribution from multinucleon knockout to the CC cross sections ðνμ;μ−Þand ð¯ νμ;μþÞaround and above the quasielastic (QE) peak region, defined by ω¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q2þm2 N p−mN, where ωis the energy transfer and qis the 3-momentum transfer. Recent ab initio calculations [12] of sum rules of weak neutral-current response functions of 12C have also stressed the importance of MECs in neutrino quasielastic scattering. The size of MEC effects is larger than that found in inclusive CC neutrino scattering from the deuteron [13]. The three existing microscopic models that have provided predictions of multinucleon knockout effects in quasielastic neutrino and antineutrino cross sections from 12C for the experimental kinematical settings are those by Martini [14–19], Nieves [10,20–22], and the superscaling analysis (SuSA) model of Refs. [11,23,24]. These three models are based on the Fermi gas, but each one contains different ingredients and approximations to face the problem. The Martini model is based on the nonrelativistic model of Ref. [25], although attempts to improve it using relativistic kinematics have been made. The model includes MEC and pionic correlation diagrams modified to account for the effective nuclear interaction. The interference between direct and exchange diagrams is neglected, in order to reduce the seven-dimensional (7D) integral over the phase space to a two-dimensional (2D) integration. The Nieves model is similar to Martini’s, but most of it is fully relativistic. In this model, the momentum of the initial nucleon in the generic WNNπvertex is fixed to an average value. Under this approximation, the Lindhard function can be factorized inside the integral, leaving only a four-dimensional integration over the momentum of one of the exchanged pions. The direct-exchange interference is neglected as well. The SuSA model includes all the interference terms at the cost of performing a seven-dimensional integration, without any approximation, but the axial part of the MEC is not yet included. It is obvious that these three models should differ numerically because they are different. But a quantitative evaluation of their differences has not been done. Furthermore, the accuracy of the approximations used in these models only can be determined by comparison with an exact calculation for some kinematics. Alternatively, phenomenological approaches have been proposed where two-particle two-hole (2p-2h) effects, estimated by a pure two-nucleon phase-space model, are fitted to the experimental cross section [26,27], while the PHYSICAL REVIEW D 90, 033012 (2014) 1550-7998=2014=90(3)=033012(23) 033012-1 © 2014 American Physical Society nucleon ejection model of Ref. [28] provides a phasespace-based algorithm to generate 2p-2h states in a Monte Carlo implementation. The present paper is a first step toward an extension of the relativistic 2p-2h model of Ref. [29] to the weak sector. We undertake this project with the final goal of including a consistent set of weak MEC in the SuSA approach to CC neutrino reactions [11,23]. The model of Ref. [29] fully described the contribution of 2p-2h states to the transverse response function in electron scattering. Based on the RFG, the model included all 2p-2h diagrams containing two pionic lines (except for nucleon correlations that were included in Ref. [30]), taking into account the quantum interferences between direct and exchange two-body matrix elements. Previous calculations of two-particle emission with MEC in ðe; e0Þinvolved nonrelativistic models [25,31–36]. The first attempts for a relativistic description were made by Dekker [37–39], followed by the model of De Pace et al. [29,40]. The extension of this model to the weak sector requires the inclusion of the axial terms of MEC. Quasielastic neutrino scattering requires one to perform an integral over the neutrino flux. This would considerably increase the computing time of the nuclear response function of Ref. [29] involving 7D integrals of thousands of terms, although improvements were made in Ref. [30] to perform the spin traces numerically. Thus, in this work, we address the problem from a different perspective, focusing first on a careful study of the 7D integral over the 2p-2h phase space as a function of the momentum and energy transfers. Our goal is to provide a comprehensive description of the angular distribution, showing that there is a divergence in the integrand for some kinematics and identifying mathematically the allowed integration intervals. At the same time, we derive a procedure to integrate the angular distribution around the divergence analytically. This procedure allows us to reduce the CPU time considerably. This program is followed first in a pure phase-space domain, without yet including the two-body current. We also sketch the future perspectives opened by this general formalism applied to the calculation of 2p-2h contributions to electroweak responses. In a forthcoming paper, we will provide a full model of weak MEC to compute the complete set of CC neutrino scattering response functions. The paper is organized as follows. In Sec. II, we define the relativistic 2p-2h response and phase space functions. In Sec. III, we review the nonrelativistic description of the 2p-2h integrals, semianalytical expressions that will be used as a check of the relativistic calculations, and some interesting properties of the phase-space integral, such as scaling and asymptotic expansion. In Sec. IV, we address the relativistic phase-space function and asymptotic expansion and show that some numerical problems arise from a straightforward calculation for high q. In Sec. V,we describe the 2p-2h angular distribution in the frozen nucleon approximation and show that this distribution has a divergence for some angles. The divergence is related to the two solutions of the energy conservation for a fixed emission angle. We give kinematical and geometrical explanations of these two solutions. In Sec. VI, we make a theoretical analysis of the angular distribution and find analytically the boundaries of the angular intervals. We get a formula, Eq. (95), for the integral around the divergent angles. In Sec. VII, we present results for the phase-space function with the new integration method. In Sec. VIII,we discuss how this formalism can be applied to the 2p-2h response functions of electron and neutrino scattering. Finally, in Sec. IX, we present our conclusions. II. 2P-2H RESPONSE FUNCTIONS When considering a lepton that scatters off a nucleus transferring 4-momentum Qμ¼ðω;qÞ, with ωthe energy transfer and qthe momentum transfer, one is involved with the hadronic tensor Wμν ¼X fhΨfjJμðQÞjΨiihΨfjJνðQÞjΨiiδðEiþω−EfÞ; ð1Þ where JμðQÞis the electroweak nuclear current operator. In this paper, we take the initial nuclear state as the RFG model ground state, jΨii¼jFi, with all states with momenta below the Fermi momentum kFoccupied. The sum over final states can be decomposed as the sum of one-particle one-hole (1p-1h) plus 2p-2h excitations plus additional channels: Wμν ¼Wμν 1p1hþWμν 2p2hþ ð2Þ In the impulse approximation, the 1p-1h channel gives the well-known response functions of the RFG. Notice that MEC also contribute to these 1p-1h responses; however, here we focus on the 2p-2h channel where the final states are of the type jΨfi¼j10;20;1−1;2−1ið3Þ ji0i¼jp0 is0 it0 iið4Þ jii¼jhisitii;i;i 0¼1;2;ð5Þ where p0 iare momenta of relativistic final nucleons above the Fermi sea, p0 i>k F, with 4-momenta P0 i¼ðE0 i;p0 iÞ, and Hi¼ðEi;hiÞare the 4-momenta of the hole states with hi<k F. The spin indices are s0 iand si, and the isospin is ti,t0 i. In this paper we study the 2p-2h channel in a fully relativistic framework. The corresponding hadronic tensor is given by I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-2 Wμν 2p-2h¼V ð2πÞ9Zd3p0 1d3p0 2d3h1d3h2 m4 N E1E2E0 1E0 2 ×rμνðp0 1;p0 2;h1;h2ÞδðE0 1þE0 2−E1−E2−ωÞ ×Θðp0 1;p 0 2;h 1;h 2Þδðp0 1þp0 2−h1−h2−qÞ; ð6Þ where mNis the nucleon mass, Vis the volume of the system, and we have defined the product of step functions Θðp0 1;p0 2;h1;h2Þ¼θðp0 2−kFÞθðp0 1−kFÞθðkF−h1ÞθðkF−h2Þ: ð7Þ The function rμνðp0 1;p0 2;h1;h2Þis the hadronic tensor for the elementary transition of a nucleon pair with the given initial and final momenta, summed up over spin and isospin, given schematically as rμνðp0 1;p0 2;h1;h2Þ¼1 4X s;t jμð10;20;1;2Þ Ajνð10;20;1;2ÞA; ð8Þ which we write in terms of the antisymmetrized two-body current matrix element jμð10;20;1;2ÞA, to be specified. The factor 1=4accounts for the antisymmetry of the 2p-2h wave function. Finally, note that the 2p-2h response is proportional to V, which is related to the number of protons or neutrons Z¼N¼A=2by V¼3π2Z=k3 F. In this work, we only consider nuclear targets with pure isospin zero. In the case of electrons, the cross section can be written as a linear combination of the longitudinal and transverse response functions defined by RL¼W00 ð9Þ RT¼W11 þW22;ð10Þ whereas additional response functions arise for neutrino scattering, due to the presence of the axial current. The generic results coming from the phase-space obtained here are applicable to all of the response functions. Integrating over p0 2using the momentum delta function, Eq. (6) becomes a nine-dimensional integral, Wμν 2p-2h¼V ð2πÞ9Zd3p0 1d3h1d3h2 m4 N E1E2E0 1E0 2 ×rμνðp0 1;p0 2;h1;h2ÞδðE0 1þE0 2−E1−E2−ωÞ ×Θðp0 1;p 0 2;h 1;h 2Þ;ð11Þ where p0 2¼h1þh2þq−p0 1. After choosing the qdirection along the zaxis, there is a global rotation symmetry over one of the azimuthal angles. We choose ϕ0 1¼0and multiply by a factor 2π. Furthermore, the energy delta function enables analytical integration over p0 1, and so the integral is reduced to seven dimensions. In general, the calculation has to be done numerically. Under some approximations [25,31,32,36], the number of dimensions can be further reduced, but this cannot be done in the fully relativistic calculation. In this paper, we study different methods to evaluate the above integral numerically and compare the relativistic and the nonrelativistic cases. In the nonrelativistic case, we reduce the hadronic tensor to a two-dimensional integral. This can be done when the function rμν only depends on the differences ki¼p0 i−hi,i¼1,2. As we want to concentrate on the numerical procedure without further complications derived from the momentum dependence of the currents, in this paper, we start by setting the elementary function to a constant rμν ¼1. Hence, we focus on the genuine kinematical effects coming from the two-particle–two-hole phase space alone. In particular, the kinematical relativistic effects arising from the energymomentum relation are contained in the energy conservation delta function that determines the analytical behavior of the hadronic tensor, where the energy-momentum relation is E¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2þm2 N p, and in the Lorentz contraction coefficients mN=Ei. Obviously, the results obtained here for constant rμν will be modified when including the two-body physical current. But as the final result is model dependent, it is not possible to disentangle whether the differences found are due to the current model employed or to the approximations (relativistic or not) used to perform the numerical evaluation of the integral. In fact all of the models of 2p-2h response functions should agree at the level of the 2p-2h phase-space integral Fðq; ωÞdefined as Fðq; ωÞ≡Zd3p0 1d3h1d3h2 m4 N E1E2E0 1E0 2 ×δðE0 1þE0 2−E1−E2−ωÞΘðp0 1;p 0 2;h 1;h 2Þ; ð12Þ with p0 2¼h1þh2þq−p0 1. Calculation of this function should be a good starting point to compare and congenialize different nuclear models. III. NONRELATIVISTIC 2P-2H PHASE SPACE A. Semianalytical integration First, we recall the semianalytical method of Ref. [32] that was used later in Refs. [25,29], for instance, to compute the nonrelativistic 2p-2h transverse response function in electron scattering. We shall use this method to check the numerical 7D quadrature both in the relativistic and nonrelativistic cases. We start with the 12-dimensional expression for the phase-space function, Eq. (6), RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-3 Fðq; ωÞ¼Zd3p0 1d3p0 2d3h1d3h2 ×δðE0 1þE0 2−E1−E2−ωÞ ×Θðp0 1;p 0 2;h 1;h 2Þδðp0 1þp0 2−h1−h2−qÞ: ð13Þ The procedure is first to perform the integral over energy. Following Ref. [32], we change variables: l1¼p0 1−h1 kF l2¼p0 2−h2 kFð14Þ x1¼p0 1þh1 2kF x2¼p0 2þh2 2kF :ð15Þ We also define the following nondimensional variables: qF≡ q kFð16Þ ν≡ mNω k2 F :ð17Þ In terms of these variables, the 2p-2h phase-space function is Fðq; ωÞ¼ð2πÞ2k7 FmNZd3l1 l3 1 d3l2 l3 2 ×δðl1þl2−qFÞAðl1;l 2;νÞ;ð18Þ where we use the Van Orden function defined as Aðl1;l2;νÞ¼ l3 1l3 2 ð2πÞ2Zd3x1d3x2δðν−l1·x1−l2·x2Þ ×θ1− x1−l1 2θ1− x2−l2 2 ×θ x1þl1 2 −1θ x2þl2 2 −1:ð19Þ This function was computed analytically in Ref. [32].In this work, we have checked that expression because we found a typo in one of the terms in the original reference (that typographical error does not affect the results of the cited reference). We give in the Appendix the correct result for future reference. Integrating now over the momentum l2, we get Fðq; ωÞ¼ð2πÞ2k7 FmNZd3l1 l3 1jqF−l1j3Aðl1;jqF−l1j;νÞ: ð20Þ The integral over the azimuthal angle ϕ1of l1gives 2π. Finally, changing to the variables x¼l1;y¼jqF−l1j;ð21Þ we obtain Fðq;ωÞ¼ð2πÞ3k7 FmN qFZxmax 0 dx x2ZqFþx jqF−xj dy y2Aðx;y;νÞ;ð22Þ where the maximum value of x(or k1=kF) is obtained from the energy conservation and momentum step functions included implicitly in the function Aðx; y; νÞ. In the Appendix, we derive the inequality x≤xmax ≡1þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1þνÞ p:ð23Þ The two-dimensional integral over the variables x,yhas to be performed numerically. B. Numerical integration The simplicity of the Fermi gas model used in this paper allows us to compute the 2p-2h hadronic tensor as a 7D integral as shown below. Note that in a more sophisticated model where the nuclear distribution details are taken into account, like shell models or the spectral function-based models, some of the numerical problems linked to the particular Jacobian appearing here and in the following section can be avoided, at the price of increasing the number of integrals or sums over shell-model states, thus making the calculations harder. The local Fermi gas used by Nieves et al. is really an average of different Fermi gases at different densities, but the basic Fermi gas equations are the same as here. The hadronic tensor for the elementary 2p-2h transition, Eq. (8), contains the direct and exchange matrix elements of the two-body current operator. If one neglects the interference between the direct and exchange terms, it is possible to express rμν as a function of x,yonly, and one can use the formalism of the above section to reduce the calculation of the 2p-2h hadronic tensor to a 2D integral. In the general case, the interference cannot be neglected. It is then necessary to evaluate a 7D integral numerically. Thus, in this work, we also compute the phase-space function, Eq. (12), numerically. This will allow us first to check the numerical procedures by comparison with the semianalytical method of the previous section; second, to determine the number of integration points needed to obtain accurate results, and, third, to optimize the computational effort. This numerical study will be very useful when including actual nuclear currents. Starting with Eq. (12), we compute the integrand for ϕ0 1¼0(the azimuthal angle of p0 1) and multiply by 2π. Then, we use the δof energies to integrate over the variable p0 1, for fixed momenta h1and h2and emission angle θ0 1.To do so, we first define the total momentum of the two particles that is fixed by momentum conservation: I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-4 p0¼p0 1þp0 2¼h1þh2þq:ð24Þ We then change from variable p0 1to variable E0: E0¼E0 1þE0 2¼p0 1 2 2mNþðp0−p0 1Þ2 2mN :ð25Þ By differentiation with respect to p0 1, we obtain  dp0 1 dE0¼mN jp0 1−p0 2·ˆ p0 1j;ð26Þ where ˆ p0 1¼p0 1=p0 1is the unit vector in the direction of the first particle. Integrating now over E0, energy conservation is obtained as E0¼E1þE2þω:ð27Þ Substituting Eq. (25) a second degree equation is obtained for p0 1: 2p0 1 2þp02−2p0·p0 1¼2mNE0:ð28Þ So, we see that there can be two values of the nucleon momentum compatible with energy conservation, for fixed a emission angle. We denote the two solutions by p0 1ðÞ ¼1 22 4vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2−4p02 2−mNE0 s3 5;ð29Þ where we have defined v≡p0·ˆ p0 1:ð30Þ Using this result, we finally evaluate the phase-space function as the 7D integral Fðq; ωÞ¼2πZd3h1d3h2dcos θ0 1 ×X α¼ p0 1 2mN jp0 1−p0 2·ˆ p0 1jΘðp0 1;p 0 2;h 1;h 2Þjp0 1¼p0 1ðαÞ ; ð31Þ where the sum inside the integral runs over the two solutions p0 1ðÞ of the energy conservation equation. C. Asymptotic expansion It is of interest to quote the limit ω→∞because it can also be used for testing the numerical integration. The most useful case applies for kF,q≪ω, when one can neglect all momenta compared with the energy transfer ω, because the phase-space integral can be performed analytically. Note that for the scattering reactions of interest, this limit is not physical (because ω<q, namely, spacelike, for real particles). It is only a mathematical property of the function F, which is well defined for all the ωvalues not only the physical ones. We start writing the momentum of the first particle, Eq. (29),as p0 1¼v 21 2ffiffiffiffi D p;ð32Þ with the discriminant D¼v2−2p02þ4mNE0:ð33Þ The limit ω→∞can be obtained by noticing that vand p0 do not depend on ωbut only on the momenta h1,h2, and q and that E0¼E1þE2þω∼ω. Then, D∼4mNω;ð34Þ and the positive solution for the momentum is p0 1∼ffiffiffiffiffiffiffiffiffiffi mNω p:ð35Þ That is, each nucleon exits the nucleus taking half of the available energy. On the other hand, using Eq. (29), we note that the denominator in Eq. (31) can be written as p0 1−p0 2·ˆ p0 1¼ ffiffiffiffi D p∼2ffiffiffiffiffiffiffiffiffiffi mNω p:ð36Þ Then, Fðq; ωÞ→ ω→∞Faðq; ωÞ ≡2πZd3h1d3h2dcos θ0 1 mN 2ffiffiffiffiffiffiffiffiffiffi mNω p ¼4π4 3πk3 F2mN 2ffiffiffiffiffiffiffiffiffiffi mNω p:ð37Þ Thus, for high energy, the nonrelativistic phase-space function increases as ffiffiffiffi ω p. We shall see in the next section a different behavior in the relativistic case. D. Nonrelativistic results In Fig. 1, we show the nonrelativistic phase-space function Fðq; ωÞas a function of ωfor three typical values of the momentum transfer q¼300, 400, and 500 MeV=c. The Fermi momentum is kF¼225 MeV=c. We compare the two computational methods: the semianalytical of Eq. (22) and the numerical 7D integration of Eq. (31). The semianalytical result is essentially exact because we can choose a very small integration step for the 2D integral (using steps of 0.02 or 0.01, the results do not change in the scale of the figure). However, the 7D integral is computationally time consuming, and the integration step cannot be very small. Here, we compute the integral with a RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-5 “straightforward”method, as an average over a grid with ntotal integration points, uniformly distributed. For large n, the straightforward integration should give results similar to the Monte Carlo methods used in previous calculations [29,32]. The number of points chosen for this calculation was 25 for the variable θ0 1(although it can safely be reduced to 16) and 16 for each one of the remaining dimensions. In total, the number of points is n¼0.42 ×109. This is well above the maximum number n¼106–107, typical of previous calculations [29,32] performed using Monte Carlo techniques. Using ten integration points in each dimension gives very similar results, except for some ωregions where the numerical error is manifested in an apparently slightly less smooth behavior. Increasing the number of points would improve the results; however, this is not practical because the inclusion of the two-body current would make the calculation too slow. The semianalytical and numerical results are quite similar, the difference between them being of a few percent. For comparison, we also show the asymptotic limit ω→∞, computed using the analytical expression in Eq. (37), which is proportional to ffiffiffiffi ω p. We see that for high ωthe function Fðq; ωÞbecomes close to the asymptotic value Faðq; ωÞ.Forq¼300 MeV=c, the asymptotic value is almost reached at the photon point ω¼q. When q increases, so does the distance to the asymptote at the photon point. IV. RELATIVISTIC 2P-2H PHASE SPACE Having two independent calculations of the phase-space function Fðq; ωÞin the nonrelativistic limit, we now consider the case of the fully relativistic calculation as given by Eq. (12). This involves adding the Lorentzcontraction mN=E factors and using relativistic kinematics in the energy δfunction. Following the scheme of the previous section, again azimuthal symmetry allows one to fix ϕ0¼0and multiply by 2π. To integrate over p0 1,we change to the variable E0¼E0 1þE0 2¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p0 1 2þm2 N qþffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðp0−p0 1Þ2þm2 N q; ð38Þ where again p0¼h1þh2þqis the final momentum for a fixed pair of holes. By differentiation, we arrive at the following Jacobian:  dp0 1 dE0¼ p0 1 E0 1 −p0 2·ˆ p0 1 E0 2 −1 :ð39Þ The nonrelativistic Jacobian of Eq. (26) is recovered for low energies E0 1≃E0 2≃mN. As before, integration over E0 gives E0¼E1þE2þω, and the phase-space function becomes Fðq; ωÞ¼2πZd3h1d3h2dθ0 1sin θ0 1 m4 N E1E2 ×X α¼ p0 1 2  p0 1 E0 1 −p0 2·ˆ p0 1 E0 2 Θðp0 1;p 0 2;h 1;h 2Þ E0 1E0 2p0 1¼p0 1ðαÞ ; ð40Þ where again the sum inside the integral runs over the two solutions p0 1ðÞ of the energy conservation equation p0 1ðÞ ¼1 ~ b~ a~ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ a2−~ bm2 N q:ð41Þ FIG. 1 (color online). Nonrelativistic phase-space function calculated for ω¼300, 400, 500 MeV, using a numerical and a semianalytical approach. The number of points used in two numerical integrations is indicated in the plot. We also show the asymptotic function for comparison. I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-6 The definitions of the quantities ~ a,~ b, and ~ vare given in the Appendix. Note that there is a difference between our Jacobian in Eq. (40) and that given in Eqs. (15–17) of Ref. [26]. The relativistic approach is more involved than the nonrelativistic one because it requires taking the square twice in the original equation to eliminate the squared roots in the energies. This can introduce spurious solutions for p0 1 depending on the kinematics, which have to be eliminated from the above sum in the numerical procedure. This is not a trivial task, and details are provided in the Appendix. The appearance of spurious solutions is a difference between the relativistic and nonrelativistic methods. A second one will be discussed below in relation to a divergence of the integrand. Therefore, the relativistic calculation is very involved, and it cannot be derived by simply extending the nonrelativistic code. We devote the rest of this section to explain in detail how to get the fully relativistic answers. A. Relativistic asymptotic expansion Although it is not possible to derive a semianalytical expression for Fðq; ωÞas in the nonrelativistic case, it is still possible to take the limit ω→∞and obtain an analytical result. As in the nonrelativistic case, we assume kF,q≪ω. If we add the condition mN≪ω, we can neglect the momenta and energies of the two holes and write E0∼ωp0∼q:ð42Þ We can also compute the quantities with tildes that appear in the solution of the energy conservation (see the Appendix), obtaining ~ a∼ω 2ð43Þ ~ v∼q·ˆ p0 1 2ωð44Þ ~ b∼1:ð45Þ Then, the discriminant of Eq. (41) becomes ~ a2−~ bm2 N∼ω2 4−m2 N∼ω2 4:ð46Þ Therefore, the allowed solution of the energy conservation equation is p0 1∼q·ˆ p0 1 4þω 2∼ω 2:ð47Þ Thus, in this limit, each nucleon carries half the total energy and momentum, E0 1∼p0 1∼E0 2∼p0 2∼ω 2:ð48Þ Now, the Jacobian, the denominator in Eq. (40), can be computed as d≡ p0 1 E0 1 −p0 2·ˆ p0 1 E0 2¼1−ðp0−p0 1Þ·ˆ p0 1 E0 2 ∼1þp0 1 E0 2 ∼2: ð49Þ Collecting Eqs. (47),(48), and (49), the integrand in Eq. (40) becomes p0 1 2 d m4 N E1E2E0 1E0 2 ∼ω2 8 m4 N m2 Nω2=4¼m2 N 2:ð50Þ Finally, performing the integral, we obtain the following asymptotic expression: Fðq; ωÞ→ ω→∞Faðq; ωÞ¼4π4 3πk3 F2m2 N 2:ð51Þ In contrast with the nonrelativistic behavior, which increases monotonically as ffiffiffiffi ω p, the relativistic result (37) goes to a constant. The Lorentz contraction factors E=mNbalance the ω2behavior coming from the phase space. This analytical result for high ωwill be useful for comparison of the numerical results for high ω. B. Relativistic straightforward calculation Before going to the high-qregion, we first check the relativistic phase-space function results by comparison with the nonrelativistic counterpart. Both should agree for low energy. We proceed by performing a straightforward numerical integration of Eq. (40) as in the nonrelativistic case. In Fig. 2, we show the results of this comparison for q¼300, 400, and 500 MeV=c. We also show both the numerical and “exact”(i.e., using the semianalytical formula) nonrelativistic function Fðq; ωÞ. A uniform distribution with ten points for each dimension is employed in the 7D integrations. As expected, relativistic and nonrelativistic results agree at low energy. The relativistic effects consist of a reduction of the strength at high energy. The amount of this reduction is very small for q¼300 MeV=c, where the nonrelativistic approximation can be safely used, and increases with q, reaching about 15% for q¼500 MeV=c. Thus, for low q, we agree that a number of ∼107points is adequate for numerical integration purposes. In Fig. 5, the asymptotic limit Faðq; ωÞ of the relativistic phase space, Eq. (51), is also shown. For these low qvalues, Fðq; ωÞis still far below the asymptote. Larger relativistic effects are expected for intermediate to large momentum transfer. In Fig. 3, we display Fðq; ωÞ RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-7 for q¼700, 1000, and 1500 MeV=c, compared with the exact nonrelativistic results. Using straightforward 7D integration, we need to increase the number of points to 16 for each dimension in order to reach some stability of the results shown in Fig. 3. However, we find that full convergence would need more points. In fact, for q¼700 MeV=c, a small deviation with respect to the exact result can be noticed at low ω. This deviation increases with qand turns into a prominent structure with a“shoulder”shape for q¼1500 MeV=c. One could be tempted to attribute this effect to relativity. But this is not the case because the same behavior is also present in a nonrelativistic numerical calculation. As we will explain below, this is just a consequence of the inadequacy of the straightforward integration method at high q. This problem affects only the inner integral over θ0 1. Below, we address this issue by a detailed analysis of the θ0 1dependence of the integrand. V. ANGULAR DISTRIBUTION A. Frozen phase-space function We start fixing a value of q¼3GeV=c that is high enough to amplify the misbehavior found above and also allows us to simplify the analysis that follows. In fact, we note that for very high q≫kF, all of the hole momenta h1, h2could safely be neglected inside the integral as a first approximation. Since this implies that the initial particles are at rest, we denote this limit the “frozen nucleon approximation.”In particular, the energies of the holes can be substituted by the nucleon mass in the δfunction, FIG. 3 (color online). Relativistic phase space function for q¼700, 1000, 1500, calculated using straightforward integration compared with the nonrelativistic calculation using the semianalytical approach. FIG. 2 (color online). Relativistic phase space function for q¼300, 400, 500, calculated using straightforward integration, compared with the nonrelativistic calculation using the semianalytical approach. I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-8 Fðq; ωÞ∼Zd3h1d3h2d3p0 1δðE0 1þE0 2−ω−2mNÞ ×Θðp0 1;p 0 2;0;0Þm2 N E0 1E0 2 ;ð52Þ where p0 2¼q−p0 1. Because the integrand does not depend on the hole momenta, one can directly integrate out those variables, Fðq; ωÞ∼4 3πk3 F2Zd3p0 1δðE0 1þE0 2−ω−2mNÞ ×Θðp0 1;p 0 2;0;0Þm2 N E0 1E0 2 :ð53Þ Now, the integral over p0 1can be done analytically as before using the delta function, with the same Jacobian evaluated for h1¼h2¼0. The integral over ϕ0 1gives again a factor 2π, Fðq; ωÞ∼2πm2 N4 3πk3 F2Zdθ0 1sin θ0 1 ×X α¼ p0 1 2  p0 1 E0 1 −p0 2·ˆ p0 1 E0 2 Θðp0 1;p 0 2;0;0Þ E0 1E0 2p0 1¼p0 1ðαÞ :ð54Þ Thus, in this approximation, the phase-space function is reduced to a one-dimensional integral over the emission angle θ0 1, which has to be performed numerically. The frozen nucleon approximation represents just a particular case of the mean-value theorem for the integral over h1,h2. We denote with a bar the quantities computed by the mean-value theorem. Thus, we define the barred phase-space function ¯ Fðq; ωÞ¼4 3πk3 F2Zd3p0 1δðE0 1þE0 2−ω−E1−E2Þ ×Θðp0 1;p 0 2;h 1;h 2Þm4 N E1E2E0 1E0 2 ;ð55Þ where p0 2¼h1þh2þq−p0 1, and ðh1;h2Þare a pair of fixed momenta below the Fermi sea. Going further, we will later turn to the question of how to choose the average nucleon momenta h1,h2for low q. For high q, we expect this function not to depend too much on the chosen values. So, at this point, we restrict our study to ¯ Fðq; ωÞin the frozen nucleon approximation, i.e., for h1¼h2¼0. B. Numerical analysis We have computed ¯ Fðq; ωÞin the frozen nucleon approximation using 100 points to perform the numerical integral over the emission angle θ0 1. Results are shown in Fig. 4. A misbehavior due to numerical error is now evident. The reason for the appearance of discontinuities by numerical integration becomes apparent by examining the angular dependence of the integrand. We define the angular distribution function, for fixed values of ðq; ωÞand h1,h2,as Φðθ0 1Þ¼sin θ0 1Zp0 1 2dp0 1δðE1þE2þω−E0 1−E0 2Þ ×Θðp0 1;p 0 2;h 1;h 2Þm4 N E1E2E0 1E0 2 ; ¼X α¼ m4 Nsin θ0 1p0 1 2Θðp0 1;p 0 2;h 1;h 2Þ E1E2E0 1E0 2 p0 1 E0 1 −p0 2·ˆ p0 1 E0 2p0 1¼p0 1ðαÞ ;ð56Þ where once more p0 2¼h1þh2þq−p0 1, such that the phase-space function is obtained by integration over the emission angle θ0 1: ¯ Fðq; ωÞ¼4 3πk3 F2 2πZπ 0 dθ0 1Φðθ0 1Þ:ð57Þ The function Φðθ0 1Þthus measures the distribution of final nucleons as a function of the angle θ0 1. This function is computed analytically, given by the integrand in Eq. (54). Results for Φðθ0 1Þare shown in Fig. 5for h1¼h2¼0; q¼3GeV=c; and for the three values of ω¼1800, 2000, FIG. 4 (color online). Phase-space function for q¼3and 0.5GeV=c, computed using the frozen nucleon approximation for fixed hole momenta h1¼h2¼0, using 100 integration points in an emission angle. RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-9 results for all energies. For high momentum transfer q¼3GeV=c, they are completely different. The nonrelativistic function is pushed toward higher energies due to the quadratic momentum dependence of the nonrelativistic kinetic energy. Thus, for q¼1.5GeV=c, the relativistic results are above (below) the nonrelativistic ones for low (high) energy. For q¼3GeV=c, the relativistic results are above for all ωvalues allowed. In all cases, Fðq; ωÞis below the asymptotic value, Eq. (51), also shown in Fig. 11. In Fig. 12, we get deeper insight into the size of relativistic effects. There, we show results for Fðq; ωÞcomputed using relativistic kinematics only, but without including the relativistic Lorentz-contraction factors mN=E in particles and holes. The results increase a lot with respect to the nonrelativistic ones. This is related to the fact pointed out after Eq. (51) for the asymptotic limit of the relativistic phase-space integral. Without the Lorentz factors, the function Fðq; ωÞwould increase as ω2. This seems to indicate that in order to “relativize”a nonrelativistic 2p-2h model, implementing only relativistic kinematics is not sufficient, since it goes in the wrong direction. In fact, results in Fig. 12 show that the effects coming solely from the relativistic kinematics lead to differences even larger than the discrepancy between the nonrelativistic and the fully relativistic calculations. Therefore, it is essential also to include the Lorentz factors mN=E. Note that the behavior of relativistic effects in the 1p-1h channel goes in the opposite direction to the one discussed FIG. 11 (color online). Total phase-space function for three values of the momentum transfer. The number of integration points in each dimension in the hole variables is indicated by n. The number of integration points over the emission angle θ0 1is indicated as m. We also show the nonrelativistic, exact result and the relativistic asymptotic value. FIG. 12 (color online). Effect of implementing relativistic kinematics in a nonrelativistic calculation of Fðq; ωÞ. Solid lines: nonrelativistic result. Thick dotted lines: relativistic kinematics only without the relativistic factors mN=E. Thin dashed lines: fully relativistic result. I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-16 here in the 2p-2h channel. In fact, in Ref. [41] it was shown that implementing relativistic kinematics without the mN=E factors in the nonrelativistic 1p-1h response function gives a result closer to the exact relativistic response function (see Fig. 33 of Ref. [41]). In Figs. 13 and 14, we present the results of a study of the validity of the frozen nucleon approximation to compute Fðq; ωÞin a range of momentum transfers. This approximation was introduced for high momentum transfer q¼3GeV=c, neglecting the momenta of the two holes inside the 7D integral, thus reducing it to a one-dimensional (1D) integral over the emission angle θ0 1. In Fig. 14, the momentum transfer is still high, and the frozen nucleon approximation remains valid. In Fig. 13, the values of qare not so large, and one could think that the frozen nucleon approximation is not valid. However, the results of Fig. 13 demonstrate that it is still a good approximation for moderate momentum transfer except for very low energy transfer, where the function Fðq; ωÞis small. This is a promising result; if the frozen nucleon approximation could be extended to the full response functions when including the nuclear current, this would mean that the 2p-2h cross section could be approximated by 1D integrals over the emission angle, which would be easy and fast to compute. In particular, calculations of this kind could be implemented in existing Monte Carlo codes. To illustrate the reasons why the frozen nucleon approximation works for moderate momentum transfer, we present Figs. 15,16, and 17. We compare Fðq; ωÞwith the barred phase-space function ¯ Fðq; ωÞ, defined in Eq. (55), computed for several ðh1;h2Þconfigurations. The averagemomentum approximation is similar to the frozen nucleon approximation in the sense that the two hole momenta h1,h2are set to a constant inside the integral. For a pair FIG. 13 (color online). Relativistic phase-space function Fðq; ωÞcompared with the frozen nucleon approximation for low to intermediate momentum transfer. FIG. 14 (color online). Relativistic phase-space function Fðq; ωÞcompared with the frozen nucleon approximation for high momentum transfer. RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-17 configuration ðh1;h2Þ, the function ¯ Fðq; ωÞgives the contribution of such a pair to the phase-space function, multiplied by V2 F, where VF¼4πk3 F=3is the volume of the Fermi sphere. The total Fðq; ωÞis the sum of the contributions from all of the pairs or, equivalently, the average of all of the barred functions ¯ Fðq; ωÞover the different pair configurations. In Fig. 18, we show the geometry for the configurations used in Figs. 15–17. For low values of the momenta h1,h2, the frozen nucleon approximation should be a good approximation to the average phase space. For larger values of the momenta, we find pairs of configurations with opposite total momentum p¼h1þh2that contribute above or below the average in approximately equal footing, so they do not change the mean value very much. In the first example, Fig. 15, we show the contribution of two pairs of nucleons with the same momentum h1¼h2¼200 MeV=c, and both parallel, pointing upward (U) and downward (D) with respect to the z axis, that is, the direction of q. The contribution of the UU configuration is smaller than average, while the DD is larger. This is so because in the UU case the total momentum p0in the final state is large. By momentum conservation, the momenta p0 1and p0 2must also be large. Therefore, these states need a large excitation energy, and they start to contribute for high ωtransfer. In the DD configuration, the total momentum p0is small, so the final momenta p0 1and p0 2can also be small, will small FIG. 15 (color online). Relativistic phase-space function Fðq; ωÞcompared with the average-momentum approximation ¯ Fðq; ωÞfor a pair of nucleons with momentum 200 MeV=c. In the UU configuration, both nucleons move along q(up). In the DD configuration, both move opposite to q(down). FIG. 16 (color online). Relativistic phase-space function Fðq; ωÞcompared with the average-momentum approximation ¯ Fðq; ωÞfor a pair of nucleons with momentum 200 MeV=c pointing in opposite directions (total momentum equal to zero). In the UD configuration one moves along q(U) and the other opposite to q(D). In the T,−Tconfiguration, one moves in the x direction and the other in the −xdirection. I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-18 excitation energy. Therefore, they start to contribute at very low ω. In the example of Fig. 16, two antiparallel configurations are shown. In the UD case, one nucleon is moving upward and the other downward the zaxis with total momentum zero of the pair. This situation is similar to that of a pair of highly correlated nucleons with large relative momentum [42]. Since the total momentum is zero, the final 2p-2h state has total momentum q, exactly the same that it would have in the frozen nucleon approximation. Therefore, the contribution of this configuration is similar to the average. The same conclusions can be drawn in the case of the configuration T,−T, with one nucleon moving along the x axis (transverse direction) and the other along −xwith opposite momentum. The contribution of this pair is exactly the same as that of the UD configuration in the total phase-space function. Finally, we show in Fig. 17 two intermediate cases that are neither parallel nor antiparallel configurations. They consist of two pairs of transverse nucleons moving along mutually perpendicular directions. In the first case, we consider a Tnucleon and a second T0nucleon moving in the yaxis out of the scattering plane. The contribution of the TT0pair is large, while the one of the opposite case, −T, −T0, is small. On the average, they are close to the total result. VIII. PERSPECTIVES ON THE CALCULATION OF 2P-2H ELECTROWEAK RESPONSE FUNCTIONS The next step in our project of an exact evaluation of the relativistic 2p-2h electroweak response functions in the Fermi gas model initiated with the approach in the present paper would be to apply it to a more realistic situation, i.e., electron and neutrino scattering. The 2p-2h states can be excited by two-body MEC operators, involving exchange of an intermediate meson between two nucleons. A complete calculation, including all of the MEC diagrams FIG. 17 (color online). Relativistic phase-space function Fðq; ωÞcompared with the average-momentum approximation ¯ Fðq; ωÞfor a pair of nucleons with momentum 200 MeV=c pointing in perpendicular directions. In the T,T0configuration, one moves along x(T) and the other along y(T0). In the −T,−T0 configuration, they move along the −xand −ydirections, respectively. FIG. 18 (color online). Geometry employed for emission of a pair of nucleons with momenta parallel (cases UU,DD), antiparallel (cases UD,T−T), and perpendicular (cases TT0,−T−T0). RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-19 with one pion exchange, is out of the scope of the present paper and will be reported in a forthcoming publication. However, in this section, we discuss the perspectives opened by the formalism presented here. One question can arise on why the integration problem related to the divergence in the angular distribution, which is the central issue in this work, does not appear in other approaches. The models developed by Martini [9] and Nieves [10] neglect the direct-exchange interference terms in the hadronic tensor. In this approximation, in the nonrelativistic case, the change of variables introduced in Eqs. (14) and (15) reduces the integration to two dimensions; the integration variables in this case are proportional to the magnitude of the transferred momenta to the two nucleons. In the relativistic model of Ref. [10], an additional approximation is made for the total WNN interaction vertex, where the dependence on the initial nucleon momentum is neglected by fixing it to an average value over the Fermi sea. This trick allows one to factorize the two Lindhard functions linked to the two nucleon loops in the many-body diagrams. Thus, two approximations are required in this case to reduce the calculation to a four-dimensional integral over the 4-momentum of one of the exchanged pions. The exact calculation including all the terms in the hadronic tensor (direct, exchange, and interference) requires the complete 7D integral. Obviously, a change of variables can be made to eliminate the divergence. One possibility is to make the change θ0 1→ffiffiffiffiffiffiffiffiffiffiffi gðθ0 1Þ p, where gðθ0 1Þis the function defined in Eq. (76). This corresponds to the change of variables made in Sec. VI B to integrate analytically around the divergence. The standard way to handle this problem in the Monte Carlo generators [28] is to compute the 2p phasespace angular distribution in the c.m. system of the final nucleons, because it is angular independent, although Pauli blocking can forbid some angular regions. A transformation to the laboratory system would give exactly the same distribution as considered in this paper. Linked to this, a further possibility that we are presently investigating would be to integrate over the c.m. emission angle instead of the Lab one considered in this work. This procedure would have the advantage of being free of the divergence coming from the Jacobian but has the drawback of requiring the performance of a boost back to the Lab system for each pair of holes ðh1;h2Þ. One should perform a full calculation with both approaches to see the advantages of each one in terms of CPU time. One of the main problems associated with a complete, exact calculation of the 2p-2h response functions is the computational time required when the full current is included. One of the outcomes of this work is the possibility opened by considering what we called the frozen nucleon approximation to compute the integral over the two holes. The validity of this approximation must be verified in the complete calculation. If the approximation is found to be accurate enough, then the calculation of the 2p-2h cross section could be done without much difficulty and could be easily implemented in Monte Carlo generators. The verification of this approximation is one of the goals of our future work. Preliminary results obtained with the seagull diagrams show that the approximation is valid for this set of diagrams. The Monte Carlo generators must not perform the integration over the outgoing final state but instead must keep these momenta explicitly because one is interested in generating a full final state to be propagated. The integration performed here is only needed for the inclusive 2p-2h cross section, which cannot be separated from the measured QE cross section if the final nucleons are not detected. With our model, there is the possibility to generate angular distributions of nucleon 2p-2h states produced by MEC, fully compatible with the inclusive 2p-2h cross sections. They could be useful for the Monte Carlo generators. IX. CONCLUSIONS We have performed a detailed study of the two-particle– two-hole phase-space function, which is proportional to the nuclear two-particle emission response function for constant current matrix elements. To obtain physically meaningful results, one should include a model for the two-body currents inside the integral. However, the knowledge gained here by disregarding the operator and focusing on the purely kinematical properties has been of great help in optimizing the computation of the 7D integral appearing in the 2p-2h response functions of electron and neutrino scattering. The frozen nucleon approximation, that is, neglecting the momenta of the initial nucleons for high momentum transfer, has allowed us to focus on the angular distribution function. We have found that this function has divergences for some angles. Our main goal has been to find the allowed angular regions and to integrate analytically around the divergent points. The CPU time of the 7D integral has been reduced significantly. The relativistic results converge to the nonrelativistic ones for low energy transfer. We are presently working on an implementation of the present method with a complete model of the MEC operators. ACKNOWLEDGMENTS This work was supported by DGI (Spain), Grants No. FIS2011-24149 and No. FIS2011-28738-C02-01; by the Junta de Andalucía (Grants No. FQM-225 and No. FQM-160); by the Spanish Consolider-Ingenio 2010 programmed CPAN, in part (M. B. B.) by the INFN project MANYBODYand in part (T. W. D.) by U.S. Department of Energy under Cooperative Agreement No. DE-FC0294ER40818. C. A. is supported by a CPAN postdoctoral contract. I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-20 APPENDIX A: CALCULATION of xmax Here, we derive the upper limit of the integral over x in Eq. (22). We first note that the function Aðx; y; νÞinside the integral contains the energy delta function δðωþE1þ E2−E0 1−E0 2Þand the step function θðkF−h1ÞθðkF−h2Þ. This implies that E0 1≤E0 1þE0 2¼ωþE1þE2≤ωþ2EF:ðA1Þ Therefore, p0 1 2 2mN ≤ωþ2k2 F 2mN :ðA2Þ Taking the square root and rearranging, one has p0 1≤kFffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2þ2mNω=k2 F q:ðA3Þ Recalling now the definition of the nondimensional variable ν¼mNω=k2 F, we have p0 1 kF ≤ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1þνÞ p:ðA4Þ Finally, using this inequality in the definition of the x variable, one finds that x¼ p0 1−h1 kF ≤p0 1þh1 kF ≤p0 1þkF kF ≤1þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1þνÞ p: ðA5Þ APPENDIX B: FUNCTION Aðl1;l2;νÞ The function Aðl1;l 2;νÞwas computed analytically in Ref. [32]. In this work, we have repeated the analytical calculation, and we have found a typographical error (a minus sign) in that reference. Although the demonstration and numerical results of Ref. [32] are correct, taking into account the given error is essential. For completeness, and because the error can mislead the reader, we write in this Appendix the correct final expression with the slightly different notation used by us. We write the function Aas the sum of 16 terms, Aðl1;l 2;νÞ¼X 4 i¼1X 4 j¼1 Aijðl1;l 2;νÞ;ðB1Þ where the Aij functions have the symmetry Aijðl1;l 2;νÞ¼Ajiðl2;l 1;νÞ:ðB2Þ Thus, we only need to give the analytical expressions for the diagonal and the upper-half off-diagonal ij elements: A11 ¼l1l2 C3 11 3! þðl1þl2ÞC4 11 4! þC5 11 5! θðC11Þ C11 ≡ν−l2 1 2−l1−l2 2 2−l2 A12 ¼l1l2 C3 12 3! þðl2−l1ÞC4 12 4! −C5 12 5! θðC12Þθðl2−2Þ C12 ≡ν−l2 1 2−l1−l2 2 2þl2 A13 ¼−l1l2 C3 13 3! þðl1þl2ÞC4 13 4! þC5 13 5! θðC13Þθð2−l2Þ C13 ≡ν−l2 1 2−l1þl2 2 2−l2 A14 ¼l1 C3 14 3! þC4 14 4! l2 2θðC14Þθð2−l2Þ C14 ≡ν−l2 1 2−l1 A22 ¼l1l2 C3 22 3! −ðl1þl2ÞC4 22 4! þC5 22 5!  ×θðC22Þθðl1−2Þθðl2−2Þ C22 ≡ν−l2 1 2þl1−l2 2 2þl2 A23 ¼−l1l2 C3 23 3! þðl1−l2ÞC4 23 4! −C5 23 5!  ×θðC23Þθðl1−2Þθð2−l2Þ C23 ≡ν−l2 1 2þl1þl2 2 2−l2 A24 ¼l1 C3 24 3! −C4 24 4! l2 2θðC24Þθðl1−2Þθð2−l2Þ C24 ≡ν−l2 1 2þl1 A33 ¼l1l2 C3 33 3! þðl1þl2ÞC4 33 4! þC5 33 5!  ×θðC33Þθð2−l1Þθð2−l2Þ C33 ≡νþl2 1 2−l1þl2 2 2−l2 A34 ¼−l1 C3 34 3! þC4 34 4! l2 2θðC34Þθð2−l1Þθð2−l2Þ C34 ≡νþl2 1 2−l1 A44 ¼l2 1l2 2 ν3 3! θðνÞθð2−l1Þθð2−l2Þ:ðB3Þ Note that the equivalent function A13 in Ref. [32] (denoted F3) has a missing global minus sign. RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-21 APPENDIX C: SOLUTIONS OF RELATIVISTIC ENERGY CONSERVATION Given q,ωand fixing the momenta of the two holes h1, h2, the total energy and momentum of the two particles is also fixed by E0¼E1þE2þω p0¼h1þh2þq: For fixed emission angles of the first particle, ϕ0 1¼0and θ0 1, the value of p0 1is restricted by momentum conservation p0 2¼p0−p0 1and energy conservation, E0 2¼E0−E0 1.In fact, taking the square of the last equation, we should solve E0 2 2¼ðE0−E0 1Þ2:ðC1Þ Having squared, we have introduced spurious solutions with E0−E0 1<0, which should be thrown away. Expanding the right-hand side, using the energymomentum relation in the squared energies, and rearranging terms, we arrive at the equivalent equation E0 1¼~ aþ~ vp0 1;ðC2Þ where we have defined ~ a¼E02−p02 2E0ðC3Þ ~ v¼p0·ˆ p0 1 E0:ðC4Þ Taking the square of Eq. (C2) and again using the energymomentum relation, we arrive at the second-degree equation for p0 1, ~ bp0 1 2−2~ a~ vp 0 1þðm2 N−~ a2Þ¼0;ðC5Þ where we have defined ~ b¼1−~ v2:ðC6Þ Note that Eq. (C2) provides an alternative way to compute the energy E0 1once p0 1is known. It also is valid as a check of the solution. However, caution is needed because taking the square of Eq. (C2) introduces spurious solutions with ~ aþ~ vp0 1<0that should be disregarded. [1] H. Gallagher, G. Garvey, and G. P. Zeller, Annu. Rev. Nucl. Part. Sci. 61, 355 (2011). [2] J. A. Formaggio and G. P. Zeller, Rev. Mod. Phys. 84, 1307 (2012). [3] J. G. Morfin, J. Nieves, and J. T. Sobczyk, Adv. High Energy Phys. 2012, 934597 (2012). [4] L. Alvarez-Ruso, Y. Hayato, and J. Nieves, New J. Phys. (to be published). [5] A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), Phys. Rev. D 81, 092005 (2010). [6] A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), Phys. Rev. D 88, 032001 (2013). [7] G. A. Fiorentini et al. (MINERvA Collaboration), Phys. Rev. Lett. 111, 022502 (2013). [8] K. Abe et al. (T2K Collaboration), Phys. Rev. D 87, 092003 (2013). [9] M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 80, 065501 (2009). [10] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev. C83, 045501 (2011). [11] J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, and C. F. Williamson, Phys. Lett. B 696, 151 (2011). [12] A. Lovato, S. Gandolfi, J. Carlson, S. C. Pieper, and R. Schiavilla, Phys. Rev. Lett. 112, 182502 (2014). [13] G. Shen, L. E. Marcucci, J. Carlson, S. Gandolfi, and R. Schiavilla, Phys. Rev. C 86, 035503 (2012). [14] M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 81, 045502 (2010). [15] M. Martini, M. Ericson, and G. Chanfray, Phys. Rev. C 84, 055502 (2011). [16] M. Martini, M. Ericson, and G. Chanfray, Phys. Rev. D 85, 093012 (2012). [17] M. Martini, M. Ericson, and G. Chanfray, Phys. Rev. D 87, 013009 (2013). [18] M. Martini and M. Ericson, Phys. Rev. C 87, 065501 (2013). [19] M. Martini and M. Ericson, arXiv:1404.1490v1. [20] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B707, 72 (2012). [21] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B721, 90 (2013). [22] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013). [23] J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. Lett. 108, 152501 (2012). [24] G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Lett. B 725, 170 (2013). [25] W. M. Alberico, M. Ericson, and A. Molinari, Ann. Phys. (N.Y.) 154, 356 (1984). [26] O. Lalakulich, K. Gallmeister, and U. Mosel, Phys. Rev. C 86, 014614 (2012). [27] O. Lalakulich, U. Mosel, and K. Gallmeister, Phys. Rev. C 86, 054606 (2012). [28] J. T. Sobczyk, Phys. Rev. C 86, 015504 (2012). [29] A. De Pace, M. Nardi, W. M. Alberico, T. W. Donnelly, and A. Molinari, Nucl. Phys. A726, 303 (2003). I. RUIZ SIMO et al. PHYSICAL REVIEW D 90, 033012 (2014) 033012-22 [30] J. E. Amaro, C. Maieron, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. C 82, 044601 (2010). [31] T. W. Donnelly, J. W. Van Orden, T. De Forest, Jr., and W. C. Hermans, Phys. Lett. 76B, 393 (1978). [32] J. W. Van Orden and T. W. Donnelly, Ann. Phys. (N.Y.) 131, 451 (1981). [33] W. M. Alberico, A. De Pace, A. Drago, and A. Molinari, Riv. Nuovo Cimento 14, 1 (1991). [34] J. E. Amaro, G. Co, and A. M. Lallena, Ann. Phys. (N.Y.) 221, 306 (1993). [35] J. E. Amaro, G. Co, and A. M. Lallena, Nucl. Phys. A578, 365 (1994). [36] A. Gil, J. Nieves, and E. Oset, Nucl. Phys. A627, 543 (1997). [37] M. J. Dekker, P. J. Brussaard, and J. A. Tjon, Phys. Lett. B 266, 249 (1991). [38] M. J. Dekker, P. J. Brussaard, and J. A. Tjon, Phys. Lett. B 289, 255 (1992). [39] M. J. Dekker, P. J. Brussaard, and J. A. Tjon, Phys. Rev. C 49, 2650 (1994). [40] A. De Pace, M. Nardi, W. M. Alberico, T. W. Donnelly, and A. Molinari, Nucl. Phys. A741, 249 (2004). [41] J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, and A. Molinari, Phys. Rep. 368, 317 (2002). [42] I. Korover et al. (The Jefferson Lab Hall A Collaboration), Phys. Rev. Lett. 113, 022501 (2014). RELATIVISTIC EFFECTS IN TWO-PARTICLE EMISSION …PHYSICAL REVIEW D 90, 033012 (2014) 033012-23