Impurity Lattice Monte Carlo for Hypernuclei

We consider the problem of including $\Lambda$ hyperons into the ab initio framework of nuclear lattice effective field theory. In order to avoid large sign oscillations in Monte Carlo simulations, we make use of the fact that the number of hyperons is typically small compared to the number of nucleons in the hypernuclei of interest. This allows us to use the impurity lattice Monte Carlo method, where the minority species of fermions in the full nuclear Hamiltonian is integrated out and treated as a worldline in Euclidean projection time. The majority fermions (nucleons) are treated as explicit degrees of freedom, with their mutual interactions described by auxiliary fields. This is the first application of the impurity lattice Monte Carlo method to systems where the majority particles are interacting. Here, we show how the impurity Monte Carlo method can be applied to compute the binding energy of the light hypernuclei. In this exploratory work we use spin-independent nucleon-nucleon and hyperon-nucleon interactions to test the computational power of the method. We find that the computational effort scales approximately linearly in the number of nucleons. The results are very promising for future studies of larger hypernuclear systems using chiral effective field theory and realistic hyperon-nucleon interactions, as well as applications to other quantum many-body systems.


Introduction
Hypernuclei are bound states of one or two hyperons together with a core composed of nucleons. They extend the nuclear chart into a third dimension, augmenting the usual two dimensions of proton number and neutron number. We will use the notation Y for a Λ or Σ hyperon and N for a nucleon. Due to the scarcity of direct hyperon-nucleon (Y N ) and hyperon-hyperon (Y Y ) scattering data, these unusual forms of baryonic matter play an important role in pinning down the fundamental baryon-baryon forces. This requires on the one hand an effective field theory (EFT) description of the underlying forces, as pioneered in Ref. [1,2], and on the other hand a numerically precise and consistent method to solve the nuclear A-body problem, such as nuclear lattice EFT (NLEFT) [3,4]. For calculations combining these chiral EFT forces at LO and NLO [5,6] with other many-body methods, see e.g. Ref. [7,8,9,10,11,12]. In view of the success of NLEFT in the description of nuclear spectra and reactions, it seems natural to extend this method to hypernuclei. However, this is not quite straightforward. While one can extend the four spin-isospin degrees of freedom comprising the nucleons to include the Λ and Σ states [13], this has not been done because there is no longer an approximate symmetry such as Wigner's SU(4) symmetry [23] that protects the Monte Carlo (MC) simulations against strong sign oscillations when using auxiliary fields. 1 The physics of hypernuclei therefore requires a different approach, and in this paper we show how the computational problems are solved using the impurity lattice Monte Carlo (ILMC) method.
The ILMC method was introduced in Ref. [15] in the context of a Hamiltonian theory of spin-up and spin-down fermions, and applied to the intrinsically non-perturbative physics of Fermi polarons in two dimensions in Ref. [16]. The ILMC method is particularly useful for the case where only one fermion (of either species) is immersed in a "sea" of the other species. Within the standard auxiliary field Monte Carlo method, such an extreme imbalance would lead to unacceptable sign oscillations in the Monte Carlo probability weight. In the ILMC method, the minority particle is "integrated out", resulting in a formalism where only the majority species fermions appear as explicit degrees of freedom, while the minority fermion is represented by a "worldline" in Euclidean projection time. The spatial position of this worldline is updated using Monte Carlo updates, while the interactions between the majority fermions are described by the auxiliary field formalism [4].
Here, we apply the ILMC method to the inclusion of hyperons into NLEFT simulations. We identify the Λ hyperon as the minority species, which we represent by a worldline in Euclidean time. This Λ worldline is treated as immersed in an environment consisting of some number of nucleons. We focus on the Monte Carlo calculation of the binding energy of light hypernuclei, by means of a simplified Y N interaction, consisting of a single contact interaction, tuned to a best description of the the empirical binding energies of the s-shell hypernuclei with A = 3, 4, 5. 2 For the N N interaction, we use a simple leading order interaction similar to that described in Ref. [17]. We benchmark our ILMC results against Lanczos calculations of transfer matrix and exact Euclidean projection calculations with initial states and number of time steps that match the ILMC calculations. We note that our Monte Carlo method is free from any approximation about the nodal structure of the many-body wave function. This is the first application of such unconstrained Monte Carlo simulations to hypernuclei. This paper is organized as follows. In Sec. 2, we present the Y N and N N interactions used for this study. In Sec. 3, we derive the impurity worldline formalism for the chosen Y N interaction, and introduce the concept of the "reduced transfer matrix", which refers to the nucleon degrees of freedom only. In Sec. 4, we discuss the Monte Carlo updating of the hyperon worldline and the auxiliary fields, which encode the interactions between nucleons. In Sec. 5, we present results for the ground state energies of the sshell nuclei and hypernuclei. In Sec. 6, we conclude with a discussion of future improvements and applications of the impurity lattice Monte Carlo method to hypernuclei and other quantum many-body systems.

Formalism
We develop the ILMC formalism following Ref. [15], who considered a system of spin-up and spin-down fermions, with a contact interaction which operates between fermions of opposite spin. The situation here is completely analogous, we have one majority species, the nucleons, and one impurity, the Λ. As usual in NLEFT, we consider positions on a spatial lattice denoted by n and lattice spacing a. We also assume that Euclidean time has been discretized, such that slices of the Euclidean time are denoted by n t with temporal lattice spacing a t . The partition function can be expressed in terms of the Grassmann path integral (1) where the subscripts N and Y refer to all nucleon and hyperon degrees of freedom, respectively. In this study we consider only Λ hyperons. In future work we will also consider Σ hyperons or account for their influence via three-baryon interactions involving a Λ and two nucleons. We also make the simplifying assumption that the hyperon-nucleon and nucleon-nucleon interaction are spinindependent and neglect Coulomb interactions.
Assuming that the exponent of the Euclidean action in Eq. (1) is treated by a Trotter decomposition, we find where the component due to the time derivative is while S Y and S N describe the kinetic energies of the hyperons and nucleons, respectively. Further, S Y N provides the Y N interaction, and S N N the N N interaction.
Derivations of Feynman rules are usually easier to perform in the Grassmann field formalism. However, actual NLEFT calculations are performed using the transfer matrix Monte Carlo method. As noted in Ref. [15], the Grassmann and transfer matrix formulations are related by where f is an arbitrary function, a † s and a s denote creation and annihilation operators for the fermion degrees of freedom, and colons signify normal ordering. We shall now consider the explicit forms of the Y N and N N interactions, and use Eq. (4) to relate expressions in the Grassmann and transfer matrix formulations.

The hyperon-nucleon interaction
For the hyperons, we take for simplicity the lowest-order (unimproved) kinetic energy where m Y is the hyperon mass, and we have defined α t ≡ a t /a as the ratio of temporal and spatial lattice spacings. The Y N interaction is given by in terms of the nucleon and hyperon densities, respectively. The tuning of the coupling constant C Y N is discussed in Section 5. Using Eq. (4), the hyperon contributions are described by the transfer matrix operator andρ are density operators for nucleons and hyperons, respectively. The a i,j (n) and a † i,j (n) are lattice annihilation and creation operators for nucleons on site n with spin i = 0, 1 (up, down) and isospin j = 0, 1 (proton, neutron).
The free Hamiltonian iŝ wherê denotes the (lowest order) kinetic energy for the hyperons in the operator formalism [3]. Here, theê l are unit vectors in lattice direction l. These lattice operators correspond to the continuum expressionŝ for the kinetic energy, and for the Y N interaction (with m Y and C Y N in physical rather than lattice units). Note that this is a simplified version of the pionless EFT calculation of Ref. [18], which also included a three-body interaction at LO. Such an interaction is sub-leading in chiral EFT approaches (such as NLEFT). See also the recent work in Ref. [19].

The nucleon-nucleon interaction
For the free Hamiltonian of the nucleon degrees of freedom, we likewise use the lowest-order expression and m N is the nucleon mass. The Wigner SU(4)-symmetric part of the leading-order (LO) N N interaction of Refs. [20,21,22] is used for the present work. This is an approximate symmetry [23] of the low-energy nucleon-nucleon interactions, where the nucleonic spin and isospin degrees of freedom can be rotated as four components of an SU(4) multiplet. Hence, we havê is the smeared nucleon density operator, and the (local) smearing function f sL is defined as and the operators a s NL † i,j (n) and a s NL i,j (n) are defined in terms of (non-local) smearing and where the values of the parameters C N N , s L and s NL for the present work are discussed in Section 5 (see also Ref. [17] for a full treatment). For ILMC calculations with the N N interaction included, we reduce the expressions quadratic in the density operators using the relation such that φ(n, n t ) is treated as a scalar auxiliary (Hubbard-Stratonovich) field. The N N term in the transfer matrix can then be written as for Euclidean time slice n t , where In the ILMC calculations, the path integral over the auxiliary field φ is evaluated using either local Metropolis algorithm updates or global lattice updates using the hybrid Monte Carlo algorithm. See Ref. [17] for details on efficient updating of the products of auxiliary-field transfer matrices.

Impurity worldline method
We shall now integrate out the hyperon degrees of freedom and derive a "reduced" transfer matrix, which refers to the nucleon degrees of freedom only. For simplicity (and without loss of generality), we shall neglect the N N interaction term for the purpose of the derivation, and consider the case of a single hyperon Y and nucleon N (which can be thought of as representing any one of the spin-isospin combinations i, j of the full theory).
Let us write down the transfer matrix element between time slices n t and n t + 1 in terms of where the χ s n t (n) count the occupation numbers for nucleons and hyperons on time slice n t and spatial lattice site n. Following the relations established in Ref. [15], we express the transfer matrix element as and are integers which assume values of either 0 or 1. Further, and are Grassmann functions (to be defined below). The impurity worldline is considered static for the purposes of this derivation, although it will be updated by the Metropolis algorithm in the actual Monte Carlo simulations. From one time slice to the next, the impurity may either remain on the same lattice site, or hop to a nearestneighbor site. For the case where the impurity remains on a given lattice site n , we have with and hence and we define / M (n t ) as the "reduced" transfer matrix. Specifically, we take with the nearest-neighbor expression Eq. (5) for the hyperon kinetic term. We also take as the Y N interaction. By evaluating the derivatives in Eq. (34), we find which we write as where the last factor, which encodes the interaction between the nucleons and the single hyperon impurity, has been exponentiated. Thus, Eq. (40) is the reduced Grassmann transfer matrix for the case where the impurity worldline remains stationary. For the case of a long-range Y N interaction, Eq. (40) should be replaced by an expression of the form whereby the hyperon impurity now also interacts with nucleons not on the same spatial lattice site. If we take then Eq. (40) for a contact interaction is recovered. Another possibility permitted by the nearest-neighbor Y N kinetic term is such that which gives for the reduced Grassmann transfer matrix, when the impurity hops to a neighboring lattice site.
Having determined the form of the reduced Grassmann transfer matrices, we may translate these to the transfer matrix formulation. The corresponding operators arê A few comments are in order about our implementation of the ILMC formalism. In our MC codes, Eq. (46) is evaluated aŝ where the prefactor (1 − 6h) and the Gaussian term from the Hubbard-Stratonovich transformation has not been written out. Note that Eq. (48) also includes the N N interaction through Eq. (25). For ILMC, the nucleons are treated as distinguishable particles, and the hyperon as a classical worldline during the Euclidean time evolution. This induces a three-body interaction when two nucleons and the hyperon occupy the same site, which is absent in Eq. (8). As we shall benchmark our ILMC codes against exact Euclidean time projection calculations of Eq. (8), we include the induced interaction to the original transfer matrix (8). This induced threebody interaction is a lattice artifact which disappears when α t → 0.

Monte Carlo calculation
We now describe how ILMC calculations are performed using the Projection Monte Carlo (PMC) method. Let us first assume that the impurity has been fixed at a given spatial lattice site, and that no "hopping" of the impurity occurs during the Euclidean time evolution. We shall then relax this constraint, and discuss a practical algorithm for updating the configuration of the hyperon worldline.

Stationary impurity
For a stationary hyperon impurity, the reduced transfer matrix is given by Eq. (46), and for the purposes of the PMC calculation, we define the Euclidean projection amplitude for a product of N t Euclidean time slices, where j and k denote different initial cluster states. As usual, this is expressed as a determinant of single-particle amplitudes, which gives where for p nucleons. By means of the projection amplitudes (51), we construct which is known as the "adiabatic transfer matrix". If we denote the eigenvalues of (53) by λ i (N t ), we find such that the low-energy spectrum is given by the "transient" energies at finite temporal lattice spacing a t . For the case of a single trial cluster state with p nucleons, Eq. (51) reduces to for the case of a single trial state. The ground-state energy is obtained from in the limit N t → ∞, where the exact low-energy spectrum of the transfer matrix will be recovered. Note that the argument N t + 1/2 is conventionally assigned to the transient energy computed from the ratio of projection amplitudes evaluated at Euclidean time steps N t + 1 and N t .
As an example, for the hypertriton we have p = 2 nucleons after the impurity hyperon has been integrated out. We start the Euclidean time projection with a single initial trial cluster state (j = k = 0) consisting of a spin-up proton, and a spin-up neutron. As there are no terms that mix spin or isospin, the other components of each single-particle state are set to zero, and remain so during the PMC calculation. For the spatial parts of the nucleon wave functions, we may choose, for example, the zero-momentum state in the notation of Ref. [15], which denotes plane-wave orbitals in a cubic box. In principle, we may also choose any other plane-wave state with non-zero momentum (see Table 1 of Ref. [15]), or any other more complicated trial state. For the heavier nuclei, it is indeed better to choose an initial state where the nucleons are clustered together.
In this case we sum over all possible translations of the cluster in order construct an initial state with zero total momentum.

Hopping impurity
If the hyperon impurity is allowed to hop between nearestneighbor sites (from one Euclidean time slice to the next), the Euclidean projection amplitude becomes a sum over hyperon worldline configurations. This gives where the product is expressed in terms of the reduced transfer matrices (46) and (47). Here, n j denotes the spatial position of the hyperon impurity (which has been integrated out) on time slice j. The expressions for the projection amplitude and determinant are generalized to where such that the determinant is now to be computed over all possible hyperon wordline configurations. We note that the worldline configuration is to be updated stochastically using a Metropolis algorithm. Thus, proposed changes in the impurity worldline are accepted or rejected by importance sampling with |Z jj (N t )| as the probability weight function. Here, j denotes one of the initial trial nucleon cluster states.

Worldline updates
The updating of the impurity worldline is handled in two steps: The generation of a new proposed worldline, and a Metropolis accept/reject step to determine whether to use the generated worldline. For this work, the worldline W (n, n t ) is a function of only the lattice site n and the Euclidean time step n t , and is equal to 1 where the impurity is present, and 0 at all other lattice points. From the expressions of the reduced transfer matrices, the worldline at two adjacent time steps, W (n , n t ) and W (n , n t + 1) must obey the relation |n − n | ≤ 1. For an illustration of the impurity (hyperon) worldline, see   Fig. 1. Illustration of the hyperon worldline. In the reduced transfer matrix formalism, the hyperon has been "integrated out", and the interaction between the hyperon and the nucleons is mediated by an effective "background field" generated by the hyperon worldline.
For the non-interacting worldline, we can generate new configurations from the free probabilities, as determined from the reduced transfer matrices. In this case, P h = h is the hopping probability, and P s = (1−6h) is the probability to remain stationary. When initializing the worldline at the beginning of the MC simulation, we may start from a configuration where the worldline is completely stationary ("cold start") or one where the worldline either hops or remains stationary at each time step according to the probabilities P h and P s ("warm start").
At the beginning of every sweep through the lattice, we propose a new worldline to use for that sweep. This is done by taking the previous worldline and choosing a random time at which we cut the worldline and regenerating it either in the forwards and backwards time direction. The new worldline is then accepted or rejected using a Metropolis accept or reject condition to preserve detailed balance associated with the absolute value of the amplitude.

Results
For the results presented in what follows, we use a spatial lattice spacing a = 1/(100 MeV) and temporal lattice spacing of a t = 1/(300 MeV). The non-local smearing parameter is chosen to be s NL = 0.2, and the local smearing parameter is set to s L = 0.0. Since we only consider s-shell nuclei and hypernuclei in this study, the local attraction provided by s L for heavier nuclei is not needed [20]. The coupling constant C N N is set to −7.5 × 10 −6 MeV −2 , and this combination of parameters yields a nucleon-nucleon scattering length a N N = 6.86 fm and effective range r N N = 1.77 fm. As stated before, in this study the spin-dependent terms of the nucleon-nucleon interaction are not accounted for.
For the Y N interaction, we set C Y N according to the best overall fit to the light hypernuclei. Fitting to the Λ separation energies for 3 Λ H, 4 Λ H/He, and 5 Λ He, we find C Y N = −1.6 × 10 −5 MeV −2 . This gives a Y N = −0.45 fm for the scattering length and r Y N = −0.45 fm for the effective range. In Table 1, we present benchmark calculations of the ILMC results for 3 Λ H in comparison with exact transfer matrix calculations. We show the results for the energy as a function of Euclidean projection time. We see that the agreement is quite good. The initial nucleon trial states for these calculations are taken to be spatially constant functions, which correspond to singleparticle states of zero momentum in a periodic cubic box.
The hyperon initial wave function is also taken be a constant function. These exact transfer matrix calculations include the induced three-baryon interaction described in Eq. (49).
In Table 2, we present exact Lanczos transfer matrix calculations of the ground state of 2 H, 3 Λ H, and separation energy B Λ , as a function of periodic box length. In this work, we also present the exact Lanczos transfer matrix calculation wherever it is computationally possible and using Monte Carlo for cases where it is not. Given the extremely small Λ separation energy, it is necessary to go to very large volumes in order to remove finite volume artifacts. Interestingly, B Λ is found to be relatively constant with the periodic box size L. This suppression of the finite volume dependence is an indication that the asymptotic normalization coefficient of the hypertriton wave function is small [24,25]. In Fig. 2, we present ILMC results for the 4 Λ H/He energy versus Euclidean time. These calculations use a periodic box size of L = 15.8 fm with up to N t = 300 Euclidean time steps. In order to extract the ground state energy, we use the extrapolation ansatz which takes into account the residual dependence of the first excited state that couples to our initial state. For this calculation, we use an initial state where the nucleon states have a spatially decaying exponential form with respect to the nucleus center of mass, while the initial hyperon wave function is a constant function. In Fig. 3, we show lattice Monte Carlo (LMC) results for the 4 He energy versus Euclidean time. As there are no hyperons in this system, these are auxiliary field Monte Carlo calculations without impurity worldlines. These calculations use a periodic box size of L = 9.9 fm with up to N t = 150 Euclidean time steps. In order to extract the ground state energy, we again use the exponential ansatz in Eq. (63). For this calculation, we again use an initial state where the nucleons have a spatially-decaying exponential form with respect to the nucleus center of mass.
In Fig. 4, ILMC results are shown for the 5 Λ He energy versus Euclidean time. These calculations use a periodic box size of L = 9.9 fm with up to N t = 250 Euclidean time steps. We again use the exponential ansatz from Eq. 63 to extract the ground state energy. Similar to the 4 Λ H/He calculation, here we use an initial state where the nucleons have a spatially decaying exponential form with respect to the nucleus center of mass, while the initial hyperon wave function is a constant function.
In Table 3, we present the lattice results for all of the s-shell nuclei and hypernuclei. The exact transfer matrix results are shown without error bars, while the ILMC and LMC results are shown with error bars that take into account stochastic errors and extrapolation errors. There is also a residual systematic error due to finite volume ef- fects. For a box size of L = 29.6 fm, the finite volume error on 2 H is 0.04 MeV, and the estimated finite volume error for 3 Λ H is also 0.04 MeV. As both corrections are in the same direction (with more binding at finite volume), the resulting finite volume error on the separation energy is < 0.002 MeV.
For a box size of L = 15.8 fm, the finite volume error on 3 H/He is 0.10 MeV, and the estimated finite volume errors for 4 Λ H/He are also 0.10 MeV. For a box size of L = 9.9 fm, the finite volume error on 4 He is 1.5 MeV, and the estimated finite volume errors for 4 Λ H/He are 2.0 MeV. Table 3. Summary of lattice results (exact transfer matrix, ILMC and LMC) for the energies of light nuclei and hypernuclei, and for separation energies. Comparisons with experimental separation energies are given where such data exists. These comparisons are averaged over Wigner SU(4) and Λ spin components. For the case of 4 Λ H/He, we average over the 0 + and 1 + separation energies for 4 Λ H and 4 Λ He weighted by number of spin components. More data can be found in the review Ref. [32].  [27,28,29] For the comparison with the experimental results, we average over Wigner SU(4) and Λ spin components where the data exists. We see that while the B exp Λ is larger than the experimental values for 3 Λ H and 5 Λ He, the separation is smaller than experimental value for 4 Λ H/He. This is an indication that there are deficiencies in our very simple treatment of the Y N and N N interactions. However, this serves as a good starting point for determining the essential features of the Y N interactions needed to describe the structure and properties of hypernuclei.

Discussion
We have shown, as a proof of principle, how state-of-theart NLEFT calculations can be extended to include hyperons. As the number of hyperons in realistic hypernuclei is small (typically one or two) relative to the number of nucleons, we have applied the ILMC method whereby the hyperon "impurity" is integrated out and represented by a hyperon "worldline", the position of which is updated during the MC calculation. Effectively, the standard NLEFT calculations for nucleons are augmented by a "background field" induced by the hyperon worldline. We have benchmarked the ILMC method by presenting preliminary MC results for the s-shell hypernuclei, using a simplified interaction similar to pionless EFT.
One of the most promising aspects of this work is the fact that the ILMC simulations scale very favorably with the number of nucleons. We have found that nearly all of the computational effort is consumed in calculating singlenucleon amplitudes as a function of the auxiliary field. As this part of the code scales linearly with the number of nucleons, it should be possible to perform calculations of hypernuclei with up to one hundred or more nucleons. We note also that the particular set of interactions that we have used here can also be directly applied to studying the properties of a bosonic impurity immersed in a superfluid Fermi gas. By modifying the included P -wave interactions of the impurity, we would also be able to describe the properties of an alpha particle immersed in a gas of superfluid neutrons. The possible applications of this method clearly go well beyond hypernuclear structure calculations and have general utility for numerous quantum many-body systems.
Returning to hypernuclear systems, the obvious next extension of this work is to include spin-dependent Y N interactions. The importance of the spin-dependence of the Y N interaction can be seen clearly in the splittings between the 0 + and 1 + states in 4 Λ H and 4 Λ He Ref. [26]. One should also include explicit ΛN -Σ 0 N transitions, see e.g. [33], as well as one-meson exchange interactions that would put the Y N interaction in the same EFT formalism [5,6] as currently used for the N N interaction in NLEFT [22].
The number of adjustable parameters in the Y N interaction will then increase. The most natural approach, in line with the treatment of the N N interaction, would be to fit such parameters to ΛN scattering phase shifts. However, due to the paucity of such data (especially at low energies), we expect to need at least the hypertriton binding energy as an additional constraint, as it is also done in continuum chiral EFT, see e.g. Ref. [6]. As the effects of ΛN -Σ 0 N transitions are included, it may be necessary to use further empirical data on other light hypernuclei to constrain the relevant LECs. A further extension concerns the extension to S = −2 hypernuclei, which on the one hand would involve the Y Y interactions [34,35,36] and on the other hand a modified ILMC algorithm for two interacting worldlines. Work along these lines is underway.