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

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 \textit{ab initio} description of various cluster phenomena in view of nucleons, pions, and many-nucleon forces.


Introduction
So far, an enormous number of studies have been carried out to understand various phenomena in nuclear systems from the view point of nuclear force. For this purpose, many-body methods based on the local Gaussian basis functions have been established and widely applied; for example, the fermionic molecular dynamics (FMD) [1,2] and antisymmetrized molecular dynamics (AMD) [3][4][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][7][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 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 Sec. 2, we give the formalism of the two-body MEs in momentum space. In Sec. 3, as an example, typical MEs are visualized. Then, Sec. 4 is devoted to a summary and future perspectives. Appendix A briefly shows how the MEs formulated in this work enter the practical many-body calculations. The detail of the chiral potentials and information relevant to the formalism are relegated to Appendices B and C.

Wave functions
We compute the two-body MEs in momentum space for the chiral interaction as the twonucleon force, while other terms of the Hamiltonian can be calculated analytically (see Ref. [26] for example). Since the main focus of this paper is the derivation of the twobody MEs of the chiral interaction by local Gaussians, we express the single-particle wave function in momentum space: where ϕ (ν) i is the Fourier transform of its coordinate-space representation: i (p n ) = 1 (2π) 3 2 2ν π 3 4 dr exp(−ip n · r) exp −ν (r − Z i ) 2 , with 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 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 Z i corresponds to the expectation value of the nucleon momentum. Indeed, if Z i is pure imaginary in Eq. (2), ϕ (ν) i localizes 3/39 at around p n ∼ 2νIm(Z i ). Therefore, complex 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 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 which can be rewritten in terms of the center-of-mass (CM) momentum P and the relative momentum p as with Here, neutrons and protons are assumed to have common mass.

General form of two-body matrix elements
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, p and p ′ , 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 loss of generality of the formalism. Thus, we derive the MEs based on the single-particle state ϕ is defined by Eq. (2) and |ξ i ⟩ is now given by Hence, the spin-up and spin-down states are respectively expressed by |m σi ⟩ = |↑ σ ⟩ and |m σi ⟩ = |↓ σ ⟩. Now the two-body states are represented by with which the MEs are formulated. 4/39 Let us move to the general form of the two-body MEs. A general two-body operatorV 2N can be written aŝ The Galilean invariance requires that the potential v 2N in association with the operator V 2N does not depend on P and P ′ . Hence, such a potential is given by which results in the two-body MEs: with Using Eq. (12), we introduce the notation of the antisymmetrized two-body MEs, For later convenience, we perform the Rayleigh expansion for the plane waves in Eq. (12): where 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 = √ z ij · z ij and z * kl = z * kl · z * kl . The spherical harmonics Y LM with the complex numbers as its arguments is given by the analytic continuation in the Cartesian 5/39 representation: The upper limit of the summation is expressed by the floor function and z = √ z · z = z 2 x + z 2 y + z 2 z is a complex number.

Two-body matrix elements of chiral interaction
Overview. Now we formulate the two-body MEs of the chiral interaction at N 3 LO. One finds that the chiral-N 3 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 , where the central, SO, tensor, and σL contributions are denoted by V , respectively. Although the σL term, the potential of which depends on V (4) σL and D 15 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 Sec. 2.3.5.

Central contributions.
As an example of the central contributions, we take C 1 term, one of the contact terms at NLO. As shown by Eq. (B9), with the nonlocal regularization, its potential is given by where C 1 is one of the low-energy constants (LECs). The transferred momentum q, as well as the average momentum Q are defined by 6/39 with the matrix U given by Eq. (7). The regulator u n characterized by the power n and the cutoff momentum Λ is nonlocal: See Appendix B for more detail of the potential form. To compute the MEs, we perform the multipole expansion (MPE) of q 2 in Eq. (19): and the MPE function is given as The bipolar spherical harmonics is defined by We adopt the abbreviationL = √ 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 plug Eq. (19) into Eq. (12), as well as employing Eqs. (16) and (17), one finds the antisymmetrized MEs as Here, the Legendre polynomial P L has the argument given by and 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 Y K (p) ⊗ Y K p ′ 00 as 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 pionexchange terms, the MPE of the q-and Q-dependent parts of the potentials involving the pion propagators is necessary, but after some manipulation, Y K (p) ⊗ Y K p ′ 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 1 /π and F (C1) L with appropriate forms, as explained in Sec. 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).

2.3.3.
Spin-orbit contributions. The C 5 term, one of the contact terms at NLO is a good example of the SO contributions. Its potential reads where C 5 is the LEC. The total spin S is given by The spin operator σ i is represented by the Pauli matrices. The operator −iS · (q × Q) can be represented in terms of p and p ′ as where the binomial coefficient is given by and the 9-j symbol is represented by the 3 × 3 matrix in the braces. 8/39 From Eqs. (12), (28), and (30), the antisymmetrized MEs can be calculated as with The expression, Eq. (30), can be applicable to other spin-orbit terms, i.e., all spin-orbit potentials can be written in terms of Y K1 (p) ⊗ Y K2 p ′ 1 ⊗ S 1 00 even if pions are exchanged. Consequently, for the MEs of the spin-orbit contributions, the structure of Eq. (32) is rather general as shown in Sec. 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 Sec. 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.

Tensor contributions.
The typical tensor contribution appears as the C 6 -contact term at NLO. As shown in Appendix B.3, its potential is written as where C 6 is the LEC. We express the operator in Eq. (34) by the irreducible-tensor representation: with The 2 × 3 matrix in the braces is the 6-j symbol. One finds that Eq. (35) with λ 0 = 0 corresponds to the central component, q 2 (σ 1 · σ 2 ) /3, while that with λ 0 = 2 is the purely tensor component, q 2 S 12 /3, with the tensor operator Note that S is defined by Eq. (29). In this paper, the terms dependent on (σ 1 · q) (σ 2 · q) are referred to as the tensor contributions, although they involve the central contributions.

10/39
From the above expressions, the antisymmetrized MEs are computed as with There is no counterpart of T (λ0) ijkl 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 Q instead of q. Every tensor potential can be expressed in terms of λ0 00 with λ 0 = 0 or 2. Therefore, again, the structure of Eq. (38) is essentially same as that of other tensor terms, and differences can be found only in the prefactors and the MPE function F (C6) λ0K . A general form of the tensor MEs is thus obtained in Sec. 2.3.5.
Note that, following the procedure for the C 6 -term MEs, one can derive the MEs of the σL terms in association with the operator, [σ 1 · (q × Q)] [σ 2 · (q × Q)], which appears in Eqs. (B27) and (B28). Indeed, this operator can be written in terms of As a result, their MEs have the structure essentially same as that of Eq. (38), i.e., λ 0 = 0 and 2 respectively correspond to the central and tensor components, although the prefactors and coefficients are much more complicated.
2.3.5. Summary of chiral two-body matrix elements. In the previous sections, the typical MEs of the central, SO, and tensor contributions are shown. Now, in this section, we 11/39 generalize them: When we focus on the central MEs, for example, the correspondence between Eqs. (25) and (40)  , which is newly defined as f (C) L , and by replacing D ijkl,L , one can obtain Eq. (40). The explicit form of the MPE functions, f λ0K , and f (σL) LqLQK are relegated to Appendix C.
Each term of the MEs involves U (X) ijkl,L , which depends on the operator form as summarized in Table 1. For instance, we explicitly show the correspondence between U (X) ijkl,L and the spin-isospin MEs of the C 1 , C 5 , and C 6 terms derived in the previous sections. They are characterized by the operators, 1, −iS · (q × Q), and (σ 1 · q) (σ 2 · q): If the interaction involves the spin-spin operator, 13/39 Table 1 The spin-isospin MEs U (X) ijkl,L of the chiral interaction. The tensor operator can be (σ 1 · q) (σ 2 · q), (σ 1 · Q) (σ 2 · Q), and [σ 1 · (q × Q)] [σ 2 · (q × Q)]. In the right-most column, the corresponding terms are listed. The N 3 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 T , The same is true for the isospin MEs, i.e., D jikl with the operator τ 1 · τ 2 . Here, τ i is the isospin operator represented by the Pauli matrices. The isospin-isospin operator, τ 1 · τ 2 , does not enter the contact terms of the usual chiral EFT since they are formulated with the choice of 1 and σ 1 · σ 2 based on the Fierz rearrangement freedom [27][28][29]. 14/39

Visualization of two-body matrix elements 3.1. Numerical details
To demonstrate the behavior of the two-body MEs formulated in Sec. 2, here we visualize typical values of them. As an example, we choose the LO (LO plus NLO) contributions, the potentials of which are given in Appendix B.2 (Appendix B.3). To visualize the MEs some constrains are necessary. First, the Gaussian-center position Z i is chosen to be real, and we select the diagonal MEs V (νν) ijij . Thus, the MEs can be computed as a function of r = z ij = 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, r = 0, and another spin-up proton moves on the x-z plane, i.e., now the three-dimensional vector r is expressed by 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, V LO and V NLO , to express the MEs V (νν) ijij of the LO and LO-plus-NLO cases, respectively. 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 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). 15/39 as the constants relevant to the 1π term (the pion mass, pion-decay constant, and axial vector coupling constant). We adopt ν = 0.26 fm −2 [25].

Matrix elements at leading order and contact-term contributions
Under the conditions described in the previous section, V LO can be represented in Figs. 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 V LO 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]. Figures 2(a)-(c) show that the spin-aligned-neutron-proton pair feels the largest attraction at 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]. Furthremore, the attractive MEs at r = 0 can be shown analytically. First, we ignore the regulator with Λ → ∞. Indeed, the largest attraction at 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 Note that the LECs at LO has the charge dependence. The superscript np stands for the LECs for the neutron-proton pair. We find that the condition C np S + C np T < 0 is satisfied by the LECs employed here. For example, C np S + C np T = −0.011374 × 10 4 GeV −2 for Λ = 450 MeV [30,31]. Thus we can show that Eq. (48) has a minimum at r = 0.
We should mention a role played by the 1π interaction at the origin. Within the configuration of Fig. 1 at 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 Figs. 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 Ref. [33] for example), 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 shortrange 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. 16 The repulsion by the 1π-contact term is moderate compared to that by the LO-contact terms. Specifically, for Λ = 450 MeV, the absolute vale 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 r = 0 can be basically explained by the contributions from the LO-contact terms.

One-pion-exchange-tensor contributions to matrix elements
By carefully watching Fig. 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 Ref. [33] for example), the OPEP in the coordinate-space representation, which is the Fourier transform of Eq. (B5), involves the tensor operator, The direction of 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 τ 1 · τ 2 appearing in the OPEP is −3 for the isoscalar state at r = 0. For finite r, the 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 (λ 0 = 2) becomes for which we notice the asymmetry with respect to the operation x ↔ z.  At r = 0, even though the radial part of the tensor component of the OPEP diverges (see Ref. [34] for example), the tensor MEs must be zero, since only the s wave contributes 18/39 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 V (T ) LO 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 Figs. 2(d)-(f) become deeper with larger ν. This is consistent with the ν dependence of tensor-force contributions to the 4 He energy reported in Ref. [35]; the energy gain by a tensor force increases when ν is greater than its typical value ν = 0.25 fm −2 for 4 He [13,15], although the total energy is saturated due to the compensation by the energy loss of the kinetic term for a such large ν.
Also one finds from Figs. 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 which can be obtained from the operator, (σ 1 · q) (σ 2 · q) q 2 + m 2 π −1 through the irreducible tensor representation as Eq. (35) with λ 0 = 2. When large q contributes (this is the case for larger Λ), the effect of O 12 is enhanced. This is consistent with the Λ dependence of V (T ) LO .

Matrix elements at next-to-leading order
Now we increase the order of the chiral EFT up to NLO to visualize the two-body MEs, V NLO , within the two-nucleon configuration same as that in the previous section.  Fig. 1. The NLO MEs of the tensor contributions denoted by V (T ) NLO are extracted from V NLO , and depicted in Fig. 3(b). The similar results for Λ = 500 and 600 MeV are shown in Figs. 4 and 5, respectively. Note that even though the C 5 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 to 5. This is because the direction of 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.
As argued in Sec. 3.2, the LO MEs at 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 r = 0. This repulsive nature stems from the high-momentum scattering described by the NLO potentials dependent on 19 NLO localizes on the z (x) axis, and vice versa for V (T ) LO . This is because the NLO potentials do not depend on the isospin-isospin operator, τ 1 · τ 2 , which induces the change of the sign of the MEs, and also because the LECs, C 6 and C 7 are negative in the present parameterization [8,30,31].

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 21/39 like AMD, which can be applicable, for instance, to efficiently describe nuclear cluster structures. We adopt the chiral potentials at N 3 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 threenucleon force into AMD is one of the important tasks, and the extension of the framework to this direction is also on going.

Acknowledgment
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.

A. Slater determinants and their superposition A.1. Slater determinants
A single-Slater determinant for a system of the mass number A is characterized by a single value of ν and a set of the generator coordinate Z n : with the A-body antisymmetrizerÂ A . Then, many-body states are obtained through the generator-coordinate method (GCM) [26,36,37] after the parity and angular-momentum projections: where the projection operatorsP ± andP J M K are given explicitly in the next sections.

22/39
The coefficient c (ν) n is obtained numerically by solving the generalized eigenvalue problem, with the energy eigenvalue E. Since the single-particle wave function expressed by the Gaussian-wave packet, |Ψ ν (Z n )⟩ can be separated into the CM and intrinsic wave functions. Thus, the MEs of the norm and HamiltonianĤ relevant to the intrinsic structure are given by where the CM-wave function Ψ (CM) ν is given by whereT is the kinetic-energy operator expressed by the sum of the one-body kinetic-energy operator andT 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 Note that Eq. (A9) is obtained under the assumption, Z G = 0, for both bra and ket states. The interaction operatorV consists of nuclear and Coulomb parts, and if it involves the chiral two-nucleon force, H (νν ′ ) nm can be written in terms of the two-body MEs given in Sec. 2.

A.2. Parity projection
First, we introduce the parity inversion operatorP π , which inverts the sign of the Gaussian center position of the single-particle states aŝ Note that the subscript n is omitted from Z for simplicity. Thus, the parity projected operator and parity projected states are respectively defined bŷ where α ± is the normalization coefficient in association with the positive (+) or negative (−) parity.

A.3. Angular-momentum projection
The angular-momentum-projected state Ψ J νM K (Z) is defined with the angularmomentum-projection operatorP J M K by 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 Dfunction D J M K and the rotation operatorR, 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 Z without the subscript n.
The operation ofR results in the rotation of the single-particle wave function. As a result, for the spatial part, we just have to rotate Z ij and z ij as 24/39 As regards the rotation of the spin wave function, the coefficients α i and β i in Eq. (3) are replaced by as well as α j and β j to be replaced by α ′ j and β ′ j , respectively.
B. Potentials derived from chiral effective field theory B.1. Overview We adopt a high-precision two-body potential based on the chiral EFT at N 3 LO [8, [10][11][12] as v 2N in Eq. (12). At this order, the potential can be written order by order as The superscript stands for the chiral-expansion power, i.e., n χ = 0, 2, 3, and 4 for leading order (LO), NLO, next-to-next-to-leading order (N 2 LO), and N 3 LO, respectively. In this work, we employ the nonlocal regularization with the regulator Here Λ is the cutoff momentum. Thus, the potential depends on the relative momenta, p and p ′ , of the initial and final channels, respectively, which are related to the average momentum Q and transferred momentum q by where U is given by Eq. (7). Note that every potential v 2N appearing in this paper involves the prefactor 1/(2π) 3 , which originates from the normalization convention, ⟨p | p ′ ⟩ = δ(p − p ′ ). 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 where g A is the axial vector coupling constant, f π is the pion-decay constant, and m π is the average pion mass, as well as the LECs, C S and C T at LO. The spin and isospin operators, σ i and τ i respectively, are represented by the Pauli matrices.

B.3. Next-to-leading order
At NLO, the potential is given by where S = (σ 1 + σ 2 )/2 and C i are the LECs at NLO. The 2π term contains with Here, the potentials involve the central terms, the SO terms, 27/39 and the tensor terms, The LECs c i enter the 2π term at N 2 LO, and m N denotes the average nucleon mass.
B.5. Next-to-next-to-next-to-leading order The LECs of the contact terms at N 3 LO are D i . 28/39 Following Ref. [8], the 2π-exchange potential at N 3 LO is categorized into several terms, i.e., c 2 i , c i /m N , and m −2 N contributions in association with one-loop diagrams, and also twoloop (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 29/39 as well as those with the isospin dependent terms, Here,d i are the LECs. The spin-spin terms are given by 30/39 and W (4) The SO terms read as well as The tensor terms are expressed by 31/39 and W (4) Here X is a representative of c 2 i , c i /m N , m −2 N , and 2L. The V (4) σL term in Eq. (B27) also behaves as a tensor force, and the potential is given by

C. Multipole-expansion function
The purpose of the MPE is to express the chiral potential, which is originally given as a function of q and Q, in terms of p and p ′ . This is relevant to the nonlocal regularization.
In this section, we derive the explicit form of the MPE functions, f λ0K , and f (σL) LqLQK . In general, the MPE function, which depends on p and p ′ , is given by integration over x, where x is defined by It appears in q and Q as q = p 2 + p ′2 − 2pp ′ x and Q = p 2 + p ′2 + 2pp ′ x/2, respectively. This integration needs to be performed numerically for the 1π and 2π terms, while it can be calculated analytically for the contact terms.

C.1. Central contributions
The central contributions are in association with the operators, 1, σ 1 · σ 2 , τ 1 · τ 2 , and (σ 1 · σ 2 ) (τ 2 · τ 2 ). The corresponding terms of the chiral interaction can be found in Table 1. It is convenient to separate f (C) L of the contact terms from that of the 2π terms. Furthermore, we put an additional symbol explicitly in the superscript of the MPE function 32/39 to distinguish each term. Thus, we find f (C) L of the contact terms as

C.2. Spin-orbit 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 f (LS;C5) where we put the corresponding LECs in the superscripts. The 2π-SO contributions are characterized by and G (LS) is given by