NuHamil : A numerical code to generate nuclear two- and three-body matrix elements from chiral effective field theory

The applicability of nuclear ab initio calculations has rapidly extended over the past decades. However, starting research projects is still challenging due to the required numerical expertise in the generation of underlying nuclear interaction matrix elements and many-body calculations. To ease the first issue, in this paper we introduce the numerical code NuHamil to generate the nucleon-nucleon (NN) and three-nucleon (3N) matrix elements expressed in a spherical harmonic-oscillator basis, inputs of many-body calculations. The ground-state energies for the selected doubly closed shell nuclei are calculated with the no-core shell-model (NCSM) and in-medium similarity renormalization group (IMSRG). The code is written in modern Fortran, and OpenMP+MPI hybrid parallelization is available for the 3N matrix-element calculations.


Introduction
The dynamics of the atomic nucleus are governed by the strong interaction, whose fundamental theory is described by quantum chromodynamics (QCD).Since the quarks are tightly confined in a nucleon, it is well established that nuclear Hamiltonians associated with the interactions between nucleons are a good starting point to understand nuclear structure and reactions.The study of nuclear interaction models has a long story starting from the pion-exchange theory [3].The quantitative understanding of nuclear interactions is an open problem.In the past decade, interactions based on chiral effective-field theory (EFT) [4,5] have become the standard starting point of ab initio many-body calculations.Chiral EFT-based interactions have several advantages over other interactions such as AV18 [6] and CD-Bonn potentials [7].For example, a systematic expansion is possible by ordering the diagrams according to the power counting, suggesting the possibility of an uncertainty quantification due to the truncation in the expansion [8][9][10][11].Further, many-nucleon interactions naturally appear at higher order, explaining the hierarchy of many-body terms.In nuclear physics, it is well known that 3N interactions play an important role (see for example Ref. [12] as a recent review).With the progress in nuclear interactions and methodological developments in many-body problems, nuclear ab initio studies are well motivated.Nowadays, the applicability extends over the nuclear chart [13] and recently reached to the heaviest known doubly-magic system 208 Pb [14], and further applications are expected.
To perform ab initio calculations, the matrix elements of nuclear Hamiltonians (and relevant operators) are essential.However, developing a numerically efficient code for the matrix-element generation requires expert knowledge and can be a barrier for those entering the field.The goal of the NuHamil code is to provide a simple way to generate NN and 3N matrix elements expressed in a spherical HO basis, applicable for basis-expansion methods, such as the nocore shell model (NCSM) [15], coupled-cluster method [16], self-consistent Green's function method [17], in-medium similarity renormalization group [18] approach, and manbody perturbation theory [19].The code currently supports the input formats of the open-source BIGSTICK [20] and imsrg++ [21] codes for the NCSM and IMSRG calculations, respectively.
This paper is organized as follows.In Sec. 2, we clarify the input NN and 3N matrix elements for the many-body calculations and briefly show how to compute them with the single-particle product state in the HO basis.Also, we discuss the free-space similarity renormalization group (SRG) prescription to soften the nuclear interactions in Sec. 3. We show some benchmark many-body calculation results with NCSM and IMSRG calculations in Sec. 4. The usage of the code and the conclusion are given in Sec. 5 and Sec.6, respectively.

Matrix elements of Hamiltonian
Here, we review how the NN and 3N matrix elements enter in the many-body problem.Our numerical goal is to solve the non-relativistic many-body Schrödinger equation H|Ψ n = E n |Ψ n , with the intrinsic Hamiltonian with up to 3N terms Here, p i is the momentum vector of ith nucleon, and m is the nucleon mass.The factor (A − 1)/A in the first term and is from the subtraction of the center-of-mass (cm) kinetic term.A is the nucleon number of the system.The terms V NN i j and V 3N i jk are the NN and 3N interactions, respectively.

Second-quantized representation
To proceed with many-body calculations with basisexpansion methods, we begin with the expression by the second quantization.To this end, we first define creation and annihilation operators for a nucleon in an HO orbit p: c † p and c p. The subscript p is a collective index specifying the HO orbit and defined as p = {n p , l p , j p , m p , t z,p }.Here, n p , l p , j p , m p , and t z,p are the nodal quantum number, orbital angular momentum, total angular momentum, z-component of j p , and the z-component of the isospin distinguishing protons and neutrons, respectively.Note that proton (neutron) states are labeled as t z = −1/2 (1/2).In nuclear physics, the use of an HO basis is particularly useful since the coordinate transformation coefficient is well known [22][23][24], as will be shown in later.The creation and annihilation operators satisfy the anticommutation relations The object δ p q is defined by products of Kronecker's delta and is written as δ p q = δ n p n q δ l p l q δ j p j q δ m p m q δ t z,p t z,q . (3) Applying the creation operators to the nucleon vacuum state |0 , one can define antisymmetrized states.For example, one-, two-, and three-nucleon states can be written as Using the creation and annihilation operators, an arbitrary n-body operator O [n] can be expressed as The object O p′ In the same way, the Hamiltonian in Eq. ( 1) can be quantized as using the matrix elements of one-body kinetic term T p′ p, two-body kinetic term T NN p′ q′ p q, NN interaction V NN p′ q′ p q, and 3N interaction V 3N p′ q′ r′ p qr .

J-coupled scheme
The number of matrix elements defined in Eq. ( 8) is greatly reduced by exploiting the rotational symmetry of the Hamiltonian.To introduce a smaller set of matrix elements, we define the J-coupled two-and three-body states using the Clebsch-Gordan coefficient C j 1 j 2 J m 1 m 2 m 1 +m 2 : |pqr : J pq JM = m p m q m r C j p j q J pq m p m q m p +m q C J pq j r J m p +m q m r M | p qr .
Here, p, q, and r are the quantum number set without m, i.e. p = {n p , l p , j p , t z,p }, and an additional Kronecker's delta product is introduced as δ pq = δ n p n q δ l p l q δ j p j q δ t z,p t z,q .
Note that the factor 1/(1 + δ pq ) in the two-body state is for the normalization so that we have pq : On the other hand, such normalization factor is not usually included in the three-body state.For the threebody state, one can define another state by a different angular momentum coupling order, which should be related with the Wigner's 6 j-symbol.In this paper, the first and second indices are always coupled first, and then the third index is coupled.In practical applications, we need permutations of the indices of the states to further reduce the storage requirement.For the two-body state, it is given by |qp : JM = −(−1) j p + j q −J |pq : JM .
Likewise, permutations of the indices in the three-body state are given by |qrp : J qr JM = − J pq (−1) j q + j r +J qr [J pq ][J qr ] × j p j q J pq j r J J qr |pqr : J pq JM , |rpq : J pr JM = − J pq (−1) j p + j q +J pq [J pq ][J pr ] × j p j q J pq J j r J pr |pqr : J pq JM , |qpr : J pq JM = −(−1) j p + j q −J pq |pqr : J pq JM , |rqp : × j p j q J pq j r J J qr |pqr : J pq JM , |prq : J pr JM = − J pq (−1) j q + j r +J pr +J pq [J pq ][J pr ] × j p j q J pq J j r J pr |pqr : J pq JM .
Here, Wigner's 6 j-symbol with the standard notation [25] and [x] = 2x + 1 are introduced.With the J-coupled states, the two-and three-body matrix elements can be introduced as p ′ q ′ : J ′ M ′ |O [2] |pq : JM and p ′ q ′ r ′ : J p ′ q ′ J ′ M ′ |O [3] |pqr : J pq JM , respectively.Because of the rotational invariance of the Hamiltonian, the Hamiltonian matrix is M-independent and diagonal with respect to J and M. Therefore, we introduce a shorthand notation for the matrix elements.
T NN,J p ′ q ′ pq = p ′ q ′ : JM|T NN |pq : JM , Since the uncoupled matrix elements with tilde indices can be computed from the J-coupled matrix elements, only calculating the J-coupled matrix elements is sufficient for many-body calculations.The relation between the uncoupled and J-coupled matrix elements are the following.

Matrix elements of kinetic terms
The one-body kinetic matrix element T p ′ p is given by Note that the one-body kinetic operator takes the tridiagonal form.Also, the two-body kinetic matrix element can be computed through the non-antisymmetrized J-coupled matrix element T NN,J p ′ q ′ pq T NN,J p ′ q ′ pq = (−1) j q ′ + j p +J 2 Am × j p ′ j q ′ J j q j p 1 p ′ ∇ p q ′ ∇ q , with the reduced matrix element of the gradient operator, which is given by Here, the HO length parameter b 2 ≡ /mω is introduced with the HO frequency ω.The antisymmetrized matrix element is obtained as with δ pq = δ n p n q δ l p l q δ j p j q δ t z,p t z,q .
The main tasks remaining are to compute the matrix elements V NN,J p ′ q ′ pq and V 3N,J p ′ q ′ J pq J p ′ q ′ r ′ pqr .

Nucleon-nucleon matrix elements
We begin with the NN matrix element.One might think that the matrix element can be calculated directly from the integral using the single-particle HO wave function.It is actually done in quantum chemistry.However, this would be a computationally expensive task since functional forms of NN interactions are complicated.Instead, the Talmi-Moshinsky transformation is widely used in nuclear physics: The quantum numbers introduced for the transformation N NN cm , L NN cm , n, l, S , and J NN rel are the NN cm radial quantum number, NN cm orbital angular momentum, relative radial quantum number, relative orbital angular momentum, total spin, and total angular momentum of the relative motion, respectively.The transformation coefficient In the above equation, 9 j-symbol is used with the standard notation [25].The symbol NLnl : Λ|n 1 l 1 n 2 l 2 : Λ d is the HO bracket defined with the notation in Ref. [24].The inner summations in Eq. ( 31) can be performed with an efficient matrix multiplication.Note that the antisymmetrization is not taken into account here.However, it is trivial and can be done by multiplying the factor f pq to Eq. ( 32): The NN matrix element in the relative HO basis V S J NN rel n ′ l ′ nl can be obtained through the integral: with the radial HO wave function: The gamma function Γ(x) and associated Laguerre polynomial (p ′ , p) stored as a function of p ′ and p.Some selected interactions are given in the input_nn_files directory.The available NN interactions are LO -N 4 LO with 500 MeV regulator cutoff by Entem-Machleidt-Nosyk [26], N 3 LO with 500 MeV regulator cutoff by Entem-Machleidt [27], N 2 LO opt [28], N 2 LO sat [29], and ∆-full EFT series by the Gothenburg-Oak Ridge collaboration [30].

Three-nucleon matrix elements
For computational reasons, the 3N matrix elements are calculated within the isospin formalism.The matrix elements with the proton-neutron basis can be obtained through the The recoupling should be done in many-body calculations, and the goal here is to obtain the JT -coupled matrix element.Note that recoupling coefficients from isopin structure have to be considered for the permutation of indices, similar to Eqs. ( 13)- (17).
Since the antisymmetrization of the 3N basis is more complicated than that of the NN basis, the 3N matrix element is cumbersome.The antisymmetrized basis is expressed as the linear combination of the nonantisymmetrized basis: where E, i, and J 3N rel are the HO principle quantum number, label distinguishing the states, and total Jacobi angular momentum, respectively.The collective index β, specifies the non-antisymmetrized basis.The quantum numbers with the subscript '12', n 12 , l 12 , s 12 , j 12 , and t 12 are for the relative motion of nucleons 1 and 2, i.e., the nodal, orbital angular momentum, spin, total angular momentum, and total isospin quantum numbers, respectively.Likewise, n 3 , l 3 , and j 3 are the quantum numbers for the nucleon 3 with respect to the cm of the nucleons 1 and 2. Note that the principle quantum number is defined as The coefficient in the linear combination c iβ can be obtained by the diagonalization of the antisymmetrizer [31,32]: with A = (1 + T 13 T 12 + T 12 T 23 − T 12 − T 13 − T 23 )/6 defined with the exchange operator T i j and the eigenvalue A i .The matrix element of the antisymmetrizer is [31,32] The eigenvalue problem can be separated into {E, J 3N rel , T } blocks.Due to the overcompleteness of the nonantisymmetrized basis, the eigenvalue A i is either 0 or 1, and A i for the physical state has to be 1.Therefore, we always see N A ≤ N NA where N A and N NA are the basis numbers in the antisymmetrized and non-antisymmetrized bases within the {E, J 3N rel , T } block, respectively.In the code, all the 3N operators are stored with the antisymmetrized basis rather than the non-antisymmetrized basis, as it is computationally easier to handle.As another option for the 3N antisymmetrization, one may apply the permutator operator to the 3N momentum state as introduced in Ref. [33].
Similarly to the NN case, the 3N matrix elements can be obtained through the three-body Talmi-Moshinsky transformation: where N 3N cm and L 3N cm denote the 3N cm radial quantum number and 3N cm orbital angular momentum, respectively.The matrix element with The 12 j-symbol of the first kind [25] is used.Note that T pqJ pq N NN cm L NN cm n 12 l 12 s 12 j 12 is defined in Eq. ( 32), and one can find a recursive relation for the N-body Talmi-Moshinsky transformation with the 12 j-symbol and HO bracket.
A typical limit of the three-body Talmi-Moshinky transformation is E 3max = 16, where E 3max is defined as max(2n p + l p + 2n q + l q + 2n r + l r ).This limit does not allow us to obtain converged results for heavier systems with A 100.Recently, the limit was extended to E 3max = 28 [34], by only computing the matrix elements relevant to the normalordered two-body (NO2B) approximation, which is widely used in basis expansion methods.For further details about the NO2B matrix elements, see Ref. [34].The code supports this format as well.
The 3N matrix element in the Jacobi HO basis can be obtained from the matrix element expressed with the non-antisymmetrized basis: with Here, the collective index α is introduced as The momentum π 2 is defined as The momentum-space matrix element includes a regulator function: The regulator function takes either a non-local form a local form or a semilocal form [35] that is not supported in the code.The code fully supports the locally regulated matrix elements at N 2 LO in chiral EFT based on Ref. [36].Also, a newly introduced local-non-local regularized form, [37], is supported.For non-local matrix elements, we have tried to implement along Ref. [38].However, we found a numerical instability in the higher angular momentum partial waves, which would be related to the discussion made in Ref. [39].For this reason, non-local matrix elements are not fully supported, and external input files are required.The code has the capability to read the momentum-space matrix element V J 3N rel T χEFT,α ′ α (p ′ , q ′ , p, q)1 from the HDF files by Hebeler et al. [12,40].In the integral (45), the cubic b-spline interpolation is used to capture the oscillating nature of the HO wave functions.

Similarity renormalization group
The momentum scales of chiral EFT interactions are significantly lower than those of the other potential models such as AV18 [6] and CD-Bonn [7].However, the momentum scale is not sufficiently low to obtain converged results in the many-body calculations.To accelerate convergence, one sometimes softens nuclear interactions.Softening procedures are well summarized for example in Ref. [43].Here, we briefly review the widely used similarity renormalization group (SRG) approach.
In the SRG, we consider a unitary transformation depending on a continuous parameter α: The SRG flow equation can be obtained by differentiating both sides: The antihermitian operator η(α) is known as the generator of the flow equation and can be chosen flexibly [44].The most widely used choice is η(α) = [T kin , H(α)] with the kinetic operator T kin , which guarantees the suppression of the coupling between low and high momenta.The flow equation ( 51) is integrated until the coupling is sufficiently suppressed.Since the unitary transformation does not change the eigenvalues, the SRG can be regarded as a reshuffling of the NN, 3N, and many-body sectors.In other words, the many-body interactions are induced by the SRG evolution even if the original interaction includes only NN interactions.
In practical applications, we extract the SRG-evolved interactions with the subtraction method (see Ref. [41] for example).For NN, the Hamiltonian H NN = T NN + V NN is evolved, and the evolved NN interaction is given as V NN (α) = H NN (α) − T kin .Note that T kin and V NN are the NN kinetic and interaction operators, respectively.From the definition of the flow equation, V NN (α) reproduces the twobody observables obtained through the original interaction V NN .However, this is not true for three-and many-body observables due to missing induced many-body forces, and the many-body observables have an artificial α dependence, showing how much the unitarity of the transformation is broken in the many-body space.To obtain more α-independent result, one has to include induced 3N interaction extracted from the 3N evolution.
For 3N, the starting Hamiltonian is H 3N = T 3N + V [3]  NN , where T 3N is the 3N kinetic operator, and V [3]  NN is the NN interaction embedded into the 3N space.Note that one can add an initial 3N term if required.The induced 3N term can be obtained as V 3N,ind (α) = H 3N (α) − T 3N − V [3]  NN (α).Here, V [3]  NN (α) is V NN (α) obtained from the NN evolution and embedded into 3N space.As seen in the NN evolution, the Table 1 The ground-state energies of selected doubly magic nuclei computed with the NCSM and IMSRG (2).The resolution scale λ is related with the end point of the SRG evolution: λ = α −1/4 in units of fm −1 .The N max and e max truncations are employed in the NCSM and IMSRG(2), respectively.The truncation defined by E 3max = max(e 1 + e 2 + e 3 ) is applied for the input 3N matrix elements.The E 3max values with asterisk indicates that the 3N matrix elements are computed within the half-precision floating-point numbers.Also, E 3max ="none" indicates that there are no input 3N matrix elements.The entry ω is the basis frequency used in the many-body calculations, and the numbers in the parentheses are the parent frequency adopted in the frequency conversion for the 3N matrix elements [41].For 4 He, 16 O, and 40 Ca, the 3N SRG evolution was done in the ramp A space defined in Ref. [41].For 132 Sn, the evolution was done in the N max = 48 space only for J 3N rel ≤ 13/2 channels.The star symbol in the E 3max column indicates that the 3N matrix elements are in half-precision floating point numbers.

Nucleus
Interaction λ (fm 3N evolution preserves three-body observables.The same procedure can be applied for many-body terms, and the many-body evolution is needed until the α-dependence of the many-body observables becomes weak enough.However, in practice, even the four-body SRG evolution is too expensive to do due to the resulting basis dimension and the cost of antisymmetrizing the basis.In the NuHamil code, the SRG evolution can be performed in the NN and 3N sectors, the current state-of-the-art.
The unitary transformation can be obtained from the flow equation for the transformation operator: However, a computationally more moderate way is used in practice, and the unitary transformation is obtained as with the eigenstates of the original and evolved Hamiltonians: Note that the relative phase of |ψ k and |ψ k (α) cannot be determined in general, which affect the sign of the matrix element of the transformation operator.Since the SRG transformation does not change the wave function drastically, the relative phase is fixed such that ψ k |ψ k (α) ≥ 0. In the same way as for the Hamiltonian, the induced three-body term of an operator can be computed in the NuHamil code as done in Ref. [45].The end point of the flow equation is usually parametrized by the momentum scale λ = α −1/4 instead of α, and we follow this convention.
In the code, the SRG evolution is done in a relativecoordinate HO space because the consistent evolution of the other operators is straightforward.This means that the evolution is done in the truncated HO space, and the N max truncation is employed in the code.The N max is defined as The UV momentum scale in the employed N max space is roughly estimated as p UV ∼ √ 2N max mω [46], and we expect that the N max should be increased until p UV is sufficiently larger than the cutoff scale of the interaction, typically 500 MeV.Although we can take sufficiently large N max for the NN evolution2 , the 3N N max can be an issue especially in heavy nuclei calculations [34] even if the frequency conversion technique [41] is used.

Many-body results
Here, we show the ground-state energies for the selected doubly magic nuclei computed with the NCSM and twobody approximated IMSRG [IMSRG (2)] as a benchmark.We do not introduce the theoretical details of the manybody calculation methods.The details can be found in Refs.[15,18,47] and references therein.The numerical codes used here are open source; the NCSM and IMSRG are done with the BIGSTICK [20] and imsrg++ [21] codes, respectively.

Program Summary and Specifications
The NuHamil code is written in modern Fortran.It requires a set of libraries, BLAS, LAPACK, GNU scientific library (gsl), zlib, and hdf5.

Installation
The source code can downloaded from GitHub: $ cd $ git clone https://github.com/Takayuki-Miy⌋ agi/NuHamil-public.git֒→ Note that downloading the code in the home directory is not mandatory, but it is recommended.One needs to download the submodules, linear algebra wrapper and b-spline interpolation: $ cd NuHamil $ git submodule init $ git submodule update The compilation can be done with the make command.The default compiler is GCC Fortran.If a user needs to use another compiler, the Makefile has to be edited appropriately.
The symbolic link will be created by the make install command. 3Once the directory is added to PATH, the code is ready to run.

How to run
A job submission can be controlled by a Python script, and some sample scripts are prepared in the exe directory.For example, the NN and 3N matrix elements can be generated with NuHamil_2BME.pyand NuHamil_3BME.py,respectively.A Python script generates the corresponding input file for NuHamil.exe and submits the job.If an user needs to run a job manually, it can be done with $ NuHamil.exeinput.txt The "input.txt" is the input file based on the Fortran namelist functionality, and the file format is given the following.

Major parameters
In the NuHamil code, there are a number of input parameters.Here, we list some of the basic input parameters that the users might need to change depending on their requirements. 4  rank: integer, particle number of the system.
hw: frequency of the HO basis in the unit of MeV.The typical range is 10 hw 40.-hw_target: target frequency for the frequency conversion technique [41].The parameter is valid if rank > 2.
To turn off the frequency conversion, set hw_target = −1.emax: e max = max(2n + l) truncation for the output labframe HO matrix element file.-e2max: e 2max = max(2n p + l p + 2n q + l q ) truncation for the output lab-frame HO matrix element file.It is recommended to use e2max = 2 × emax.-e3max: E 3max = max(2n p + l p + 2n q + l q + 2n r + l r ) truncation for the output 3N lab-frame HO matrix element file.A typical limit is e3max = 16, and it will not work for the larger e3max because of the memory requirements.
If only the matrix elements relevant for the NO2B approximation are needed [34], e3max = 24 would be a typical choice without MPI parallelization.-file_name_nn: file name of the output lab-frame NN HO matrix elements.-file_name_3n: file name of the output lab-frame 3N HO matrix elements.renorm: renormalization method; "bare", "srg", and "Vlowk" are available.Note that the code does not support the 3N evolution for "Vlowk" option. 4Users can contact the author for the parameters not listed here.
-renorm_space2: NN interaction renormalization space; "ho" and "mom" are available.The NN renormalization procedure is done in HO-("ho") or momentum-("mom")space.The default is "ho".-input_nn_file: file name of the NN interaction represented in the relative momentum space.The files are in the input_nn_files directory.-NNInt: name of the NN interaction.
-N2max: N max truncation for the NN system.
-only_no2b_element: If it is set True, only the matrix elements relevant for the NO2B approximation will be computed [34].-jmax3: maximum value of the 3N Jacobi angular momentum taken into account.This has to be an integer and twice the actual angular momentum.-genuine_3bf: set True if the bare 3N interaction needs to be included.If it is set False and renorm="srg", the SRG induced 3N interaction will be computed.49), in units of MeV.

File format of input NN interactions in relative momentum space
As mentioned in Sec.2.4, some selected NN interaction files are prepared in the input_nn_files directory.Furthermore, one can use their own momentum-space NN interaction.The file needs to be written with the binary 5 .The file should begin with listing the following variables:

Number_of_mesh_points J_max Number_of_relative_coordinate_channels
where the 32-bit integers Number_of_mesh_points, J_max, and Number_of_relative_coordinate_channels are the size of momentum mesh points, maximum total angular momentum in the relative coordinate, i.e., max(J NN rel ), and the number of [J NN rel , (−1) l , S , t z,p + t z,q ] combinations written in the file, respectively.Then, one needs to write the momentum mesh points and corresponding weights for a quadrature method, which are number_of_mesh_points-dimensional arrays with the 64-bit float: momentum_mesh_points weights Finally, the momentum-space matrix for each [J NN rel , (−1) l , S , t z,p + t z,q ] block should be written in the following way: Here, the 32-bit integers Angular_momentum, Parity, Spin, Isospin_z_component, and Matrix_dimension correspond to J NN rel , (−1) l , S , t z,p + t z,q , and the size of momentum-space matrix, respectively.Note that t z,p + t z,q can take either of −1 (proton-proton), 0 (proton-neutron), or 1 (neutron-neutron).The Momentum_space_matrix is the flattened ("Matrix_dimension") 2 -dimensional array with the 64-bit float, corresponding to V S J NN rel l ′ l (p ′ , p) 6 .Notice that Number_of_mesh_points and Matrix_dimension are not always the same because Matrix_dimension is twice of Number_of_mesh_points for spin-triplet coupled channels.

File format of non-local 3N matrix elements in Jacobi momentum space
As mentioned in Sec.2.5, external input HDF5 files are needed for the non-local 3N matrix elements.The HDF5 files need to be prepared for each [J 3N  rel , (−1) l 12 +l 3 , T ] partial waves and placed in the directory (directory path)/T3_2T /J3_2J 3N rel /PAR_(−1) l 12 +l 3 , e.g., HOME/3NF_m atrix_elements_nonlocal_V/T3_1/J3_1/PAR_1 for the 3 H and 3 He ground-state channel.Each HDF5 file must include the entries Nalpha, Np, Nq, p mesh, q mesh, pw channels, and matrix elements, written as the dataset type supported in the HDF5 format.The Nalpha is the 32bit integer corresponding to the number of α channels.The objects Np and Nq are also 32-bit integers and the number of momentum mesh points for p and q, respectively.Note that α is defined in Eq. (46).The dataset p mesh should at least have mesh point and mesh weight entries, corresponding to the p-momentum mesh points and associated weights for a quadrature method, respectively.The objects mesh point and mesh weight should be Npdimensional array with the 64-bit float.The dataset q mesh is the same as p mesh except that it is for q.Regarding the dataset pw channels, it should at least include L_12, S_12, J_12, T_12, l_3, and 2*j_3 entries, which are l 12 , s 12 , j 12 , t 12 , l 3 , and 2 j 3 , respectively.Each entry has to be "Nalpha"-dimensional array with the 32-bit integer.Finally, the "matrix elements" corresponds to the 3N matrix stored as [Nalpha, Nq, Np, Nalpha, Nq, Np]-dimensional array with the 32-bit float 7 .

Summary and future perspective
We introduce the NuHamil code to generate NN and 3N matrix elements.The jobs can be managed with a simple Python script.The available NN interactions are LO -N 4 LO with 500 MeV regulator cutoff by Entem-Machleidt-Nosyk [26], N 3 LO with 500 MeV regulator cutoff by Entem-Machleidt [27], N 2 LO opt [28], N 2 LO sat [29], and ∆-full EFT series by Gothenburg-Oak Ridge collaboration [30].The code can generate locally regulated 3N interactions [36].Additional input files are needed for nonlocally regulated 3N interactions.The code also supports the free-space NN and 3N SRG evolution, and the consistent evolution of the other operators is implemented.The output files can be used for the NCSM calculations with the BIGSTICK code [20] and IMSRG calculations with the imsrg++ code [21].
For a comprehensive understanding of nuclear structure, interactions between a nucleus and external field should be addressed.For example the electromagnetic observables are results of the nucleus-photon interaction and are related with the multipole components of the electromagnetic current operators.Although we know that higher-order contributions are essential, see Ref. [49,50] for example, the LO current is used in most calculations due to the complexity of the matrix element calculations.As a future development, we plan to implement the higher-order current operators, including the two-body contributions.
consistent with the definition of the HO bracket.Note that π 1 is different from the usual relative momentum definition p = | p 1 − p 2 |/2.The NuHamil code requires the input file for V S J NN rel l ′ l