Implementation of chiral two-nucleon forces to nuclear many-body methods with Gaussian-wave packets

Abstract

Many-body methods that use Gaussian-wave packets to describe nucleon spatial distribution have been widely employed for depicting various phenomena in nuclear systems, in particular clustering. So far, however, the chiral effective field theory, a state-of-the-art theory of nuclear force, has not been applied to such methods. In this paper, we give the formalism to calculate the two-body matrix elements of the chiral two-nucleon forces using the Gaussian-wave packets. We also visualize the matrix elements and investigate the contributions of the central and tensor forces. This work is a foothold towards an ab initio description of various cluster phenomena in view of nucleons, pions, and many-nucleon forces.

1. Introduction

So far, an enormous number of studies have been carried out to understand various phenomena in nuclear systems from the viewpoint of nuclear force. For this purpose, many-body methods based on the local Gaussian basis functions have been established and widely applied; e.g. the fermionic molecular dynamics (FMD) [1,2] and antisymmetrized molecular dynamics (AMD) [3–5]. These methods employ the Gaussian-wave packet as the spatial part of the single-nucleon wave function, and have an advantage in describing localized cluster phenomena. However, a way to combine these methods with the chiral effective field theory (EFT) [6–8], a state-of-the-art theory of nuclear force, has not been proposed.

The chiral EFT derives nuclear forces from low-energy quantum chromodynamics (QCD), where nucleons and pions are the effective degree of freedom under the chiral symmetry of QCD. The chiral EFT succeeded phenomenologically, i.e. it gives realistic nuclear forces that describe nucleon–nucleon phase shifts with high precision (see, e.g. Ref. [9]). Furthermore, the chiral EFT has an advantage over previous theories in consistently handling many-nucleon forces. Therefore, nowadays, the potentials derived from the chiral EFT are regarded as a standard input of ab initio calculations.

The purpose of this work is to lay the foundation of implementing the chiral interaction to the many-body method based on the Gaussian-wave packets. We put the main focus on deriving the two-body matrix elements (MEs) of the chiral interaction by local Gaussians, and then, investigate their behavior. We expect that results of this work will be a milestone to understand how nuclear clusters emerge from the fundamental degree of freedom, namely, nucleons and pions. In this article, we adopt the two-nucleon force at next-to-next-to-next-to-leading order (N³LO). At this order, the chiral interaction becomes realistic, i.e. the χ²/datum for the reproduction of nucleon–nucleon scattering data reaches about unity [8,10–12]. The higher-order contributions and many-nucleon forces will be sought in forthcoming works.

Bare nucleon–nucleon interaction, such as the chiral interaction, leads to the strong tensor and short-range correlations. How to handle these correlations within the many-body methods with the Gaussian-wave packets has been intensively studied. However, this point lies outside the scope of this paper, since the main target is now the formulation of the two-body MEs, as mentioned above. Therefore, here, we just mention previous studies on the correlations with the Gaussian-wave packets.

As regards the tensor correlation, recent works with AMD [13,14] and the microscopic cluster model [15] showed that the application of the imaginary Gaussian center and spatially compact wave packets to these methods is an efficient way to take into account the high-momentum components of two-nucleon pairs induced by the tensor interaction. The formalism of this paper allows us to use complex numbers as the Gaussian center and vary the wave-packet size.

Also, there are many papers addressing how to manage the short-range correlation with the local Gaussian many-body methods. For instance, FMD was combined with the realistic Argonne interaction [16] transformed into a phase-shift equivalent low-momentum interaction by means of the unitary correlation operator method [17–19], and it was applied to the investigation of nuclear clusters [20,21] and the simulation of the astrophysical reactions, ³He(α, γ)⁷Be and ³H(α, γ)⁷Li [22]. Furthermore, the short-range correlation was investigated within the framework of AMD [23] and the quasi-cluster model [24] by directly using nucleon–nucleon potentials with a repulsive core. In future, the findings by these studies will be helpful to clarify the role of the short-range correlation induced by the chiral interaction within the local Gaussian many-body methods.

We put a comment on the first attempt to combine the cluster model and the chiral interaction. In our recent study [25], the effective interaction relevant to the Brink model [26] was obtained from the chiral interaction, but the noncentral forces were all missing due to presumed α clusters and a phenomenological prescription was introduced. Such shortcomings will be overcome by an ab initio calculation based on the formalism of this work.

This article is constructed as follows. In Sect. 2, we give the formalism of the two-body MEs in momentum space. In Sect. 3, as an example, typical MEs are visualized. Then, Sect. 4 is devoted to a summary and future perspectives. Appendix  A briefly shows how the MEs formulated in this work enter practical many-body calculations. The detail of the chiral potentials and information relevant to the formalism are relegated to Appendices  B and  C.

2. Formalism

2.1. Wave functions

We compute the two-body MEs in momentum space for the chiral interaction as the two-nucleon force, while other terms of the Hamiltonian can be calculated analytically (see, e.g. Ref. [26]). Since the main focus of this paper is the derivation of the two-body MEs of the chiral interaction by local Gaussians, we express the single-particle wave function in momentum space:

$$\begin{eqnarray} \left|\varphi _i^{(\nu )}\right\rangle &= \left|\phi _i^{(\nu )}\chi _i\right\rangle , \end{eqnarray}$$

(1)

$$\begin{eqnarray} \left\langle {\boldsymbol{p}_n}\,\vert \,{\phi _i^{(\nu )}}\right\rangle = \phi _i^{(\nu )}(\boldsymbol{p}_n) = \frac{1}{\left(2\pi \nu \right)^{\frac{3}{4}}} \exp \!\left[-\frac{1}{4\nu }(\boldsymbol{p}_n+2i\nu \boldsymbol{Z}_i)^2 \right]\exp \!\left[-\nu Z_i^2\right], \end{eqnarray}$$

(2)

$$\begin{eqnarray} \left|\chi _i\right\rangle &= \left[\alpha _i\left|\uparrow _{\sigma }\right\rangle +\beta _i\left|\downarrow _{\sigma }\right\rangle \right] \left|m_{\tau _i}\right\rangle , \end{eqnarray}$$

(3)

where |$\phi _i^{(\nu )}$| is the Fourier transform of its coordinate-space representation:

$$\begin{eqnarray} \phi _i^{(\nu )}(\boldsymbol{p}_n) = \frac{1}{(2\pi )^{\frac{3}{2}}} \left(\frac{2\nu }{\pi }\right)^{\frac{3}{4}} \int d\boldsymbol{r} \exp \!\left(-i\boldsymbol{p}_n\cdot \boldsymbol{r}\right) \exp \!\left[-\nu \left(\boldsymbol{r}-\boldsymbol{Z}_i\right)^2\right], \end{eqnarray}$$

(4)

with |$\boldsymbol{p}_n$| the momentum of the nth nucleon. In this paper we use natural units such that ℏ = c = 1. The Gaussian-wave packet is characterized by the range parameter ν and the Gaussian center |$\boldsymbol{Z}_i$|⁠, which is in general complex. The spin-isospin state |χ_i〉 is expressed as a superposition of the spin-up and spin-down states, |↑_σ〉 and |↓_σ〉, respectively. The weights α_i and β_i are determined by many-body calculations through variational processes. The isospin state can be |m_τ〉 = |↑_τ〉 (neutron) or |m_τ〉 = |↓_τ〉 (proton).

As argued in Refs. [13,15], the imaginary part of |$\boldsymbol{Z}_i$| corresponds to the expectation value of the nucleon momentum. Indeed, if |$\boldsymbol{Z}_i$| is pure imaginary in Eq. (2), |$\phi _i^{(\nu )}$| localizes at around |$\boldsymbol{p}_n\sim 2\nu \mathrm{Im}(\boldsymbol{Z}_i)$|⁠. Therefore, complex |$\boldsymbol{Z}_i$| enables us to efficiently take into account the high-momentum components caused by the tensor correlation. Also, inclusion of the spatially compact wave packets by varying ν can improve calculations [14]. Our formalism allows us to vary ν and we explain it in Appendix  A, where many-body states are defined using the single-particle wave function of Eq. (1).

From Eq. (1), the two-body state in momentum space reads

$$\begin{eqnarray} \left\langle { \boldsymbol{p}_1, \boldsymbol{p}_2 }\,\vert \,{\phi _i^{(\nu )}\phi _j^{(\nu )}}\right\rangle = \frac{1}{\left(2\pi \nu \right)^{\frac{3}{2}}} \exp \!\left[-\frac{1}{4\nu }\left\lbrace (\boldsymbol{p}_1+2i\nu \boldsymbol{Z}_i)^2+(\boldsymbol{p}_2+2i\nu \boldsymbol{Z}_j)^2\right\rbrace -\nu \left(Z_i^2+Z_j^2\right)\right], \end{eqnarray}$$

(5)

which can be rewritten in terms of the center-of-mass (CM) momentum |$\boldsymbol{P}$| and the relative momentum |$\boldsymbol{p}$| as

$$\begin{eqnarray} \left\langle { \boldsymbol{p}_1, \boldsymbol{p}_2 }\,\vert \,{\phi _i^{(\nu )}\phi _j^{(\nu )}}\right\rangle = \frac{1}{\left(2\pi \nu \right)^{\frac{3}{2}}} \exp \!\left[-\frac{P^2}{8\nu }-i\boldsymbol{Z}_{ij}\cdot \boldsymbol{P}\right] \exp \!\left[-\frac{p^2}{2\nu }-i\boldsymbol{z}_{ij}\cdot \boldsymbol{p}\right], \end{eqnarray}$$

(6)

with

$$\begin{eqnarray} {\begin{pmatrix}\boldsymbol{Z}_{ij} \\ \boldsymbol{z}_{ij} \end{pmatrix}} = \mathcal {U} {\begin{pmatrix}\boldsymbol{Z}_i \\ \boldsymbol{Z}_j \end{pmatrix}}, \qquad {\begin{pmatrix}\boldsymbol{P} \\ \boldsymbol{p} \end{pmatrix}} = \left(\mathcal {U}^{-1}\right)^T {\begin{pmatrix}\boldsymbol{p}_1 \\ \boldsymbol{p}_2 \end{pmatrix}}, \qquad \mathcal {U} = {\begin{pmatrix}\frac{1}{2} && \,\,\,\, \frac{1}{2} \\ 1 && -1 \end{pmatrix}}. \end{eqnarray}$$

(7)

Here, neutrons and protons are assumed to have common mass.

2.2. General form of two-body MEs

In this section we derive the general expression of the two-body MEs in momentum space with a nonlocal potential dependent on the relative momenta, |$\boldsymbol{p}$| and |$\boldsymbol{p}^{\prime }$|⁠, of the initial and final channels, respectively. The nonlocality originates from the regularization scheme of the chiral EFT (see Appendix  B).

For simplicity, the nucleon spin is assumed to be either |↑_σ〉 or |↓_σ〉. Note that this does not undermine the formalism by loss of generality . Thus, we derive the MEs based on the single-particle state |$\left|\phi _i^{(\nu )}\xi _i\right\rangle$|⁠, where |$\phi _i^{(\nu )}$| is defined by Eq. (2) and |ξ_i〉 is now given by

$$\begin{eqnarray} \left|\xi _i\right\rangle &= \left|m_{\sigma _i}\right\rangle \left|m_{\tau _i}\right\rangle . \end{eqnarray}$$

(8)

Hence, the spin-up and spin-down states are respectively expressed by |$\left|m_{\sigma _i}\right\rangle =\left|\uparrow _{\sigma }\right\rangle$| and |$\left|m_{\sigma _i}\right\rangle =\left|\downarrow _{\sigma }\right\rangle$|⁠. Now the two-body states are represented by

$$\begin{eqnarray} \left|i j,\nu \right\rangle =\left|\phi _i^{(\nu )}\phi _j^{(\nu )}\xi _i\xi _j\right\rangle , \end{eqnarray}$$

(9)

with which the MEs are formulated.

Let us move to the general form of the two-body MEs. A general two-body operator |$\hat{V}_{2N}$| can be written as

$$\begin{eqnarray} \hat{V}_{2N} = \int \!\!\!\!\int \!\!\!\!\int \!\!\!\!\int d\boldsymbol{p} d\boldsymbol{p}^{\prime } d\boldsymbol{P} d\boldsymbol{P}^{\prime}\left.|\boldsymbol{p}^{\prime },\boldsymbol{P}^{\prime }\right\rangle \left \langle \boldsymbol{p}^{\prime }, \boldsymbol{P}^{\prime }\left | \hat{V}_{2N} \right | \boldsymbol{p}, \boldsymbol{P} \right \rangle \left\langle \boldsymbol{p},\boldsymbol{P} \right. | . \end{eqnarray}$$

(10)

The Galilean invariance requires that the potential v_2N in association with the operator |$\hat{V}_{2N}$| does not depend on |$\boldsymbol{P}$| and |$\boldsymbol{P}^{\prime }$|⁠. Hence, such a potential is given by

$$\begin{eqnarray} \left\langle \boldsymbol{p}^{\prime }, \boldsymbol{P}^{\prime } \left| \left. \hat{V}_{2N}\right. \right| \boldsymbol{p}, \boldsymbol{P} \right\rangle = v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) \delta (\boldsymbol{P}-\boldsymbol{P}^{\prime }), \end{eqnarray}$$

(11)

which results in the two-body MEs:

$$\begin{eqnarray} && \left\langle kl,\nu ^{\prime } \left| \left. \hat{V}_{2N}\right. \right| ij,\nu \right\rangle \\ && \quad = \int \!\!\!\!\int \!\!\!\!\int d\boldsymbol{p} d\boldsymbol{p}^{\prime } d\boldsymbol{P} \left\langle kl,\nu ^{\prime } \left| \boldsymbol{p}^{\prime }, \boldsymbol{P} \right\rangle \right. v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) \left\langle \boldsymbol{p}, \boldsymbol{P} \left| ij,\nu \right\rangle \right. \\ && \quad = A_{ijkl}^{(\nu \nu ^{\prime })} \int \!\!\!\!\int d\boldsymbol{p} d\boldsymbol{p}^{\prime } \left\langle \xi _k\xi _l \left| \left. v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) \right. \right| \xi _i\xi _j \right\rangle g_{\nu \nu ^{\prime }}(p,p^{\prime }) \exp \!\left[-i\boldsymbol{z}_{ij}\cdot \boldsymbol{p}+i\boldsymbol{z}_{kl}^{*}\cdot \boldsymbol{p}^{\prime }\right], \end{eqnarray}$$

(12)

with

$$\begin{eqnarray} g_{\nu \nu ^{\prime }}(p,p^{\prime }) = \exp \!\left[-\frac{p^2}{2\nu }-\frac{p^{\prime 2}}{2\nu ^{\prime }}\right], \end{eqnarray}$$

(13)

$$\begin{eqnarray} A_{ijkl}^{(\nu \nu ^{\prime })} = \left[\frac{2}{\pi (\nu +\nu ^{\prime })}\right]^{\!\frac{3}{2}} \exp \!\left[-\frac{2\nu \nu ^{\prime }}{\nu +\nu ^{\prime }}\left(\boldsymbol{Z}_{ij}-\boldsymbol{Z}_{kl}^{*}\right)^2\right]. \end{eqnarray}$$

(14)

Using Eq. (12), we introduce the notation of the antisymmetrized two-body MEs,

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime })}= \left\langle kl,\nu ^{\prime } \left| \left. \hat{V}_{2N}\right. \right| ij,\nu \right\rangle -\left\langle kl,\nu ^{\prime } \left| \left. \hat{V}_{2N}\right. \right| ji,\nu \right\rangle . \end{eqnarray}$$

(15)

For later convenience, we perform the Rayleigh expansion for the plane waves in Eq. (12):

$$\begin{eqnarray} \exp \!\left[-i\boldsymbol{z}_{ij}\cdot \boldsymbol{p}\right] = 4\pi \sum _{L} (-)^L i^{L} j_{L}\!\left(z_{ij}p\right) \sum _{M}(-)^M Y_{LM}\!\left(\hat{\boldsymbol{p}}\right) Y_{L,-M}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right), \end{eqnarray}$$

(16)

$$\begin{eqnarray} \exp \!\left[i\boldsymbol{z}_{kl}^{*}\cdot \boldsymbol{p}^{\prime }\right] = 4\pi \sum _{L^{\prime }} i^{L^{\prime }} j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right) \sum _{M^{\prime }} Y_{L^{\prime }M^{\prime }}^{*}\!\left(\hat{\boldsymbol{p}}^{\prime }\right) Y_{L^{\prime }M^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right), \end{eqnarray}$$

(17)

where |$\boldsymbol{z}_{ij}=(z_{ij,x},z_{ij,y},z_{ij,z})$|⁠. Note that z_{ij, x}, z_{ij, y}, and z_{ij, z} are complex numbers. Hence, the arguments of the spherical Bessel function j_L are also complex numbers as |${z_{ij}=\sqrt{\boldsymbol{z}_{ij}\cdot \boldsymbol{z}_{ij}}}$| and |$z_{kl}^{*}=\sqrt{\boldsymbol{z}_{kl}^{*}\cdot \boldsymbol{z}_{kl}^{*}}$|⁠. The spherical harmonics Y_LM with the complex numbers as its arguments is given by the analytic continuation in the Cartesian representation:

$$\begin{eqnarray} Y_{lm}(z_x,z_y,z_z) = &&(-)^m \sum _{\lambda =0}^{[(l-m)/2]} \left[\frac{(2l+1)(l-m)!}{4\pi (l+m)!}\right]^{\frac{1}{2}} \frac{(l+m)!}{2^m} \frac{(-)^\lambda }{4^\lambda \lambda !(m+\lambda )!(l-m-2\lambda )!} \\ &&\times \left(\frac{z_z}{z}\right)^{l-m-2\lambda } \left(\frac{z_x+iz_y}{z}\right)^{m+\lambda } \left(\frac{z_x-iz_y}{z}\right)^\lambda \qquad (m\ge 0). \end{eqnarray}$$

(18)

The upper limit of the summation is expressed by the floor function and |${z=\sqrt{\boldsymbol{z}\cdot \boldsymbol{z}}=\sqrt{z_x^2+z_y^2+z_z^2}}$| is a complex number.

2.3. Two-body MEs of chiral interaction

2.3.1. Overview

Now we formulate the two-body MEs of the chiral interaction at N³LO. One finds that the chiral-N³LO potential v_2N consists of the central, spin-orbit (SO), and tensor terms. Therefore, the MEs can also be decomposed into these terms as |$V_{ijkl}^{(\nu \nu ^{\prime })}=V_{ijkl}^{(\nu \nu ^{\prime }:C)}+V_{ijkl}^{(\nu \nu ^{\prime }:LS)}+V_{ijkl}^{(\nu \nu ^{\prime }:T)}+V_{ijkl}^{(\nu \nu ^{\prime }:\sigma L)}$|⁠, where the central, SO, tensor, and σL contributions are denoted by |$V_{ijkl}^{(\nu \nu ^{\prime }:C)}$|⁠, |$V_{ijkl}^{(\nu \nu ^{\prime }:LS)}$|⁠, |$V_{ijkl}^{(\nu \nu ^{\prime }:T)}$|⁠, and |$V_{ijkl}^{(\nu \nu ^{\prime }:\sigma L)}$|⁠, respectively. Although the σL term, the potential of which depends on |$V_{\sigma L}^{(4)}$| and D₁₅ in Eqs. (B27) and (B28), respectively, is one of the tensor contributions, for simplicity, we treat it separately from standard tensor terms.

In the following sections, first we derive the typical MEs of the central, SO, and tensor contributions to explain basic procedures of the calculations. The MEs of the contact terms at next-to-leading order (NLO) are chosen as examples since they consist of the central, SO, and tensor components, and their MEs are relatively simple due to the absence of the pion propagators. Then, unified expressions of the MEs are given in Sect. 2.3.5.

2.3.2. Central contributions

As an example of the central contributions, we take the C₁ term, one of the contact terms at NLO. As shown by Eq. (B9), with nonlocal regularization, its potential is given by

$$\begin{eqnarray} v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) = \frac{1}{(2\pi )^3}C_1 q^2 u_n\!\left(p,p^{\prime },\Lambda \right), \end{eqnarray}$$

(19)

where C₁ is one of the low-energy constants (LECs). The transferred momentum |$\boldsymbol{q}$|⁠, as well as the average momentum |$\boldsymbol{Q}$|⁠, are defined by

$$\begin{eqnarray} {\begin{pmatrix}\boldsymbol{Q} \\ \boldsymbol{q} \end{pmatrix}} = \mathcal {U} {\begin{pmatrix}\boldsymbol{p}^{\prime } \\ \boldsymbol{p} \end{pmatrix}}, \end{eqnarray}$$

(20)

with the matrix |$\mathcal {U}$| given by Eq. (7). The regulator u_n characterized by the power n and the cutoff momentum Λ is nonlocal:

$$\begin{eqnarray} u_n\!\left(p,p^{\prime },\Lambda \right) = \exp \!\left[-\!\left(\frac{p}{\Lambda }\right)^{\!2n}-\!\left(\frac{p^{\prime }}{\Lambda }\right)^{\!2n}\right]. \end{eqnarray}$$

(21)

See Appendix  B for more detail of the potential form. To compute the MEs, we perform the multipole expansion (MPE) of q² in Eq. (19):

$$\begin{eqnarray} q^2 = 4\pi \sum _{K=0,1} \frac{(-)^K}{\hat{K}}\mathcal {F}^{(C_1)}_K\!\left(p,p^{\prime }\right) \left[Y_K\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_K\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{00}, \end{eqnarray}$$

(22)

and the MPE function is given as

$$\begin{eqnarray} \mathcal {F}^{(C_1)}_K\!\left(p,p^{\prime }\right) = \left(p^2+p^{\prime 2}\right)\delta _{K 0} -2pp^{\prime }\delta _{K 1}. \end{eqnarray}$$

(23)

The bipolar spherical harmonics is defined by

$$\begin{eqnarray} \left[Y_L\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{L^{\prime }}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{\lambda \mu } = \sum _{M=-L}^{L} \sum _{M^{\prime }=-L^{\prime }}^{L^{\prime }} \left( L M L^{\prime } M^{\prime } | \lambda \mu \right) Y_{LM}\!\left(\hat{\boldsymbol{p}}\right) Y_{L^{\prime }M^{\prime }}\!\left(\hat{\boldsymbol{p}}^{\prime }\right). \end{eqnarray}$$

(24)

We adopt the abbreviation |$\hat{L}=\sqrt{2L+1}$| and the Clebsch–Gordan coefficient is denoted by (⋅⋅⋅ · · · | · ·). In Appendix  C, the detail of the MPE can be found.

Using these expressions and plugging Eq. (19) into Eq. (12), as well as employing Eqs. (16) and (17), one finds the antisymmetrized MEs as

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:C)} =&& \frac{2}{\pi } C_1 A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{L=0,1} \left[D_{ijkl}^{(\sigma )} D_{ijkl}^{(\tau )}-(-)^LD_{jikl}^{(\sigma )} D_{jikl}^{(\tau )}\right] P_L(z_{ijkl}) \\ && \quad \times \int \!\!\!\!\int dp dp^{\prime } p^{2}p^{\prime 2} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L}\!\left(z_{kl}^{*}p^{\prime }\right) \mathcal {F}^{(C_1)}_L\!\left(p,p^{\prime }\right). \end{eqnarray}$$

(25)

Here, the Legendre polynomial P_L has the argument given by

$$\begin{eqnarray} z_{ijkl} = \left\lbrace \begin{aligned} &1 && \qquad (\text{$i=j$ and/or $k=l$}),\\ &\frac{\boldsymbol{z}_{ij}\cdot \boldsymbol{z}_{kl}^{*}}{z_{ij}z_{kl}^{*}} && \qquad (\text{$i\ne j$ and $k\ne l$}), \end{aligned} \right. \end{eqnarray}$$

(26)

and |$D_{ijkl}^{(\sigma )}$| is defined by

$$\begin{eqnarray} D_{ijkl}^{(\sigma )} &= \delta _{m_{\sigma _i} m_{\sigma _k}} \delta _{m_{\sigma _j} m_{\sigma _l}}. \end{eqnarray}$$

(27)

The superscripts (σ) and (τ) stand for the spin and isospin MEs, respectively.

We emphasize that any central term of the chiral potential can be represented in terms of |$\left[Y_K\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_K\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{00}$| as per Eq. (22). This is because the spatial and spin parts of the central operators are decoupled, i.e. they both form a scalar operator individually. For the pion-exchange terms, the MPE of the q- and Q-dependent parts of the potentials involving the pion propagators is necessary, but after some manipulation, |$\left[Y_K\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_K\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{00}$| must appear. As a result, the structure of the central MEs is always expressed by Eq. (25). Indeed, for other central terms, we just have to replace the prefactor 2C₁/π and |$\mathcal {F}^{(C_1)}_L$| with appropriate forms, as explained in Sect. 2.3.5 and Appendix  C. Note that, in principle, the summation over L runs up to infinity for general cases, and the spin-isospin MEs need to be modified depending on the operator forms (see Table 1).

Table 1.

The spin-isospin MEs |$U_{ijkl,L}^{(X)}$| of the chiral interaction. The tensor operator can be |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)$|⁠, |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{Q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{Q}\right)$|⁠, and |$\left[\boldsymbol{\sigma }_1\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \left[\boldsymbol{\sigma }_2\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right]$|⁠. In the rightmost column, the corresponding terms are listed. The N³LO potentials are characterized by the superscripts, for which c_i, m_N, and 2L denote the LEC, the average nucleon mass, and the two-loop contributions, respectively (see Appendix B.5).

Operator type	|$U_{ijkl,L}^{(X)}$|	Chiral interaction
|$\mathbb {1}$|	|$D_{ijkl}^{(\sigma )} D_{ijkl}^{(\tau )} -(-)^LD_{jikl}^{(\sigma )} D_{jikl}^{(\tau )}$|	CS, C1, C2, |$V_C^{(3)}$|
		D1, D2, D3, D4
		|$V_C^{(c_i^2)}$|⁠, |$V_C^{(c_i/m_N)}$|⁠, |$V_C^{(m_N^{-2})}$|⁠, |$V_C^{(\mathrm{2L})}$|
|$\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2$|	|$X_{ijkl}^{(\sigma )} D_{ijkl}^{(\tau )} -(-)^LX_{jikl}^{(\sigma )} D_{jikl}^{(\tau )}$|	CT, C3, C4, |$V_S^{(2)}$|⁠, |$V_S^{(3)}$|
		D5, D6, D7, D8, |$V_S^{(m_N^{-2})}$|⁠, |$V_S^{(\mathrm{2L})}$|
|$\boldsymbol{\tau }_1\cdot \boldsymbol{\tau }_2$|	|$D_{ijkl}^{(\sigma )} X_{ijkl}^{(\tau )} -(-)^LD_{jikl}^{(\sigma )} X_{jikl}^{(\tau )}$|	|$W_C^{(2)}$|⁠, |$W_C^{(3)}$|
		|$W_C^{(c_i/m_N)}$|⁠, |$W_C^{(m_N^{-2})}$|⁠, |$W_C^{(\mathrm{2L})}$|
|$(\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2)(\boldsymbol{\tau }_2\cdot \boldsymbol{\tau }_2)$|	|$X_{ijkl}^{(\sigma )} X_{ijkl}^{(\tau )} -(-)^LX_{jikl}^{(\sigma )} X_{jikl}^{(\tau )}$|	|$W_S^{(3)}$|
		|$W_S^{(c_i^2)}$|⁠, |$W_S^{(c_i/m_N)}$|⁠, |$W_S^{(m_N^{-2})}$|⁠, |$W_S^{(\mathrm{2L})}$|
|$-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)$|	|$S_{ijkl} \left[D_{ijkl}^{(\tau )}-(-)^LD_{ijkl}^{(\tau )}\right]$|	C5, |$V_{LS}^{(3)}$|
		D9, D10, |$V_{LS}^{(c_i/m_N)}$|⁠, |$V_{LS}^{(m_N^{-2})}$|
|$-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)$|	|$S_{ijkl} \left[X_{ijkl}^{(\tau )}-(-)^LX_{jikl}^{(\tau )}\right]$|	|$W_{LS}^{(3)}$|⁠, |$W_{LS}^{(c_i/m_N)}$|⁠, |$W_{LS}^{(m_N^{-2})}$|
|$\times \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2$|
Tensor	|$T_{ijkl}^{(\lambda _0)}D_{ijkl}^{(\tau )}-(-)^LT_{jikl}^{(\lambda _0)}D_{jikl}^{(\tau )}$|	C6, C7, |$V_T^{(2)}$|⁠, |$V_T^{(3)}$|
		D11, D12, D13, D14, D15
		|$V_T^{(m_N^{-2})}$|⁠, |$V_{\sigma L}^{(m_N^{-2})}$|⁠, |$V_T^{(\mathrm{2L})}$|
(Tensor)|$\boldsymbol{\tau }_1\cdot \boldsymbol{\tau }_2$|	|$T_{ijkl}^{(\lambda _0)}X_{ijkl}^{(\tau )}-(-)^LT_{jikl}^{(\lambda _0)}X_{jikl}^{(\tau )}$|	1π, |$W_T^{(3)}$|
		|$W_T^{(c_i^2)}$|⁠, |$W_T^{(c_i/m_N)}$|⁠, |$W_T^{(m_N^{-2})}$|⁠, |$W_T^{(\mathrm{2L})}$|

2.3.3. SO contributions

The C₅ term, one of the contact terms at NLO, is a good example of the SO contributions. Its potential reads

$$\begin{eqnarray} v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) = \frac{1}{(2\pi )^3}C_5 u_n\!\left(p,p^{\prime },\Lambda \right) \left[-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right], \end{eqnarray}$$

(28)

where C₅ is the LEC. The total spin |$\boldsymbol{S}$| is given by

$$\begin{eqnarray} \boldsymbol{S}&=\frac{1}{2}\left(\boldsymbol{\sigma }_1+\boldsymbol{\sigma }_2\right). \end{eqnarray}$$

(29)

The spin operator |$\boldsymbol{\sigma }_i$| is represented by the Pauli matrices. The operator |$-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)$| can be represented in terms of |$\boldsymbol{p}$| and |$\boldsymbol{p}^{\prime }$| as

$$\begin{eqnarray} -i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right) =&& 4\pi \sqrt{\frac{3}{2}} \sum _{\lambda _q=0,1} \sum _{\lambda _Q=0,1} (-)^{\lambda _q}\widehat{1\!-\!\lambda _q}\widehat{1\!-\!\lambda _Q} \left[\binom{3}{2\lambda _q}\binom{3}{2\lambda _Q}\right]^{\!\frac{1}{2}} p^{\lambda _q+\lambda _Q} p^{\prime \, 2-\lambda _q-\lambda _Q} \\ &&\times \sum _{K_1K_2} \left( \lambda _q 0 \lambda _Q 0 | K_1 0 \right) \left( 1\!-\!\lambda _q, 0, 1\!-\!\lambda _Q, 0 | K_2 0 \right) {\left\lbrace \begin{array}{*{10}c}\lambda _q & \quad 1\!-\!\lambda _q & \quad 1 \\ \lambda _Q & \quad 1\!-\!\lambda _Q & \quad 1 \\ K_1 & \quad K_2 & \quad 1 \end{array}\right\rbrace } \\ &&\times \left[ \left[Y_{K_1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{K_2}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{1} \otimes S_1 \right]_{00}, \end{eqnarray}$$

(30)

where the binomial coefficient is given by

$$\begin{eqnarray} \binom{2n+1}{2m} = \frac{\left(2n+1\right)!}{\left(2n+1-2m\right)!\left(2m\right)!}, \end{eqnarray}$$

(31)

and the 9-j symbol is represented by the 3 × 3 matrix in the braces.

From Eqs. (12), (28), and (30), the antisymmetrized MEs can be calculated as

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:LS)} &=& 4\sqrt{3} C_5 A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{\lambda _q\lambda _Q} \sum _{L=0}^\infty \sum _{L^{\prime }=0}^\infty (-)^{\lambda _Q} i^{L+L^{\prime }} \widehat{1\!-\!\lambda _q}\widehat{1\!-\!\lambda _Q} \left[\binom{3}{2\lambda _q}\binom{3}{2\lambda _Q}\right]^{\!\frac{1}{2}} \\ && \times S_{ijkl}\left[D_{ijkl}^{(\tau )}-(-)^LD_{jikl}^{(\tau )}\right] \\ && \times \left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{1,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ &&\times \left( \lambda _q 0 \lambda _Q 0 | L 0 \right) \left( 1\!-\!\lambda _q, 0, 1\!-\!\lambda _Q, 0 | L^{\prime } 0 \right) {\left\lbrace \begin{array}{*{10}c}\lambda _q & \quad 1\!-\!\lambda _q & \quad 1 \\ \lambda _Q & \quad 1\!-\!\lambda _Q & \quad 1 \\ L & \quad L^{\prime } & \quad 1 \end{array}\right\rbrace } \\ && \times \int \!\!\!\!\int dp dp^{\prime } p^{\lambda _q+\lambda _Q+2}p^{\prime 4-\lambda _q-\lambda _Q} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right), \end{eqnarray}$$

(32)

with

$$\begin{eqnarray} S_{ijkl} =&& \left(\left. \frac{1}{2} m_{\sigma _i} \frac{1}{2} m_{\sigma _j} \right| 1, m_{\sigma _i}\!+m_{\sigma _j} \right) \left(\left. \frac{1}{2} m_{\sigma _k} \frac{1}{2} m_{\sigma _l} \right| 1, m_{\sigma _k}\!+m_{\sigma _l} \right) \\ &&\times \left( 1, m_{\sigma _i}\!+m_{\sigma _j}\!-m_{\sigma _k}\!-m_{\sigma _l}, 1, m_{\sigma _k}\!+m_{\sigma _l} | 1, m_{\sigma _i}\!+m_{\sigma _j} \right). \end{eqnarray}$$

(33)

The expression, Eq. (30), can be applicable to other SO terms, i.e. all SO potentials can be written in terms of |$\left[ \left[Y_{K_1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{K_2}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{1} \otimes S_1 \right]_{00}$| even if pions are exchanged. Consequently, for the MEs of the SO contributions, the structure of Eq. (32) is rather general as shown in Sect. 2.3.5. In general cases, several prefactors and angular-momentum-coupling coefficients in Eq. (32) are packed into a single function, which is the MPE function of the SO term described in Appendix C.2.

Note that Eq. (32) can be further simplified since only L = L^′ = 1 is allowed and other configurations are forbidden by the angular-momentum-coupling coefficients involved. However, we do not show explicitly such simplified MEs because this section is intended to demonstrate the derivation of the SO MEs and the general form, Eq. (32), is useful rather than showing such specific MEs.

2.3.4. Tensor contributions

The typical tensor contribution appears as the C₆-contact term at NLO. As shown in Appendix B.3, its potential is written as

$$\begin{eqnarray} v_{2N}\!\left(\boldsymbol{p}^{\prime },\boldsymbol{p}\right) = \frac{1}{(2\pi )^3}C_6u_n\!\left(p,p^{\prime },\Lambda \right) \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right), \end{eqnarray}$$

(34)

where C₆ is the LEC. We express the operator in Eq. (34) by the irreducible-tensor representation:

$$\begin{eqnarray} \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) &=& 4\pi \sum _{\lambda _0=0,2}\sum _{\lambda _1=0}^{\lambda _0} \frac{(-)^{\lambda _1}}{\hat{\lambda }_1} \binom{2\lambda _0+1}{2\lambda _1}^{\!\frac{1}{2}} \left( 1 0 1 0 | \lambda _0 0 \right) q^{2-\lambda _0} p^{\lambda _1}p^{\prime \lambda _0-\lambda _1} \\ && \times \left[ \left[Y_{\lambda _1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{\lambda _0-\lambda _1}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{\lambda _0} \otimes \left[\sigma _1(1)\otimes \sigma _1(2)\right]_{\lambda _0} \right]_{00} \\ & =& 4\pi \sum _{\lambda _0\lambda _1}\sum _{K=0,1} \binom{2\lambda _0+1}{2\lambda _1}^{\!\frac{1}{2}} \widehat{\lambda _0\!-\!\lambda _1} \left( 1 0 1 0 | \lambda _0 0 \right) \mathcal {F}_{\lambda _0 K}^{(C_6)}\!\left(p,p^{\prime }\right) p^{\lambda _1}p^{\prime \lambda _0-\lambda _1} \\ && \times \sum _{K_1K_2}(-)^{K_2} \left( \lambda _1 0 K 0 | K_1 0\right) \left( \lambda _0\!-\!\lambda _1, 0 K 0 | K_2 0\right) {\left\lbrace \begin{array}{*{10}c}\lambda _0\!-\!\lambda _1 & \quad \lambda _1 & \quad \lambda _0 \\ K_1 & \quad K_2 & \quad K \end{array}\right\rbrace } \\ && \times \left[ \left[Y_{K_1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{K_2}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{\lambda _0} \otimes \left[\sigma _{1}(1)\otimes \sigma _{1}(2)\right]_{\lambda _0} \right]_{00}, \end{eqnarray}$$

(35)

with

$$\begin{eqnarray} \mathcal {F}_{\lambda _0 K}^{(C_6)}\!\left(p,p^{\prime }\right) = \left[ \left(p^2+p^{\prime 2}\right)\delta _{K 0}-2pp^{\prime }\delta _{K 1} \right]\delta _{\lambda _0 0} + \delta _{\lambda _0 2}\delta _{K 0}. \end{eqnarray}$$

(36)

The 2 × 3 matrix in the braces is the 6-j symbol.

One finds that Eq. (35) with λ₀ = 0 corresponds to the central component, |$q^2 \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2\right)/3$|⁠, whereas that with λ₀ = 2 is the purely tensor component, |$q^2\mathcal {S}_{12}/3$|⁠, with the tensor operator

$$\begin{eqnarray} \mathcal {S}_{12}\!\left(\hat{\boldsymbol{q}}\right) &=& 3\left(\boldsymbol{\sigma }_1\cdot \hat{\boldsymbol{q}}\right) \left(\boldsymbol{\sigma }_2\cdot \hat{\boldsymbol{q}}\right) -\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2 \\ &=& 2\left[3\left(\boldsymbol{S}\cdot \hat{\boldsymbol{q}}\right)^2 -\boldsymbol{S}^2\right]. \end{eqnarray}$$

(37)

Note that |$\boldsymbol{S}$| is defined by Eq. (29). In this paper, the terms dependent on |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)$| are referred to as the tensor contributions, although they involve the central contributions.

From the above expressions, the antisymmetrized MEs are computed as

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:T)} &=& 24 C_6 A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{\lambda _0\lambda _1} \binom{2\lambda _0+1}{2\lambda _1}^{\!\frac{1}{2}} \frac{\widehat{\lambda _0\!-\!\lambda _1}}{\hat{\lambda }_0}\left( 1 0 1 0 | \lambda _0 0 \right) \\ && \times \sum _{LL^{\prime }K} i^{L+L^{\prime }} \left[T_{ijkl}^{(\lambda _0)} D_{ijkl}^{(\tau )}-(-)^LT_{jikl}^{(\lambda _0)} D_{jikl}^{(\tau )}\right] \\ && \times \left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{\lambda _0,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ && \times \left( \lambda _1 0 K 0 | L 0\right) \left( \lambda _0\!-\!\lambda _1, 0 K 0 | L^{\prime } 0\right) {\left\lbrace \begin{array}{*{10}c}\lambda _0\!-\!\lambda _1 && \lambda _1 && \lambda _0 \\ L && L^{\prime } && K \end{array}\right\rbrace } \\ && \times \int \!\!\!\!\int dp dp^{\prime } p^{\lambda _1+2}p^{\prime \lambda _0-\lambda _1+2} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right) \mathcal {F}_{\lambda _0 K}^{(C_6)}\!\left(p,p^{\prime }\right), \\ \end{eqnarray}$$

(38)

with

$$\begin{eqnarray} T_{ijkl}^{(\lambda _0)} =&& (-)^{m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \left( 1, m_{\sigma _i}\!-m_{\sigma _k}, 1, m_{\sigma _j}\!-m_{\sigma _l} | \lambda _0, m_{\sigma _i}\!+m_{\sigma _j}\!-m_{\sigma _k}\!-m_{\sigma _l} \right) \\ &&\times \left(\left. \frac{1}{2} m_{\sigma _i} 1, m_{\sigma _k}\!-m_{\sigma _i} \right| \frac{1}{2} m_{\sigma _k} \right) \left(\left. \frac{1}{2} m_{\sigma _j} 1, m_{\sigma _l}\!-m_{\sigma _j} \right| \frac{1}{2} m_{\sigma _l} \right). \end{eqnarray}$$

(39)

There is no counterpart of |$T_{ijkl}^{(\lambda _0)}$| for the isospin indices, and therefore, the superscript (σ) is not necessary in Eq. (39).

The point of the calculations of the tensor MEs is that the irreducible-tensor representation by Eq. (35) is valid also for other tensor terms, no matter whether they have the pion propagators or the operator is given by |$\boldsymbol{Q}$| instead of |$\boldsymbol{q}$|⁠. Every tensor potential can be expressed in terms of |$\big [\! \left[Y_{K_1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{K_2}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{\lambda _0} \otimes \left[\sigma _{1}(1)\otimes \sigma _{1}(2)\right]_{\lambda _0} \!\big ]_{00}$| with λ₀ = 0 or 2. Therefore, again, the structure of Eq. (38) is essentially the same as that of other tensor terms, and differences can be found only in the prefactors and the MPE function |$\mathcal {F}_{\lambda _0 K}^{(C_6)}$|⁠. A general form of the tensor MEs is thus obtained in Sect. 2.3.5.

Note that, following the procedure for the C₆-term MEs, one can derive the MEs of the σL terms in association with the operator, |$\left[\boldsymbol{\sigma }_1\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \left[\boldsymbol{\sigma }_2\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right]$|⁠, which appears in Eqs. (B27) and (B28). Indeed, this operator can be written in terms of |$\big [\! \left[Y_{K_1}\!\left(\hat{\boldsymbol{p}}\right)\otimes Y_{K_2}\!\left(\hat{\boldsymbol{p}}^{\prime }\right)\right]_{\lambda _0} \otimes \left[\sigma _{1}(1)\otimes \sigma _{1}(2)\right]_{\lambda _0} \!\big ]_{00}$|⁠. As a result, their MEs have a structure essentially the same as that of Eq. (38), i.e. λ₀ = 0 and 2 correspond to the central and tensor components respectively, although the prefactors and coefficients are much more complicated.

2.3.5. Summary of chiral two-body MEs

In the previous sections, the typical MEs of the central, SO, and tensor contributions are shown. Now, in this section, we generalize them:

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:C)} &=& A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{L=0}^\infty U_{ijkl,L}^{(C)} P_L(z_{ijkl}) \\ &&\times \int \!\!\!\!\int dp dp^{\prime } p^{2}p^{\prime 2} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L}\!\left(z_{kl}^{*}p^{\prime }\right) {f}_{L}^{(C)}\!\left(p,p^{\prime }\right),\\ V_{ijkl}^{(\nu \nu ^{\prime }:LS)} &=& A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{\lambda _q=0,1}\sum _{\lambda _Q=0,1} \sum _{L=0}^\infty \sum _{L^{\prime }=0}^\infty (-)^{\lambda _Q} i^{L+L^{\prime }} \widehat{1\!-\!\lambda _q}\widehat{1\!-\!\lambda _Q} \left[\binom{3}{2\lambda _q}\binom{3}{2\lambda _Q}\right]^{\!\frac{1}{2}} \\ && \times \, U_{ijkl,L}^{(LS)} \left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{1,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ && \times \sum _{L_1L_2K} \left( \lambda _q 0 \lambda _Q 0 | L_1 0 \right) \left( 1\!-\!\lambda _q, 0, 1\!-\!\lambda _Q, 0 | L_2 0 \right) \\ && \times \int \!\!\!\!\int dp dp^{\prime } p^{\lambda _q+\lambda _Q+2}p^{\prime 4-\lambda _q-\lambda _Q} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right) {f}_{LL^{\prime }L_1L_2K}^{(LS)}\!\left(p,p^{\prime }\right), \\ \end{eqnarray}$$

(40),(41)

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:T)} &=& A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{\lambda _0=0,2} \sum _{\lambda _1=0}^{\lambda _0} \binom{2\lambda _0+1}{2\lambda _1}^{\!\frac{1}{2}} \frac{\widehat{\lambda _0\!-\!\lambda _1}}{\hat{\lambda }_0}\left( 1 0 1 0 | \lambda _0 0 \right) \\ && \times \sum _{L=0}^\infty \sum _{L^{\prime }=0}^\infty \sum _{K} i^{L+L^{\prime }} \left( \lambda _1 0 K 0 | L 0\right) \left( \lambda _0\!-\!\lambda _1, 0 K 0 | L^{\prime } 0\right) {\left\lbrace \begin{array}{*{10}c}\lambda _0-\lambda _1\, && \lambda _1\, && \lambda _0 \\ L\, && L^{\prime }\, && K \end{array}\right\rbrace } \\ && \times \, U_{ijkl,L}^{(T)} \left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{\lambda _0,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ && \times \int \!\!\!\!\int dp dp^{\prime } p^{\lambda _1+2}p^{\prime \lambda _0-\lambda _1+2} g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right) { f}_{\lambda _0 K}^{(T)}\!\left(p,p^{\prime }\right)\!,\, \, \, \,\,\,\, \end{eqnarray}$$

(42)

$$\begin{eqnarray} V_{ijkl}^{(\nu \nu ^{\prime }:\sigma L)} =&& A_{ijkl}^{(\nu \nu ^{\prime })} \sum _{\lambda _0=0,2} \sum _{L_q=0,2} \sum _{L_Q=0,2} \left(\frac{1}{2}\right)^{\!L_Q} \hat{L}_q \hat{L}_Q \left( 1 0 1 0 | L_q 0 \right) \left( 1 0 1 0 | L_Q 0 \right) {\left\lbrace \begin{array}{*{10}c}1 & \quad 1 & \quad 1 \\ 1 & \quad 1 & \quad 1 \\ L_q & \quad L_Q & \quad \lambda _0 \end{array}\right\rbrace } \\ &&\times \sum _{\lambda _q=0}^{L_q} \sum _{\lambda _Q=0}^{L_Q} (-)^{\lambda _Q} \widehat{L_q\!-\!\lambda _q}\widehat{L_Q\!-\!\lambda _Q} \left[\binom{2L_q+1}{2\lambda _q}\binom{2L_Q+1}{2\lambda _Q}\right]^{\!\frac{1}{2}} \\ &&\times \sum _{L_1L_2} \hat{L}_1 \hat{L}_2 \left( \lambda _q 0 \lambda _Q 0 | L_1 0 \right) \left( L_q\!-\!\lambda _q, 0, L_Q\!-\!\lambda _Q, 0 | L_2 0 \right) {\left\lbrace \begin{array}{*{10}c}\lambda _q && \quad L_q\!-\!\lambda _q && \quad L_q \\ \lambda _Q && \quad L_Q\!-\!\lambda _Q && \quad L_Q \\ L_1 && \quad L_2 && \quad \lambda _0 \end{array}\right\rbrace } \\ &&\times \sum _{L=0}^\infty \sum _{L^{\prime }=0}^\infty \sum _{K} i^{L+L^{\prime }} \left( L_1 0 K 0 | L 0\right) \left( L_2 0 K 0 | L^{\prime } 0\right) {\left\lbrace \begin{array}{*{10}c}L && L^{\prime } && \lambda _0 \\ L_2 && L_1 && K \end{array}\right\rbrace } \\ &&\times \, U_{ijkl,L}^{(T)} \left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{\lambda _0,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ &&\times \int \!\!\!\!\int dp dp^{\prime } p^{\lambda _q+\lambda _Q+2} p^{\prime \, L_q+L_Q-\lambda _q-\lambda _Q+2} \\ &&\times \, g_{\nu \nu ^{\prime }}(p,p^{\prime }) u_n\!\left(p,p^{\prime },\Lambda \right) j_{L}\!\left(z_{ij}p\right) j_{L^{\prime }}\!\left(z_{kl}^{*}p^{\prime }\right) { f}_{L_qL_QK}^{(\sigma L)}\!\left(p,p^{\prime }\right). \end{eqnarray}$$

(43)

When we focus on the central MEs, e.g. the correspondence between Eqs. (25) and (40) is clear: By absorbing the prefactors of Eq. (25) into |$\mathcal {F}^{(C_1)}_L$|⁠, which is newly defined as |${ f}_{L}^{(C)}$|⁠, and by replacing |$D_{ijkl}^{(\sigma )} D_{ijkl}^{(\tau )}-(-)^LD_{jikl}^{(\sigma )} D_{jikl}^{(\tau )}$| with |$U_{ijkl,L}^{(C)}$|⁠, one can obtain Eq. (40). The explicit forms of the MPE functions, |${ f}_{L}^{(C)}$|⁠, |${ f}_{LL^{\prime }L_1L_2K}^{(LS)}$|⁠, |${ f}_{\lambda _0 K}^{(T)}$|⁠, and |${ f}_{L_qL_QK}^{(\sigma L)}$|⁠, are relegated to Appendix  C.

Each term of the MEs involves |$U_{ijkl,L}^{(X)}$|⁠, which depends on the operator form as summarized in Table 1. For instance, we explicitly show the correspondence between |$U_{ijkl,L}^{(X)}$| and the spin-isospin MEs of the C₁, C₅, and C₆ terms derived in the previous sections. They are characterized by the operators, |$\mathbb {1}$|⁠, |$-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)$|⁠, and |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)$|⁠:

$$\begin{eqnarray} U_{ijkl,L}^{(C)} &= D_{ijkl}^{(\sigma )} D_{ijkl}^{(\tau )} -(-)^L D_{jikl}^{(\sigma )} D_{jikl}^{(\tau )}, \qquad \left[\mathrm{for}~\mathbb {1}\right], \end{eqnarray}$$

(44)

$$\begin{eqnarray} U_{ijkl,L}^{(LS)} &= S_{ijkl} \left[D_{ijkl}^{(\tau )}-(-)^L D_{jikl}^{(\tau )}\right] \qquad \left[\mathrm{for}~-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right], \end{eqnarray}$$

(45)

$$\begin{eqnarray} U_{ijkl,L}^{(T)} &= T_{ijkl}^{(\lambda _0)}D_{ijkl}^{(\tau )}-(-)^LT_{jikl}^{(\lambda _0)}D_{jikl}^{(\tau )} \quad \left[\mathrm{for}~\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)\right]. \end{eqnarray}$$

(46)

If the interaction involves the spin-spin operator, |$\boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2$|⁠, we have to replace |$D_{ijkl}^{(\sigma )}$| and |$D_{jikl}^{(\sigma )}$| respectively with |$X_{ijkl}^{(\sigma )}$| and |$X_{jikl}^{(\sigma )}$| defined by

$$\begin{eqnarray} X_{ijkl}^{(\sigma )} &= 2\delta _{m_{\sigma _i} m_{\sigma _l}} \delta _{m_{\sigma _j} m_{\sigma _k}} -\delta _{m_{\sigma _i} m_{\sigma _k}} \delta _{m_{\sigma _j} m_{\sigma _l}}. \end{eqnarray}$$

(47)

The same is true for the isospin MEs, i.e. |$D_{ijkl}^{(\tau )}$| and |$D_{jikl}^{(\tau )}$| respectively become |$X_{ijkl}^{(\tau )}$| and |$X_{jikl}^{(\tau )}$| with the operator |$\boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2$|⁠. Here, |$\boldsymbol{\tau }_i$| is the isospin operator represented by the Pauli matrices. The isospin-isospin operator, |$\boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2$|⁠, does not enter the contact terms of the usual chiral EFT since they are formulated with the choice of |$\mathbb {1}$| and |$\boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2$| based on the Fierz rearrangement freedom [27–29].

3. Visualization of two-body MEs

3.1. Numerical details

To demonstrate the behavior of the two-body MEs formulated in Sect. 2, here we visualize typical values of them. As an example, we choose the leading order (LO) (LO plus NLO) contributions, the potentials of which are given in Appendix B.2 (Appendix B.3). To visualize the MEs some constraints are necessary. First, the Gaussian-center position |$\boldsymbol{Z}_i$| is chosen to be real, and we select the diagonal MEs |$V_{ijij}^{(\nu \nu )}$|⁠. Thus, the MEs can be computed as a function of |$\boldsymbol{r}=\boldsymbol{z}_{ij}=\boldsymbol{z}_{kl}$|⁠, which is the relative distance between two nucleons. Next, as depicted in Fig. 1, a spin-up neutron is settled at the origin, |$\boldsymbol{r}=0$|⁠, and another spin-up proton moves on the x-z plane, i.e. now the three-dimensional vector |$\boldsymbol{r}$| is expressed by |$\boldsymbol{r}=(x,0,z)$|⁠. Note that the spin direction is aligned with the z axis, and therefore, the effect of the tensor contributions originating from the 1π exchange at LO, as well as from 2π exchange and contacts at NLO, can be seen on the x-z plane. Hereafter we use the shorthand notations, |$\left\langle {\hat{V}_{\mathrm{LO}}}\right\rangle$| and |$\left\langle {\hat{V}_{\mathrm{NLO}}}\right\rangle$|⁠, to express the MEs |$V_{ijij}^{(\nu \nu )}$| of the LO and LO-plus-NLO cases, respectively.

Fig. 1.

The geometric configuration for the visualization of the two-body MEs. A spin-up neutron (spin-up proton) is settled at the origin (on the x-z plane).

The parameters we employ here are summarized in Ref. [25], where the LECs are originally taken from [8,30,31]. One can also find the regulator parameters, Λ and n, as well as the constants relevant to the 1π term (the pion mass, pion-decay constant, and axial vector coupling constant). We adopt ν = 0.26 fm⁻² [25].

3.2. MEs at LO and contact-term contributions

Under the conditions described in the previous section, |$\left\langle {\hat{V}_{\mathrm{LO}}}\right\rangle$| can be represented in Fig. 2(a)–(c). The results with three different values of Λ are shown (Λ = 450, 500, and 600 MeV). One finds that the larger Λ is, the less attractive |$\left\langle {\hat{V}_{\mathrm{LO}}}\right\rangle$| is. This is due to the Λ dependence of the LECs at LO, i.e. C_S (C_T), which is responsible for the attraction (repulsion), becomes smaller (larger) as Λ increases [8,30,31].

Fig. 2.

The two-body MEs of the chiral interaction at LO with the cutoff (a) Λ = 450, (b) 500, and (c) 600 MeV, with x and z corresponding to the axes shown in Fig. 1. The MEs of the 1π-exchange potential (OPEP)-tensor contributions are also displayed with (d) Λ = 450, (e) 500, and (f) 600 MeV. See the text for detail.

Figure 2(a)–(c) shows that the spin-aligned neutron–proton pair feels the largest attraction at |$\boldsymbol{r}=0$|⁠. This is because the LO potential is designed to simulate nucleon–nucleon scattering at very low momentum, where the neutron–proton interaction of the triplet-s state is attractive as deduced from the scattering phase shift [32]. Furthermore, the attractive MEs at |$\boldsymbol{r}=0$| can be shown analytically. First, we ignore the regulator with Λ → ∞. Indeed, the largest attraction at |$\boldsymbol{r}=0$| can be seen independently of Λ, and therefore, the regulator does not play an essential role for the present discussion. Then, the MEs of the LO-contact term within the configuration of Fig. 1 can be simplified as

$$\begin{eqnarray} \left\langle {\hat{V}_{\textrm {LO-ct}}}\right\rangle \xrightarrow {\Lambda \rightarrow \infty } \left(\frac{\nu }{\pi }\right)^{\!\frac{3}{2}} (C_S^{np}+C_T^{np})\exp \!\left[-\nu r^2\right]. \end{eqnarray}$$

(48)

Note that the LECs at LO have the charge dependence. The superscript np stands for the LECs for the neutron–proton pair. We find that the condition |$C_S^{np}+C_T^{np} \lt 0$| is satisfied by the LECs employed here. For example, |$C_S^{np}+C_T^{np} =-0.011374\times 10^4$| GeV⁻² for Λ = 450 MeV [30,31]. Thus we can show that Eq. (48) has a minimum at |$\boldsymbol{r}=0$|⁠.

We should mention a role played by the 1π interaction at the origin. Within the configuration of Fig. 1 at |$\boldsymbol{r}=0$|⁠, the system is a triplet-even state, and hence, the central term of the 1π-exchange interaction is repulsive, while the 1π-tensor term has no contributions there [see Fig. 2(d)–(f)].

Note that the 1π-exchange potential (OPEP) considered here contains the short-range delta function in coordinate space. This short-range term is often subtracted by hand since 1π exchange should be responsible for the long-range part of the two-nucleon force (see, e.g. Ref. [33]), and thus, the central term of this “subtracted” OPEP is attractive for the triplet-even state. Instead, in the chiral EFT, the short-range term of the OPEP remains included but the addition of the LO-contact term as a counterterm suppresses the short-range repulsion of the OPEP. Consequently, for the triplet-even state, the OPEP-central term of the chiral EFT alone is repulsive, but the whole LO potential is attractive.

The repulsion by the 1π-contact term is moderate compared to that by the LO-contact terms. Specifically, for Λ = 450 MeV, the absolute value of the 1π-central ME at the origin is less than |$30\%$| of that of the contacts. Therefore, as addressed above, the minimum of the MEs at |$\boldsymbol{r}=0$| can be basically explained by the contributions from the LO-contact terms.

3.3. One-pion-exchange-tensor contributions to MEs

By carefully observing Figs. 2(a)–(c), one finds that the MEs are not symmetric with respect to the z = x line. For example, the contour line of −1 MeV in Fig. 2(a) crosses the z and x axes at ∼4 and ∼3 fm, respectively. Hence, the MEs are more attractive in the z > x region and vice versa. This asymmetry is due to the tensor contributions of the 1π-exchange term.

Analytically, the asymmetry can be understood as follows. The tensor operator defined by Eq. (37) has anisotropy. As well known (see, e.g. Ref. [33]), the OPEP in the coordinate-space representation, which is the Fourier transform of Eq. (B5), involves the tensor operator,

$$\begin{eqnarray} \mathcal {S}_{12}(\hat{\boldsymbol{r}}) &= 2\left[3\left(\boldsymbol{S}\cdot \hat{\boldsymbol{r}}\right)^2 -\boldsymbol{S}^2\right]. \end{eqnarray}$$

(49)

The direction of |$\boldsymbol{S}$| is now aligned with the z axis, and therefore, the attraction by the tensor contributions becomes stronger for the z > x region on the x-z plane. Note that the expectation value of the isospin operator |$\boldsymbol{\tau }_1\cdot \boldsymbol{\tau }_2$| appearing in the OPEP is −3 for the isoscalar state at |$\boldsymbol{r}=0$|⁠. For finite |$\boldsymbol{r}$|⁠, isoscalar-isovector mixing occurs and it is taken into account in the present calculations.

Moreover, one can realize that the bipolar spherical harmonics, which is involved in the tensor MEs given by Eq. (42), is asymmetric with respect to the exchange of x and z. Indeed, under the present conditions, the bipolar spherical harmonics of the tensor component (λ₀ = 2) becomes

$$\begin{eqnarray} &&\left[Y_{L}\!\left(z_{ij,x},z_{ij,y},z_{ij,z}\right) \otimes Y_{L^{\prime }}\!\left(z_{kl,x}^{*},z_{kl,y}^{*},z_{kl,z}^{*}\right)\right]_{\lambda _0,m_{\sigma _i}+m_{\sigma _j}-m_{\sigma _k}-m_{\sigma _l}} \\ &&\quad \rightarrow \frac{1}{8\pi }\hat{L}\hat{L}^{\prime }\left(L0L^{\prime }0|20\right) \frac{2z^2-x^2}{x^2+z^2}, \end{eqnarray}$$

(50)

for which we notice the asymmetry with respect to the operation x↔z.

Numerically, the asymmetry of the tensor MEs |$\left\langle {\hat{V}_{\mathrm{LO}}^{(T)}}\right\rangle$| can be clearly seen in Fig. 2(d)–(f), where only the tensor component λ₀ = 2 of Eq. (42) for the OPEP is taken into account. Figure 2(d)–(f) displays that the tensor contributions are attractive (repulsive) on the z (x) axis, as expected from the structure of the tensor operator, Eq. (49).

At |$\boldsymbol{r}=0$|⁠, even though the radial part of the tensor component of the OPEP diverges (see, e.g. Ref. [34]), the tensor MEs must be zero, since only the s wave contributes when two nucleons contact with each other. Note that the divergence of the OPEP-tensor part does not matter in the present calculations since the MEs are computed in momentum space, where the corresponding high-momentum component of the OPEP is suppressed by the regulator. The position of the extrema of |$\left\langle {\hat{V}_{\mathrm{LO}}^{(T)}}\right\rangle$| depends on ν. If we compute the tensor MEs with larger ν, the minimum (maximum) point moves in the direction that z (x) becomes smaller, and the attractive pockets in Fig. 2(d)–(f) become deeper with larger ν. This is consistent with the ν dependence of tensor-force contributions to the ⁴He energy reported in Ref. [35]; the energy gain by a tensor force increases when ν is greater than its typical value ν = 0.25 fm⁻² for ⁴He [13,15], although the total energy is saturated due to the compensation by the energy loss of the kinetic term for such a large ν.

Also one finds from Fig. 2(d)–(f) that both attraction and repulsion are enhanced as Λ increases. This results from the operator form of the tensor potential. In momentum space, the tensor component of the OPEP can be written as

$$\begin{eqnarray} \mathcal {O}_{12}(\boldsymbol{q}) = \frac{1}{3\left\lbrace 1+(m_\pi /q)^2\right\rbrace }\mathcal {S}_{12}(\hat{\boldsymbol{q}}), \end{eqnarray}$$

(51)

which can be obtained from the operator, |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)\left(q^2+m_\pi ^2\right)^{-1}$| through the irreducible-tensor representation as Eq. (35) with λ₀ = 2. When large q contributes (this is the case for larger Λ), the effect of |$\mathcal {O}_{12}$| is enhanced. This is consistent with the Λ dependence of |$\left\langle {\hat{V}_{\mathrm{LO}}^{(T)}}\right\rangle$|⁠.

3.4. MEs at NLO

Now we increase the order of the chiral EFT up to NLO to visualize the two-body MEs, |$\left\langle {\hat{V}_{\mathrm{NLO}}}\right\rangle$|⁠, within the same two-nucleon configuration as that in the previous section. Figure 3(a) shows |$\left\langle {\hat{V}_{\mathrm{NLO}}}\right\rangle$| for the regulator cutoff Λ = 450 MeV, where the x and z axes are given in Fig. 1. The NLO MEs of the tensor contributions denoted by |$\left\langle {\hat{V}_{\mathrm{NLO}}^{(T)}}\right\rangle$| are extracted from |$\left\langle {\hat{V}_{\mathrm{NLO}}}\right\rangle$|⁠, and depicted in Fig. 3(b). Similar results for Λ = 500 and 600 MeV are shown in Figs. 4 and 5, respectively. Note that even though the C₅ term, one of the NLO contact terms, enters the NLO potential as the SO contributions, it plays no roles on the MEs visualized in Figs. 3–5. This is because the direction of |$\boldsymbol{S}$| is aligned with the z axis and also the Gaussian center is real. Indeed one can show analytically that the bipolar spherical harmonics involved in Eq. (32) becomes zero in the x-z plane.

Fig. 3.

The two-body MEs of the chiral nucleon–nucleon force at NLO for (a) the full contributions and (b) tensor contributions only, where the x and z axes are given in Fig. 1. The regulator cutoff is Λ = 450 MeV.

Fig. 4.

Same as Fig. 3 but for Λ = 500 MeV.

Fig. 5.

Same as Fig. 3 but for Λ = 600 MeV.

As argued in Sect. 3.2, the LO MEs at |$\boldsymbol{r}=0$| are attractive since the LO potential is tailored to very-low-momentum scattering of two nucleons. In contrast to the LO MEs, the NLO MEs in Figs. 3(a), 4(a), and 5(a) are all repulsive at |$\boldsymbol{r}=0$|⁠. This repulsive nature stems from the high-momentum scattering described by the NLO potentials dependent on the square of the momentum (see Appendix B.3). Furthermore, one finds that the more Λ increases, the more the repulsion of the MEs enhances. This Λ dependence of the MEs reflects the property of the neutron–proton interaction, which turns from attractive to repulsive around 400 MeV/c of the relative momentum as indicated by the neutron–proton scattering phase shift of the triplet-s state [32].

In Figs. 3(b), 4(b), and 5(b), the tensor MEs of the LO-plus-NLO contributions, |$\left\langle {\hat{V}_{\mathrm{NLO}}^{(T)}}\right\rangle$|⁠, clearly show asymmetry with respect to the z = x line. The origin of this asymmetry is the same as that in the OPEP-tensor MEs, as already discussed in Sect. 3.3. However, |$\left\langle {\hat{V}_{\mathrm{NLO}}^{(T)}}\right\rangle$| has the opposite sign compared to that of |$\left\langle {\hat{V}_{\mathrm{LO}}^{(T)}}\right\rangle$|⁠, i.e. the repulsive peak (attraction pocket) of |$\left\langle {\hat{V}_{\mathrm{NLO}}^{(T)}}\right\rangle$| localizes on the z (x) axis, and vice versa for |$\left\langle {\hat{V}_{\mathrm{LO}}^{(T)}}\right\rangle$|⁠. This is because the NLO potentials do not depend on the isospin-isospin operator, |$\boldsymbol{\tau }_1\cdot \boldsymbol{\tau }_2$|⁠, which induces the change of the sign of the MEs, and also because the LECs, C₆ and C₇, are negative in the present parameterization [8,30,31].

4. Conclusions and perspectives

We have presented the formalism of the two-body MEs of the two-nucleon forces in momentum space derived from the chiral EFT based on the single-nucleon wave functions expressed by the Gaussian-wave packet. Such MEs are relevant to many-body calculations like AMD, which can be applicable, for instance, to efficiently describe nuclear cluster structures. We adopt the chiral potentials at N³LO based on the nonlocal regularization.

We have visualized the MEs formulated in this paper by selecting the spin-up neutron and spin-up proton pair, for which the tensor contributions can be seen. As an example, the chiral potentials at LO have been chosen. We have addressed the cutoff dependence of the MEs and the origin of the ME extrema, as well as the individual contributions of the central and tensor forces.

As a next step of this work, we are now implementing this formalism to AMD. Then, benchmark calculations will be performed for light nuclei. The inclusion of the chiral three-nucleon force into AMD is one of the important tasks, and the extension of the framework in this direction is also ongoing.

Acknowledgement

The author thanks N. Itagaki and M. Kimura for fruitful discussions and useful comments. This work was supported by Japan Society for the Promotion of Science KAKENHI with Grant Number JP21K13919. The calculations were carried out using the computer facilities at Yukawa Institute for Theoretical Physics, Kyoto University.

Appendix A. Slater determinants and their superposition

A.1. Slater determinants

A single-Slater determinant for a system of mass number A is characterized by a single value of ν and a set of the generator coordinate |$\mathcal {Z}_n$|⁠:

$$\begin{eqnarray} \left|\Psi _\nu (\mathcal {Z}_n)\right\rangle & = \sqrt{A!}\hat{\mathcal {A}}_{A} \left|\varphi _1^{(\nu )} \varphi _2^{(\nu )} \cdots \varphi _i^{(\nu )} \cdots \varphi _A^{(\nu )}\right\rangle , \end{eqnarray}$$

(A1)

$$\begin{eqnarray} \mathcal {Z}_n&=\left\lbrace \boldsymbol{Z}_1,\boldsymbol{Z}_{2},\cdots ,\boldsymbol{Z}_{A}\right\rbrace _n, \end{eqnarray}$$

(A2)

with the A-body antisymmetrizer |$\hat{\mathcal {A}}_{A}$|⁠. Then, many-body states are obtained through the generator-coordinate method (GCM) [26,36,37] after the parity and angular-momentum projections:

$$\begin{eqnarray} \left|\Psi _{\mathrm{GCM}}\right\rangle = \sum _{n\nu }c_n^{(\nu )} \left|\Psi _{\nu MK}^{J\pm }(\mathcal {Z}_n)\right\rangle , \end{eqnarray}$$

(A3)

$$\begin{eqnarray} \left|\Psi _{\nu MK}^{J\pm }(\mathcal {Z}_n)\right\rangle &= \hat{\mathcal {P}}_{MK}^J \hat{\mathcal {P}}^{\pm } \left|\Psi _{\nu }(\mathcal {Z}_n)\right\rangle , \end{eqnarray}$$

(A4)

where the projection operators |$\hat{\mathcal {P}}^{\pm }$| and |$\hat{\mathcal {P}}_{MK}^J$| are given explicitly in the next sections.

The coefficient |$c_n^{(\nu )}$| is obtained numerically by solving the generalized eigenvalue problem,

$$\begin{eqnarray} \sum _{m\nu ^{\prime }}\left(H_{nm}^{(\nu \nu ^{\prime })}-EN_{nm}^{(\nu \nu ^{\prime })}\right)c_m^{(\nu ^{\prime })}=0, \end{eqnarray}$$

(A5)

with the energy eigenvalue E. The single-particle wave function expressed by the Gaussian-wave packet, |${\left|\Psi _\nu (\mathcal {Z}_n)\right\rangle }$|⁠, can be separated into the CM and intrinsic wave functions. Thus, the MEs of the norm and Hamiltonian |$\hat{H}$| relevant to the intrinsic structure are given by

$$\begin{eqnarray} {\left\lbrace \begin{array}{*{10}c}N_{nm}^{(\nu \nu ^{\prime })}\\ H_{nm}^{(\nu \nu ^{\prime })} \end{array}\right\rbrace } &= \frac{ \left\langle {\Psi _{\nu MK}^{J\pm }(\mathcal {Z}_n) \left| {\left\lbrace \begin{array}{*{10}c}\mathbb {1}\\ \hat{H} \end{array}\right\rbrace } \right| \Psi _{\nu ^{\prime } MK}^{J\pm }(\mathcal {Z}_m) }\right\rangle }{\left\langle {\Psi _\nu ^{(\mathrm{CM})}}\,\vert \,{\Psi _{\nu ^{\prime }}^{(\mathrm{CM})}}\right\rangle }, \end{eqnarray}$$

(A6)

where the CM-wave function |$\left|\Psi _\nu ^{(\mathrm{CM})}\right\rangle$| is given by

$$\begin{eqnarray} \Psi _\nu ^{(\mathrm{CM})}\!\left(\boldsymbol{r}_{\mathrm{G}}-\boldsymbol{Z}_{\mathrm{G}}\right) = \left(\frac{2\nu }{\pi }\right)^{\!\frac{3}{4}} \exp \!\left[-A\nu \left(\boldsymbol{r}_{\mathrm{G}}-\boldsymbol{Z}_{\mathrm{G}}\right)^2\right], \end{eqnarray}$$

(A7)

with |$\boldsymbol{r}_{\mathrm{G}}=\sum _{i=1}^A \boldsymbol{r}_i/A$| and |$\boldsymbol{Z}_{\mathrm{G}}=\sum _{i=1}^A \boldsymbol{Z}_i/A$|⁠.

The Hamiltonian |$\hat{H}$| is written as

$$\begin{eqnarray} \hat{H}=\hat{T}-\hat{T}_\mathrm{G}+\hat{V} \end{eqnarray}$$

(A8)

where |$\hat{T}$| is the kinetic-energy operator expressed by the sum of the one-body kinetic-energy operator and |$\hat{T}_\mathrm{G}$| is the kinetic-energy operator of the CM system. The expectation values of these operators can be calculated analytically. For example, as shown in Ref. [14], one finds

$$\begin{eqnarray} \frac{\left\langle {\Psi _\nu (\mathcal {Z}_n)\left|\hat{T}_\mathrm{G}\right|\Psi _{\nu ^{\prime }}(\mathcal {Z}_m)}\right\rangle }{\left\langle {\Psi _\nu (\mathcal {Z}_n)}\,\vert \,{\Psi _{\nu ^{\prime }}(\mathcal {Z}_m)}\right\rangle } = \frac{3\hbar ^2\nu \nu ^{\prime }}{m_N(\nu +\nu ^{\prime })}. \end{eqnarray}$$

(A9)

Note that Eq. (A9) is obtained under the assumption, |$\boldsymbol{Z}_{\mathrm{G}}=0$|⁠, for both bra and ket states. The interaction operator |$\hat{V}$| consists of nuclear and Coulomb parts, and if it involves the chiral two-nucleon force, |$H_{nm}^{(\nu \nu ^{\prime })}$| can be written in terms of the two-body MEs given in Sect. 2.

A.2. Parity projection

First, we introduce the parity inversion operator |$\hat{\mathcal {P}}_\pi$|⁠, which inverts the sign of the Gaussian center position of the single-particle states as

$$\begin{eqnarray} \hat{\mathcal {P}}_\pi \left|\Psi _\nu (\mathcal {Z})\right\rangle &=\left|\Psi _\nu (-\mathcal {Z})\right\rangle , \end{eqnarray}$$

(A10)

$$\begin{eqnarray} -\mathcal {Z} &=\left\lbrace -\boldsymbol{Z}_1,-\boldsymbol{Z}_{2},\cdots ,-\boldsymbol{Z}_{A}\right\rbrace . \end{eqnarray}$$

(A11)

Note that the subscript n is omitted from |$\mathcal {Z}$| for simplicity. Thus, the parity projected operator and parity projected states are respectively defined by

$$\begin{eqnarray} \hat{\mathcal {P}}^{\pm } = \frac{1\pm \hat{\mathcal {P}}_\pi }{2\alpha _\pm }, \end{eqnarray}$$

(A12)

$$\begin{eqnarray} \left|\Psi _\nu ^{\pm }\right\rangle =\hat{\mathcal {P}}^{\pm } \left|\Psi _\nu (\mathcal {Z})\right\rangle = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{2\alpha _+}\left[\left|\Psi _\nu (\mathcal {Z})\right\rangle +\left|\Psi _\nu (-\mathcal {Z})\right\rangle \right], \\ \frac{1}{2\alpha _-}\left[\left|\Psi _\nu (\mathcal {Z})\right\rangle -\left|\Psi _\nu (-\mathcal {Z})\right\rangle \right], \end{array}\right. \end{eqnarray}$$

(A13)

where α_± is the normalization coefficient in association with the positive (+) or negative (−) parity.

A.3. Angular-momentum projection

The angular-momentum-projected state |$\left|\Psi _{\nu MK}^{J}(\mathcal {Z})\right\rangle$| is defined with the angular-momentum-projection operator |$\hat{\mathcal {P}}_{MK}^J$| by

$$\begin{eqnarray} \left|\Psi _{\nu MK}^{J}(\mathcal {Z})\right\rangle &= \hat{\mathcal {P}}_{MK}^J \left|\Psi _{\nu }(\mathcal {Z})\right\rangle , \end{eqnarray}$$

(A14)

$$\begin{eqnarray} \hat{\mathcal {P}}_{MK}^J = \frac{2J+1}{8\pi ^2} \int d\Omega D_{MK}^{J*}(\Omega )\hat{R}(\Omega ). \end{eqnarray}$$

(A15)

Here, M is the z-component of the total angular momentum J in the laboratory frame, where the rotational symmetry is restored, while K is that for the body-fixed (intrinsic) frame. The three-dimensional Euler angle Ω appears as an argument of the Wigner D-function |$D_{MK}^{J}$| and the rotation operator |$\hat{R}$|⁠, which is associated with the rotation in spatial and spin spaces. The integration over Ω in Eq. (A15) can be performed numerically. Also we simply denote |$\mathcal {Z}$| without the subscript n.

The operation of |$\hat{R}$| results in the rotation of the single-particle wave function. As a result, for the spatial part, we just have to rotate |$\boldsymbol{Z}_{ij}$| and |$\boldsymbol{z}_{ij}$| as

$$\begin{eqnarray} \boldsymbol{Z}_{ij} &\rightarrow \boldsymbol{Z}_{ij}^{(R)} = \frac{1}{2}\left[ \hat{R}(\Omega )\boldsymbol{Z}_i +\hat{R}(\Omega )\boldsymbol{Z}_j \right], \end{eqnarray}$$

(A16)

$$\begin{eqnarray} \boldsymbol{z}_{ij} &\rightarrow \boldsymbol{z}_{ij}^{(R)} = \hat{R}(\Omega )\boldsymbol{Z}_i -\hat{R}(\Omega )\boldsymbol{Z}_j. \end{eqnarray}$$

(A17)

As regards the rotation of the spin wave function, the coefficients α_i and β_i in Eq. (3) are replaced by

$$\begin{eqnarray} {\begin{pmatrix}\alpha _i\\ \beta _i \end{pmatrix}} \rightarrow {\begin{pmatrix}\alpha _i^{\prime }\\ \beta _i^{\prime } \end{pmatrix}} = {\begin{pmatrix}D_{\frac{1}{2}\frac{1}{2}}^{\frac{1}{2}}(\Omega ) &&\,\,\, D_{\frac{1}{2},-\frac{1}{2}}^{\frac{1}{2}}(\Omega )\\ D_{-\frac{1}{2},\frac{1}{2}}^{\frac{1}{2}}(\Omega ) &&\,\,\, D_{-\frac{1}{2},-\frac{1}{2}}^{\frac{1}{2}}(\Omega ) \end{pmatrix}} {\begin{pmatrix}\alpha _i\\ \beta _i \end{pmatrix}}, \end{eqnarray}$$

(A18)

as well as α_j and β_j being replaced by |$\alpha _j^{\prime }$| and |$\beta _j^{\prime }$|⁠, respectively.

Appendix B. Potentials derived from chiral EFT

B.1. Overview

We adopt a high-precision two-body potential based on the chiral EFT at N³LO [8,10–12] as v_2N in Eq. (12). At this order, the potential can be written order by order as

$$\begin{eqnarray} v_{2N}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &= v_{2N}^{(\mathrm{LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) +v_{2N}^{(\mathrm{NLO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) +v_{2N}^{(\mathrm{N^2LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) +v_{2N}^{(\mathrm{N^3LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right). \end{eqnarray}$$

(B1)

The superscript stands for the chiral-expansion power, i.e. n_χ = 0, 2, 3, and 4 for LO, NLO, next-to-next-to-leading order (N²LO), and N³LO, respectively.

In this work, we employ the nonlocal regularization with the regulator

$$\begin{eqnarray} u_n\!\left(p,p^{\prime },\Lambda \right) = \exp \!\left[-\!\left(\frac{p}{\Lambda }\right)^{\!2n}-\!\left(\frac{p^{\prime }}{\Lambda }\right)^{\!2n}\right]. \end{eqnarray}$$

(B2)

Here Λ is the cutoff momentum. Thus, the potential depends on the relative momenta, |$\boldsymbol{p}$| and |$\boldsymbol{p}^{\prime }$|⁠, of the initial and final channels, respectively, which are related to the average momentum |$\boldsymbol{Q}$| and transferred momentum |$\boldsymbol{q}$| by

$$\begin{eqnarray} {\begin{pmatrix}\boldsymbol{Q} \\ \boldsymbol{q} \end{pmatrix}} = \mathcal {U} {\begin{pmatrix}\boldsymbol{p}^{\prime } \\ \boldsymbol{p} \end{pmatrix}}, \end{eqnarray}$$

(B3)

where |$\mathcal {U}$| is given by Eq. (7). Note that every potential v_2N appearing in this paper involves the prefactor 1/(2π)³, which originates from the normalization convention, |${\left\langle {\boldsymbol{p}}\,\vert \,{\boldsymbol{p}^{\prime }}\right\rangle =\delta (\boldsymbol{p}-\boldsymbol{p}^{\prime })}$|⁠. A similar prefactor for a potential of the chiral three-nucleon force can be found in Refs. [38,39].

At each order, the potentials consist of the 1π exchange term v_1π, the two-pion (2π) exchange term |$v_{2\pi }^{(n_\chi )}$|⁠, and the contact term |$v_{\mathrm{ct}}^{(n_\chi )}$|⁠. We see these contributions in the next sections.

B.2. LO

The LO potential reads

$$\begin{eqnarray} v_{2N}^{(\mathrm{LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& v_{1\pi }\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right)+v_{\mathrm{ct}}^{(0)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right), \\ v_{1\pi }\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& -\frac{1}{(2\pi )^3}\!\left(\frac{g_A}{2f_\pi }\right)^{\!2} u_n\!\left(p,p^{\prime },\Lambda \right) \frac{\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right)}{q^2+m_\pi ^2} \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2,\\ v_{\mathrm{ct}}^{(0)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& \frac{1}{(2\pi )^3} u_n\!\left(p,p^{\prime },\Lambda \right) \left(C_S+C_T\, \boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2\right), \end{eqnarray}$$

(B4),(B5),(B6)

where g_A is the axial vector coupling constant, f_π is the pion-decay constant, m_π is the average pion mass, and the LECs are C_S and C_T at LO. The spin and isospin operators, |$\boldsymbol{\sigma }_i$| and |$\boldsymbol{\tau }_i$|⁠, respectively, are represented by the Pauli matrices.

B.3. NLO

At NLO, the potential is given by

$$\begin{eqnarray} v_{2N}^{(\mathrm{NLO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& v_{2\pi }^{(2)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right)+v_{\mathrm{ct}}^{(2)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right), \hphantom{\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right)\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right)} \\ v_{2\pi }^{(2)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& u_n\!\left(p,p^{\prime },\Lambda \right) \Big [ W_C^{(2)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2 +V_S^{(2)}(q,Q) \, \boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2 \\ &&+V_T^{(2)}(q,Q) \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) \Big ],\\ v_{\mathrm{ct}}^{(2)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& \frac{1}{(2\pi )^3}u_n\!\left(p,p^{\prime },\Lambda \right) \Big [C_1 q^2 +C_2 Q^2 +\left(C_3 q^2 +C_4 Q^2\right)\boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2 \\ && +\,C_5\left[-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] +C_6\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) +C_7\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{Q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{Q}\right)\Big ], \end{eqnarray}$$

(B7),(B8),(B9)

where |$\boldsymbol{S}=(\boldsymbol{\sigma }_1+\boldsymbol{\sigma }_2)/2$| and C_i are the LECs at NLO. The 2π term contains

$$\begin{eqnarray} W_C^{(2)}(q,Q) &=& -\frac{L(q)}{3(4\pi )^5 f_\pi ^4} \\ && \times \left[4m_\pi ^2\left(5g_A^4-4g_A^2-1\right) +q^2\left(23g_A^4-10g_A^2-1\right) +48g_A^4 m_\pi ^4w^{-2}\right],\\ V_S^{(2)}(q,Q) &=& \frac{6g_A^4}{(4\pi )^5 f_\pi ^4}q^2L(q), \\ V_T^{(2)}(q,Q) &=& -\frac{6g_A^4}{(4\pi )^5 f_\pi ^4}L(q), \end{eqnarray}$$

(B10),(B11),(B12)

with

$$\begin{eqnarray} w = \sqrt{4m_\pi ^2 +q^2}, \end{eqnarray}$$

(B13)

$$\begin{eqnarray} L(q) = \frac{w}{q}\ln \frac{w+q}{2m_\pi }. \end{eqnarray}$$

(B14)

B.4. N²LO

At N²LO, we have the potential

$$\begin{eqnarray} v_{2N}^{(\mathrm{N^2LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& v_{2\pi }^{(3)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) \\ &=& u_n\!\left(p,p^{\prime },\Lambda \right) \Bigg [V_C^{(3)}(q,Q)+W_C^{(3)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2 \\ && +\,\left[V_S^{(3)}(q,Q)+W_S^{(3)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2 \right] \boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2 \\ && +\,\left[V_{LS}^{(3)}(q,Q) +W_{LS}^{(3)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2\right] \left[-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \\ && +\,\left[V_{T}^{(3)}(q,Q) +W_T^{(3)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2\right] \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) \Bigg ]. \end{eqnarray}$$

(B15)

Here, the potentials involve the central terms,

$$\begin{eqnarray} V_C^{(3)}(q,Q) &=& \frac{6g_A^2}{(4\pi )^4 f_\pi ^4} \Bigg [ \frac{g_A^2 m_\pi ^5}{16 m_N w^2} -\left\lbrace 2m_\pi ^2\left(2c_1-c_3\right) -q^2\left(c_3 +\frac{3g_A^2}{16m_N}\right)\right\rbrace \tilde{w}^2 A(q) \\ &&-\frac{g_A^2}{16 m_N}\left\lbrace m_\pi w^2 +\tilde{w}^4 A(q)\right\rbrace \Bigg ], \\ W_C^{(3)}(q,Q) &=& \frac{g_A^2}{4(4\pi )^4 m_N f_\pi ^4} \Bigg [ 3g_A^2 m_\pi ^5 w^{-2} -\left\lbrace 4m_\pi ^2 +2q^2-g_A^2\left(4m_\pi ^2+3q^2\right)\right\rbrace \tilde{w}^2 A(q) \\ &&+\,\,g_A^2\left\lbrace m_\pi w^2 +\tilde{w}^4 A(q)\right\rbrace \Bigg ],\\ V_S^{(3)}(q,Q) &=& -\frac{3g_A^4}{(8\pi )^4 m_N f_\pi ^4} \Big [ 3q^2 \tilde{w}^2 A(q) +q^2 \left\lbrace m_\pi +w^2 A(q)\right\rbrace \Big ], \\ W_S^{(3)}(q,Q) &=& \frac{g_A^2}{(4\pi )^4 f_\pi ^4} \Bigg [ q^2 A(q) \left\lbrace \left(c_4+\frac{1}{4m_N}\right)w^2 -\frac{g_A^2}{8m_N}\left(10m_\pi ^2+3q^2\right)\right\rbrace \\ &&+\,\,\frac{g_A^2}{8m_N}q^2 \left\lbrace m_\pi +w^2 A(q)\right\rbrace \Bigg ], \end{eqnarray}$$

(B16),(B17),(B18),(B19)

the SO terms,

$$\begin{eqnarray} V_{LS}^{(3)}(q,Q) = \frac{3g_A^4}{(4\pi )^4 m_N f_\pi ^4} \tilde{w}^2 A(q), \end{eqnarray}$$

(B20)

$$\begin{eqnarray} W_{LS}^{(3)}(q,Q) = \frac{g_A^2\left(1-g_A^2\right)}{(4\pi )^4 m_N f_\pi ^4} w^2 A(q), \end{eqnarray}$$

(B21)

and the tensor terms,

$$\begin{eqnarray} V_T^{(3)}(q,Q) = -\frac{1}{q^2} V_S^{(3)}(q,Q), \end{eqnarray}$$

(B22)

$$\begin{eqnarray} W_T^{(3)}(q,Q) = -\frac{1}{q^2} W_S^{(3)}(q,Q), \end{eqnarray}$$

(B23)

with

$$\begin{eqnarray} \tilde{w} = \sqrt{2m_\pi ^2 +q^2}, \end{eqnarray}$$

(B24)

$$\begin{eqnarray} A(q) = \frac{1}{2q}\arctan \frac{q}{2m_\pi }. \end{eqnarray}$$

(B25)

The LECs c_i enter the 2π term at N²LO, and m_N denotes the average nucleon mass.

B.5. N³LO

At N³LO, the potentials are represented as

$$\begin{eqnarray} v_{2N}^{(\mathrm{N^3LO})}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& v_{2\pi }^{(4)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right)+v_{\mathrm{ct}}^{(4)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right), \\ v_{2\pi }^{(4)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) &=& u_n\!\left(p,p^{\prime },\Lambda \right) \Big [ V_{C}^{(4)}(q,Q)+W_{C}^{(4)}(q,Q) \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2 \\ && +\,V_{S}^{(4)}(q,Q) \boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2 +W_{S}^{(4)}(q,Q) \left(\boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2\right) \left(\boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2\right) \\ && +\,\left[V_{LS}^{(4)}(q,Q) +W_{LS}^{(4)}(q,Q) \, \boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2\right] \left[-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \\ && +\,\left[V_T^{(4)}(q,Q)+W_T^{(4)}(q,Q)\boldsymbol{\tau }_1 \cdot \boldsymbol{\tau }_2\right] \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) \\ && +\,V_{\sigma L}^{(4)}(q,Q) \left[\boldsymbol{\sigma }_1\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \left[\boldsymbol{\sigma }_2\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \Big ], \end{eqnarray}$$

(B26),(B27)

$$\begin{eqnarray} v_{\mathrm{ct}}^{(4)}\!\left(\boldsymbol{p},\boldsymbol{p}^{\prime }\right) =&& \frac{1}{(2\pi )^3} u_n\!\left(p,p^{\prime },\Lambda \right) \Big [ D_1 q^4 +D_2 Q^4 +D_3q^2Q^2 +D_4 \left(\boldsymbol{q}\times \boldsymbol{Q}\right)^2 \\ &&+ \left[D_5 q^4 +D_6 Q^4 +D_7q^2Q^2 +D_8 \left(\boldsymbol{q}\times \boldsymbol{Q}\right)^2\right] \boldsymbol{\sigma }_1 \cdot \boldsymbol{\sigma }_2 \\ &&+ \left[D_9 q^2 +D_{10} Q^2\right] \left[-i\boldsymbol{S}\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \\ &&+ \left[D_{11} q^2 +D_{12} Q^2\right] \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{q}\right) \\ &&+ \left[D_{13} q^2 +D_{14} Q^2\right] \left(\boldsymbol{\sigma }_1\cdot \boldsymbol{Q}\right)\left(\boldsymbol{\sigma }_2\cdot \boldsymbol{Q}\right) \\ &&+\, D_{15} \left[\boldsymbol{\sigma }_1\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right] \left[\boldsymbol{\sigma }_2\cdot \left(\boldsymbol{q}\times \boldsymbol{Q}\right)\right]\Big]. \end{eqnarray}$$

(B28)

The LECs of the contact terms at N³LO are D_i.

Following Ref. [8], the 2π-exchange potential at N³LO is categorized into several terms, i.e. |$c_i^2$|⁠, c_i/m_N, and |$m_N^{-2}$| contributions in association with one-loop diagrams, and also two-loop (2L) contributions, where c_i stands for one of the LECs. Thus, potentials in Eq. (B27) can be decomposed further. The purely central terms of |$v_{2\pi }^{(4)}$| read

$$\begin{eqnarray} V_C^{(4)}(q,Q) &=& V_C^{(c_i^2)}(q,Q) +V_C^{(c_i/m_N)}(q,Q) +V_C^{(m_N^{-2})}(q,Q) +V_C^{(\mathrm{2L})}(q,Q),\\ V_C^{(c_i^2)}(q,Q) &=& \frac{3}{4(2\pi )^5 f_\pi ^4} L(q) \left[\left(\frac{c_2}{6}w^2+c_3\tilde{w}^2 -4c_1m_\pi ^2\right)^2 +\frac{c_2^2}{45}w^4\right], \end{eqnarray}$$

(B29),(B30)

$$\begin{eqnarray} V_C^{(c_i/m_N)}(q,Q) &=& -\frac{g_A^2}{8(2\pi )^5 m_N f_\pi ^4} L(q) \Big [ \left\lbrace (c_2-6c_3)q^4 +4(6c_1+c_2-3c_3)q^2m_\pi ^2 \right. \\ &&\left. +\,6(c_2-2c_3)m_\pi ^4 +24(2c_1+c_3)m_\pi ^6w^{-2}\right\rbrace \Big ], \\ V_C^{(m_N^{-2})}(q,Q) &=& -\,\frac{g_A^4}{8(2\pi )^5 m_N^2 f_\pi ^4} \Big [ L(q) \left(2m_\pi ^8 w^{-4} +8m_\pi ^6w^{-2}-q^4-2m_\pi ^4\right) \\ &&+\,\frac{1}{2}m_\pi ^6w^{-2} \Big ], \\ V_C^{(\mathrm{2L})}(q,Q) &=& \frac{3g_A^4}{8(4\pi )^5 f_\pi ^6} \tilde{w}^2 A(q) \left[\left(m_\pi ^2+2q^2\right)\left(2m_\pi +\tilde{w}^2 A(q)\right) +4g_A^2 m_\pi \tilde{w}^2\right], \end{eqnarray}$$

(B31),(B32),(B33)

as well as those with the isospin-dependent terms,

$$\begin{eqnarray} W_C^{(4)}(q,Q) &=& W_C^{(c_i/m_N)}(q,Q) +W_C^{(m_N^{-2})}(q,Q) +W_C^{(\mathrm{2L})}(q,Q),\\ W_C^{(c_i/m_N)}(q,Q) &=& -\frac{1}{48(2\pi )^5 m_N f_\pi ^4} c_4 q^2 L(q) \left[g_A^2\left(8m_\pi ^2+5q^2\right)+w^2\right], \\ W_C^{(m_N^{-2})}(q,Q) &=& -\frac{1}{6(4\pi )^5 m_N^2 f_\pi ^4} L(q) \\ && \times \Bigg [ \bigg [ 8g_A^2\left\lbrace \frac{3}{2}q^4+3m_\pi ^2q^2+3m_\pi ^4-6m_\pi ^6w^{-2} -Q^2\left(8m_\pi ^2+5q^2\right)\right\rbrace \\ && +\,4g_A^4\left\lbrace Q^2\left(20m_\pi ^2+7q^2-16m_\pi ^4w^{-2}\right) \right. \\ && \left. +16m_\pi ^8w^{-4}+12m_\pi ^6w^{-2} -4m_\pi ^4q^2w^{-2}-5q^4-6m_\pi ^2q^2-6m_\pi ^4 \right\rbrace -4Q^2w^2 \bigg ] \\ && +\,16g_A^4m_\pi ^6w^{-2} \Bigg ], \\ W_C^{(\mathrm{2L})}(q,Q) &=& \frac{1}{9(4\pi )^7 f_\pi ^6} L(q) \\ && \times \Bigg [ 192\pi ^2f_\pi ^2w^2\bar{d}_3 \left\lbrace 2g_A^2\tilde{w}^2-\frac{3}{5}\left(g_A^2-1\right) w^2\right\rbrace +\left\lbrace 6g_A^2\tilde{w}^2-\left(g_A^2-1\right)w^2\right\rbrace \\ && \times \bigg [ 384\pi ^2f_\pi ^2 \left\lbrace \tilde{w}^2 \left(\bar{d}_1 +\bar{d}_2\right)+4m_\pi ^2 \bar{d}_5 \right\rbrace +L(q)\left\lbrace 4m_\pi ^2\left(1+2g_A^2\right)+q^2\left(1+5g_A^2\right)\right\rbrace \\ && \left. -\frac{1}{3}q^2\left(5+13g_A^2\right) -8m_\pi ^2\left(1+2g_A^2\right) \right]\Bigg ]. \end{eqnarray}$$

(B34),(B35),(B36),(B37)

Here, |$\bar{d}_i$| are the LECs. The spin-spin terms are given by

$$\begin{eqnarray} V_S^{(4)}(q,Q) &=& V_S^{(m_N^{-2})}(q,Q) +V_S^{(\mathrm{2L})}(q,Q), \\ V_S^{(m_N^{-2})}(q,Q) &=& -\frac{g_A^4}{8(2\pi )^5 m_N^2 f_\pi ^4} q^2 L(q) \left[Q^2+\frac{5}{8}q^2+m_\pi ^4w^{-2}\right], \\ V_S^{(\mathrm{2L})}(q,Q) &=& \frac{g_A^2}{8(2\pi )^5 f_\pi ^4} \left(\bar{d}_{14}-\bar{d}_{15}\right) q^2 w^2 L(q), \end{eqnarray}$$

(B38),(B39),(B40)

and

$$\begin{eqnarray} W_S^{(4)}(q,Q) &=& V_S^{(c_i^2)}(q,Q) +V_S^{(c_i/m_N)}(q,Q) +V_S^{(m_N^{-2})}(q,Q) +V_S^{(\mathrm{2L})}(q,Q),\\ W_S^{(c_i^2)}(q,Q) &=& -\frac{4}{3(4\pi )^5 f_\pi ^4} c_4^2 q^2 w^2 L(q), \\ W_S^{(c_i/m_N)}(q,Q) &=& \frac{1}{48(2\pi )^5 m_N f_\pi ^4} c_4 q^2 L(q) \left[g_A^2\left(16m_\pi ^2+7q^2\right)-w^2\right],\\ W_S^{(m_N^{-2})}(q,Q) &=& -\frac{1}{12(4\pi )^5 m_N^2 f_\pi ^4} q^2 L(q) \\ && \times \left[ 4g_A^4\left(7m_\pi ^2+\frac{17}{4}q^2+4m_\pi ^4w^{-2}\right) -32g_A^2\left(m_\pi ^2+\frac{7}{16}q^2\right)+w^2 \right],\\ W_S^{(\mathrm{2L})}(q,Q) &=& -\frac{g_A^4}{16(4\pi )^5 f_\pi ^6} q^2 w^2 A(q) \left[ w^2 A(q) +2m_\pi \left(1+2g_A^2\right) \right]. \end{eqnarray}$$

(B41),(B42),(B43),(B44),(B45)

The SO terms read

$$\begin{eqnarray} V_{LS}^{(4)}(q,Q) &=& V_{LS}^{(c_i/m_N)}(q,Q) +V_{LS}^{(m_N^{-2})}(q,Q),\\ V_{LS}^{(c_i/m_N)}(q,Q) &=& \frac{g_A^2}{2(2\pi )^5 m_N f_\pi ^4} c_2 w^2 L(q), \\ V_{LS}^{(m_N^{-2})}(q,Q) &=& \frac{g_A^4}{(2\pi )^5 m_N^2 f_\pi ^4} L(q) \left[\frac{11}{32}q^2+m_\pi ^4w^{-2}\right], \end{eqnarray}$$

(B46),(B47),(B48)

as well as

$$\begin{eqnarray} W_{LS}^{(4)}(q,Q) &=& W_{LS}^{(c_i/m_N)}(q,Q) +W_{LS}^{(m_N^{-2})}(q,Q),\\ W_{LS}^{(c_i/m_N)}(q,Q) &=& -\frac{1}{12(2\pi )^5 m_N f_\pi ^4} c_4 L(q) \left[g_A^2\left(8m_\pi ^2+5q^2\right)+w^2\right], \\ W_{LS}^{(m_N^{-2})}(q,Q) &=& \frac{1}{2(4\pi )^5 m_N^2 f_\pi ^4} L(q) \\ && \times \left[ 16g_A^2\left(m_\pi ^2+\frac{3}{8}q^2\right) +\frac{4}{3}g_A^4\left(4m_\pi ^4w^{-2}-\frac{11}{4}q^2-9m_\pi ^2\right) -w^2\right]. \end{eqnarray}$$

(B49),(B50),(B51)

The tensor terms are expressed by

$$\begin{eqnarray} V_T^{(4)}(q,Q) &=& V_T^{(m_N^{-2})}(q,Q) +V_T^{(\mathrm{2L})}(q,Q),\\ V_T^{(X)}(q,Q) &=&-\frac{1}{q^2} V_S^{(X)}(q,Q), \end{eqnarray}$$

(B52),(B53)

and

$$\begin{eqnarray} W_T^{(4)}(q,Q) &=& W_T^{(c_i^2)}(q,Q) +W_T^{(c_i/m_N)}(q,Q) +W_T^{(m_N^{-2})}(q,Q) +W_T^{(\mathrm{2L})}(q,Q),\\ W_T^{(X)}(q,Q) &=& -\frac{1}{q^2} W_S^{(X)}(q,Q). \end{eqnarray}$$

(B54),(B55)

Here, X is a representative of |$c_i^2$|⁠, c_i/m_N, |$m_N^{-2}$|⁠, and 2L. The |$V_{\sigma L}^{(4)}$| term in Eq. (B27) also behaves as a tensor force, and the potential is given by

$$\begin{eqnarray} V_{\sigma L}^{(4)}(q,Q) = V_{\sigma L}^{(m_N^{-2})}(q,Q) = \frac{g_A^4}{8(2\pi )^5 m_N^2 f_\pi ^4} L(q). \end{eqnarray}$$

(B56)

Appendix C. MPE function

The purpose of the MPE is to express the chiral potential, which is originally given as a function of |$\boldsymbol{q}$| and |$\boldsymbol{Q}$|⁠, in terms of |$\boldsymbol{p}$| and |$\boldsymbol{p}^{\prime }$|⁠. This is relevant to the nonlocal regularization.

In this section, we derive the explicit form of the MPE functions, |${ f}_{L}^{(C)}$|⁠, |${ f}_{LL^{\prime }L_1L_2K}^{(LS)}$|⁠, |${ f}_{\lambda _0 K}^{(T)}$|⁠, and |${ f}_{L_qL_QK}^{(\sigma L)}$|⁠. In general, the MPE function, which depends on p and p^′, is given by integration over x, where x is defined by

$$\begin{eqnarray} x=\frac{\boldsymbol{p}\cdot \boldsymbol{p}^{\prime }}{pp^{\prime }}. \end{eqnarray}$$

(C1)

It appears in q and Q as |$q=\sqrt{p^2+p^{\prime 2}-2pp^{\prime }x}$| and |$Q=\sqrt{p^2+p^{\prime 2}+2pp^{\prime }x}/2$|⁠, respectively. This integration needs to be performed numerically for the 1π and 2π terms, whereas it can be calculated analytically for the contact terms.

C.1. Central contributions

The central contributions are in association with the operators, |$\mathbb {1}$|⁠, |$\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2$|⁠, |$\boldsymbol{\tau }_1\cdot \boldsymbol{\tau }_2$|⁠, and |$\left(\boldsymbol{\sigma }_1\cdot \boldsymbol{\sigma }_2\right)\left(\boldsymbol{\tau }_2\cdot \boldsymbol{\tau }_2\right)$|⁠. The corresponding terms of the chiral interaction can be found in Table 1. It is convenient to separate |${ f}_{L}^{(C)}$| of the contact terms from that of the 2π terms. Furthermore, we put an additional symbol explicitly in the superscript of the MPE function to distinguish each term. Thus, we find |${ f}_{L}^{(C)}$| of the contact terms as

$$\begin{eqnarray} { f}_{L}^{(C;C_S)}\!\left(p,p^{\prime }\right) &=& \frac{C_S}{C_T} { f}_{L}^{(C;C_T)}\!\left(p,p^{\prime }\right) = \frac{2}{\pi }C_S \delta _{L0},\\ { f}_{L}^{(C;C_1)}\!\left(p,p^{\prime }\right) &=& \frac{4C_1}{C_2}(-)^L { f}_{L}^{(C;C_2)}\!\left(p,p^{\prime }\right) = \frac{C_1}{C_3} { f}_{L}^{(C;C_3)}\!\left(p,p^{\prime }\right) = \frac{4C_1}{C_4}(-)^L {f}_{L}^{(C;C_4)}\!\left(p,p^{\prime }\right) \\ &=& \frac{2}{\pi } C_1\left[\left(p^2+p^{\prime 2}\right)\delta _{L 0}-2pp^{\prime }\delta _{L 1}\right], \\ { f}_{L}^{(C;D_1)}\!\left(p,p^{\prime }\right) &=& \frac{16D_1}{D_2}(-)^L { f}_{L}^{(C;D_2)}\!\left(p,p^{\prime }\right) = \frac{D_1}{D_5} { f}_{L}^{(C;D_5)}\!\left(p,p^{\prime }\right) = \frac{16D_1}{D_6}(-)^L { f}_{L}^{(C;D_6)}\!\left(p,p^{\prime }\right) \\ &=& \frac{2}{\pi }D_1\left[\left\lbrace \left(p^2+p^{\prime 2}\right)^2+\frac{4}{3}p^2p^{\prime 2}\right\rbrace \delta _{L 0} -4 \left(p^2+p^{\prime 2}\right)pp^{\prime }\delta _{L 1} +\frac{8}{3} p^2p^{\prime 2}\delta _{L 2}\right],\\ { f}_{L}^{(C;D_3)}\!\left(p,p^{\prime }\right) &=& \frac{D_3}{D_7} { f}_{L}^{(C;D_7)}\!\left(p,p^{\prime }\right) \\ &=& \frac{1}{2\pi }D_3 \left[\left\lbrace \left(p^2+p^{\prime 2}\right)^2 -\frac{4}{3}p^2p^{\prime 2}\right\rbrace \delta _{L 0} -\frac{8}{3}p^2p^{\prime 2}\delta _{L 2} \right],\\ { f}_{L}^{(C;D_4)}\!\left(p,p^{\prime }\right) &=& \frac{2D_4}{3D_3} { f}_{L}^{(C;D_3)}\!\left(p,p^{\prime }\right) +D_4 \mathcal {G}_{L}\!\left(p,p^{\prime }\right), \\ { f}_{L}^{(C;D_8)}\!\left(p,p^{\prime }\right) &=& \frac{2D_8}{3D_7} { f}_{L}^{(C;D_7)}\!\left(p,p^{\prime }\right) +D_8 \mathcal {G}_{L}\!\left(p,p^{\prime }\right),\\ \mathcal {G}_{L}\!\left(p,p^{\prime }\right) &=& -\frac{1}{15\pi } \sum _{\Lambda _q=0}^2 \sum _{\Lambda _Q=0}^2 (-)^{\Lambda _Q} \widehat{2\!-\!\Lambda _q}\widehat{2\!-\!\Lambda _Q} \left[\binom{5}{2\Lambda _q}\binom{5}{2\Lambda _Q}\right]^{\!\frac{1}{2}} \\ && \times \left( \Lambda _q 0 \Lambda _Q 0 | L 0 \right) \left( 2\!-\!\Lambda _q, 0, 2\!-\!\Lambda _Q, 0 | L 0 \right) {\left\lbrace \begin{array}{*{10}c}\Lambda _q && 2\!-\!\Lambda _q && 2 \\ 2\!-\!\Lambda _Q && \Lambda _Q && L \end{array}\right\rbrace } \\ && \times \, p^{\Lambda _q+\Lambda _Q} p^{\prime 4-\Lambda _q-\Lambda _Q}. \end{eqnarray}$$

(C2),(C3),(C4),(C5),(C6),(C7),(C8)

As regards the pion-exchange terms, |${ f}_{L}^{(C)}$| can be always written as

$$\begin{eqnarray} { f}_{L}^{(C)}\!\left(p,p^{\prime }\right) &= \frac{(4\pi )^2}{2}\hat{L}^2 \int _{-1}^1 dx P_L(x) G^{(C)}(q,Q), \end{eqnarray}$$

(C9)

where G^(C)(q, Q) is a representative for each term defined in Appendix  B:

$$\begin{eqnarray} G^{(C)}(q,Q) = \left\lbrace \begin{aligned} &W_C^{(2)}(q,Q) &&\qquad (\text{$\tau $ exchange at NLO}),\\ &V_S^{(2)}(q,Q) &&\qquad (\text{$\sigma $ exchange at NLO}),\\ &V_C^{(3)}(q,Q) &&\qquad (\text{Purely central at N$^2$LO}),\\ &V_S^{(3)}(q,Q) &&\qquad (\text{$\sigma $ exchange at N$^2$LO}),\\ &W_C^{(3)}(q,Q) &&\qquad (\text{$\tau $ exchange at N$^2$LO}),\\ &W_S^{(3)}(q,Q) &&\qquad (\text{$\sigma $/$\tau $ exchange at N$^2$LO}),\\ &V_C^{(c_i^2)}(q,Q) &&\qquad (\text{Purely central for N$^3$LO-$c_i^2$ term}),\\ &V_C^{(c_i/m_N)}(q,Q) &&\qquad (\text{Purely central for N$^3$LO-$c_i/m_N$ term}),\\ &V_C^{(m_N^{-2})}(q,Q) &&\qquad (\text{Purely central for N$^3$LO-$m_N^{-2}$ term}),\\ &V_C^{(\mathrm{2L})}(q,Q) &&\qquad (\text{Purely central for N$^3$LO-2L term}),\\ &W_C^{(c_i/m_N)}(q,Q) &&\qquad (\text{$\tau $ exchange for N$^3$LO-$c_i/m_N$ term}),\\ &W_C^{(m_N^{-2})}(q,Q) &&\qquad (\text{$\tau $ exchange for N$^3$LO-$m_N^{-2}$ term}),\\ &W_C^{(\mathrm{2L})}(q,Q) &&\qquad (\text{$\tau $ exchange for N$^3$LO-2L term}),\\ &V_S^{(m_N^{-2})}(q,Q) &&\qquad (\text{$\sigma $ exchange for N$^3$LO-$m_N^{-2}$ term}),\\ &V_S^{(\mathrm{2L})}(q,Q) &&\qquad (\text{$\sigma $ exchange for N$^3$LO-2L term}),\\ &W_S^{(c_i^2)}(q,Q) &&\qquad (\text{$\sigma $/$\tau $ exchange for N$^3$LO-$c_i^2$ term}),\\ &W_S^{(c_i/m_N)}(q,Q) &&\qquad (\text{$\sigma $/$\tau $ exchange for N$^3$LO-$c_i/m_N$ term}),\\ &W_S^{(m_N^{-2})}(q,Q) &&\qquad (\text{$\sigma $/$\tau $ exchange for N$^3$LO-$m_N^{-2}$ term}),\\ &W_S^{(\mathrm{2L})}(q,Q) &&\qquad (\text{$\sigma $/$\tau $ exchange for N$^3$LO-2L term}). \end{aligned} \right. \end{eqnarray}$$

(C10)

C.2. SO contributions

Here again, the contact and 2π terms of the SO contributions are separately formulated. The MPE function of the SO contributions of the contact terms reads

$$\begin{eqnarray} { f}_{LL^{\prime }L_1L_2K}^{(LS;C_5)}\!\left(p,p^{\prime }\right) &=& 4\sqrt{3} C_5 {\left\lbrace \begin{array}{*{10}c}\lambda _q & 1\!-\!\lambda _q & 1 \\ \lambda _Q & 1\!-\!\lambda _Q & 1 \\ L & L^{\prime } & 1 \end{array}\right\rbrace } \delta _{LL_1} \delta _{L^{\prime }L_2} \delta _{K1}, \\ { f}_{LL^{\prime }L_1L_2K}^{(LS;D_9)}\!\left(p,p^{\prime }\right) &=& \frac{4D_9}{D_{10}}(-)^K { f}_{LL^{\prime }L_1L_2K}^{(LS;D_{10})}\!\left(p,p^{\prime }\right) \\ &=& -\,4\sqrt{3} D_9 \hat{L}_1\hat{L}_2 \left( L_1 0 K 0 | L 0 \right) \left( L_2 0 K 0 | L^{\prime } 0 \right) \\ && \times {\left\lbrace \begin{array}{*{10}c}L && L^{\prime } && 1 \\ L_2 && L_1 && K \end{array}\right\rbrace } {\left\lbrace \begin{array}{*{10}c}\lambda _q && 1\!-\!\lambda _q && 1 \\ \lambda _Q && 1\!-\!\lambda _Q && 1 \\ L_1 && L_2 && 1 \end{array}\right\rbrace } \left[ \left(p^2+p^{\prime 2}\right)\delta _{K 0} -2pp^{\prime }\delta _{K 1}\right], \end{eqnarray}$$

(C11),(C12)

where we put the corresponding LECs in the superscripts.

The 2π-SO contributions are characterized by

$$\begin{eqnarray} { f}_{LL^{\prime }L_1L_2K}^{(LS)}\!\left(p,p^{\prime }\right) &=& -\frac{\sqrt{3}(4\pi )^3}{4} \hat{L}^2 \left( L_1 0 K 0 | L 0 \right) \left( L_2 0 K 0 | L^{\prime } 0 \right) \\ &&\times {\left\lbrace \begin{array}{*{10}c}L && L^{\prime } && 1 \\ L_2 && L_1 && K \end{array}\right\rbrace } {\left\lbrace \begin{array}{*{10}c}\lambda _q && 1\!-\!\lambda _q && 1 \\ \lambda _Q && 1\!-\!\lambda _Q && 1 \\ L_1 && L_2 && K \end{array}\right\rbrace } \\ &&\times \int _{-1}^1 dx P_K(x) G^{(LS)}(q,Q), \end{eqnarray}$$

(C13)

and G^(LS) is given by

$$\begin{eqnarray} G^{(LS)}(q,Q) = \left\lbrace \begin{aligned} &V_{LS}^{(3)}(q,Q) &&\qquad (\text{SO at N$^2$LO}),\\ &W_{LS}^{(3)}(q,Q) &&\qquad (\text{SO with $\tau $ exchange at N$^2$LO}),\\ &V_{LS}^{(c_i/m_N)}(q,Q) &&\qquad (\text{SO for N$^3$LO-$c_i/m_N$ term}),\\ &V_{LS}^{(m_N^{-2})}(q,Q) &&\qquad (\text{SO for N$^3$LO-$m_N^{-2}$ term}),\\ &W_{LS}^{(c_i/m_N)}(q,Q) &&\qquad (\text{SO with $\tau $ exchange for N$^3$LO-$c_i/m_N$ term}),\\ &W_{LS}^{(m_N^{-2})}(q,Q) &&\qquad (\text{SO with $\tau $ exchange for N$^3$LO-$m_N^{-2}$ term}). \end{aligned} \right. \end{eqnarray}$$

(C14)

See Appendix  B for the explicit form of G^(LS).

C.3. Tensor contributions

By explicitly putting LECs in the superscript of |${ f}_{\lambda _0K}^{(T)}$| to distinguish each term of the chiral interaction, we can write the MPE function for the tensor contributions of the contact terms:

$$\begin{eqnarray} { f}_{\lambda _0K}^{(T;C_6)}\!\left(p,p^{\prime }\right) &=& \frac{4C_6}{C_7} (-)^K { f}_{\lambda _0K}^{(T;C_7)}\!\left(p,p^{\prime }\right) \\ &=& 24C_6 \left[ \left\lbrace \left(p^2+p^{\prime 2}\right)\delta _{K 0}-2pp^{\prime }\delta _{K 1} \right\rbrace \delta _{\lambda _0 0} + \delta _{\lambda _0 2}\delta _{K 0} \right],\\ { f}_{\lambda _0K}^{(T;D_{11})}\!\left(p,p^{\prime }\right) &=& 24D_{11} \Bigg [ \left[\left\lbrace \left(p^2+p^{\prime 2}\right)^2 +\frac{4}{3}p^2p^{\prime 2}\right\rbrace \delta _{K0} -4\left(p^2+p^{\prime 2}\right)pp^{\prime }\delta _{K1} +\frac{8}{3}p^2p^{\prime 2}\delta _{K2} \right]\delta _{\lambda _0 0} \\ &&+\left[ \left(p^2+p^{\prime 2}\right)\delta _{K0}-2pp^{\prime }\delta _{K1} \right]\delta _{\lambda _0 2} \Bigg ], \\ { f}_{\lambda _0K}^{(T;D_{12})}\!\left(p,p^{\prime }\right) &=& 6D_{12} \Bigg [ \left[\left\lbrace \left(p^2+p^{\prime 2}\right)^2 -\frac{4}{3}p^2p^{\prime 2}\right\rbrace \delta _{K0} -\frac{8}{3}p^2p^{\prime 2}\delta _{K2} \right]\delta _{\lambda _0 0} \\ &&+\left[ \left(p^2+p^{\prime 2}\right)\delta _{K0}+2pp^{\prime }\delta _{K1} \right]\delta _{\lambda _0 2} \Bigg ],\\ { f}_{\lambda _0K}^{(T;D_{13})}\!\left(p,p^{\prime }\right) &=& D_{13} \left[ \frac{4}{D_{12}}{ f}_{\lambda _0K}^{(T;D_{12})}\!\left(p,p^{\prime }\right)\delta _{\lambda _0 0} + \frac{1}{D_{11}}{ f}_{\lambda _0K}^{(T;D_{11})}\!\left(p,p^{\prime }\right)\delta _{\lambda _0 2} \right],\\ { f}_{\lambda _0K}^{(T;D_{14})}\!\left(p,p^{\prime }\right) &=& D_{14} \left[ \frac{1}{16D_{11}}{ f}_{\lambda _0K}^{(T;D_{11})}\!\left(p,p^{\prime }\right)\delta _{\lambda _0 2} + \frac{4}{D_{12}}{ f}_{\lambda _0K}^{(T;D_{12})}\!\left(p,p^{\prime }\right)\delta _{\lambda _0 0} \right]. \end{eqnarray}$$

(C15),(C16),(C17),(C18),(C19)

A form similar to Eq. (C9) can be found for |${ f}_{\lambda _0K}^{(T)}$| of the pion-exchange terms as

$$\begin{eqnarray} { f}_{\lambda _0K}^{(T)}\!\left(p,p^{\prime }\right) &= \frac{3(4\pi )^3}{2}\hat{L}^2 \int _{-1}^1 dx P_K(x) G^{(T)}(q,Q), \end{eqnarray}$$

(C20)

with

$$\begin{eqnarray} G^{(T)}(q,Q) = \left\lbrace \begin{aligned} &-\frac{1}{(2\pi )^3}\!\left(\frac{g_A}{2f_\pi }\right)^{\!2} &&\qquad (\text{$1\pi $}),\\ &V_{T}^{(2)}(q,Q) &&\qquad (\text{Tensor at NLO}),\\ &V_{T}^{(3)}(q,Q) &&\qquad (\text{Tensor at N$^2$LO}),\\ &W_{T}^{(3)}(q,Q) &&\qquad (\text{Tensor with $\tau $ exchange at N$^2$LO}),\\ &V_T^{(m_N^{-2})}(q,Q) &&\qquad (\text{Tensor for N$^3$LO-$m_N^{-2}$ term}),\\ &V_T^{(\mathrm{2L})}(q,Q) &&\qquad (\text{Tensor for N$^3$LO-2L term}),\\ &W_T^{(c_i^2)}(q,Q) &&\qquad (\text{Tensor with $\tau $ exchange for N$^3$LO-$c_i^2$ term}),\\ &W_T^{(c_i/m_N)}(q,Q) &&\qquad (\text{Tensor with $\tau $ exchange for N$^3$LO-$c_i/m_N$ term}),\\ &W_T^{(m_N^{-2})}(q,Q) &&\qquad (\text{Tensor with $\tau $ exchange for N$^3$LO-$m_N^{-2}$ term}),\\ &W_T^{(\mathrm{2L})}(q,Q) &&\qquad (\text{Tensor with $\tau $ exchange for N$^3$LO-2L term}). \end{aligned} \right. \end{eqnarray}$$

(C21)

Again, the explicit forms of G^(T) can be found in Appendix  B.

The σL terms in association with D₁₅ and |$V_{\sigma L}^{(m_N^{-2})}$| are also the tensor contributions. The former has the MPE function,

$$\begin{eqnarray} { f}_{L_qL_QK}^{(\sigma L;D_{15})}\!\left(p,p^{\prime }\right) &=& -\,864 D_{15} \Bigg [\frac{1}{4}\left[\left\lbrace \left(p^2+p^{\prime 2}\right)^2 -\frac{4}{3}p^2p^{\prime 2}\right\rbrace \delta _{K0} -\frac{8}{3}p^2p^{\prime 2}\delta _{K2} \right]\delta _{L_q 0}\delta _{L_Q 0} \\ &&+\,\left[ \left(p^2+p^{\prime 2}\right)\delta _{K0}-2pp^{\prime }\delta _{K1} \right]\delta _{L_q 0}\delta _{L_Q 2} \\ &&+\,\frac{1}{4}\left[ \left(p^2+p^{\prime 2}\right)\delta _{K0}+2pp^{\prime }\delta _{K1} \right]\delta _{L_q 2}\delta _{L_Q 0} +\delta _{K0}\delta _{L_q 2}\delta _{L_Q 2} \Bigg ], \end{eqnarray}$$

(C22)

while that for the latter reads

$$\begin{eqnarray} { f}_{L_qL_QK}^{(\sigma L;m_N^{-2})}\!\left(p,p^{\prime }\right) &= -2(12\pi )^3\hat{K}^2 \int _{-1}^1 dx P_K(x) q^{2-L_q} Q^{2-L_Q} V_{\sigma L}^{(m_N^{-2})}(q,Q), \end{eqnarray}$$

(C23)

where |$V_{\sigma L}^{(m_N^{-2})}$| is defined by Eq. (B56).