Hubbard-Stratonovich-like Transformations for Few-Body Interactions

Through the development of many-body methodology and algorithms, it has become possible to describe quantum systems composed of a large number of particles with great accuracy. Essential to all these methods is the application of auxiliary fields via the Hubbard-Stratonovich transformation. This transformation effectively reduces two-body interactions to interactions of one particle with the auxiliary field, thereby improving the computational scaling of the respective algorithms. The relevance of collective phenomena and interactions grows with the number of particles. For many theories, e.g. Chiral Perturbation Theory, the inclusion of three-body forces has become essential in order to further increase the accuracy on the many-body level. In this proceeding, the analytical framework for establishing a Hubbard-Stratonovich-like transformation, which allows for the systematic and controlled inclusion of contact three- and more-body interactions, is presented.


Introduction
The Hubbard-Stratonovitch (HS) transformation [1,2] is a common tool in many areas of theoretical physics, ranging from condensed matter to nuclear physics [3,4]. Though the basic idea-completing the square-is quite simple, the practical benefits can be numerous. It is used to linearize the Hamiltonian in terms of its density or number operator by the introduction of an auxiliary field, which eventually is integrated out. Thus, a many-body system can effectively be described as a one-body system which interacts with a background field, drastically simplifying the numerical description of many-body problems [5][6][7][8].
In this proceeding, which closely follows [9], we generalize the HS transformation to linearize arbitrary many-body systems by the introduction of just one auxiliary field, which is distributed according to specific probability distributions that allow the inclusion of different auxiliary field interactions. The presentation is organized as follows: in section 2 we review the basic idea of the original HS transformation from a more diagrammatic perspective and explore how it could possibly be generalized. In section 3 we discuss the formalism which enables the identification of a linearized Hamiltonian including auxiliary fields which eventually corresponds to a specific many-body Hamiltonian after integrating out the auxiliary fields. A case study for an analytically solvable two-site model is presented in section 4 and we finally discuss our results in section 5.

The Basic Idea
In the original Hubbard-Stratonovich transformation one linearizes the action by the introduction of a background field φ ∈ R, which couples to the usual fermionic (or bosonic) fields ψ. From now on, we denote the fermionic (or bosonic) bilinear by the density operatorρ ≡ψψ. After integrating over this auxiliary field, one recovers the original interaction From a diagrammatical point of view, the auxiliary field can be interpreted as a bosonic spin-zero field with a Yukawa-like fermion interaction (and a constant propagator). If one interprets the right-hand-site of equation (1) as a path integral, the diagrammatic interpretation would be the sum of products of two-fermion diagrams which exchange an auxiliary field. Integrating out the auxiliary field, which diagrammatically corresponds to contracting the auxiliary field lines, effectively generates the desired two-body fermion interaction (see figure 1). The generalization of this diagrammatic interpretation is the basic idea to include three-body forces. We start with a possible diagram which consists of one-body fermion and any possible Nbody auxiliary field interaction.  Figure 2: Allowed diagrams and induced forces due to transformation with a quadratic auxiliary field, one-body fermion interaction.
The simplest example is the case of a one-body fermion, two-body auxiliary field interaction (figure 2a) which can generate a three-body fermion interaction (figure 2b). In order to guarantee a convergent integrand, the highest power of auxiliary fields appearing in the exponential should be at least φ 4 . After contracting the auxiliary field lines, one obtains the desired three-fermion force (figure 2c). However, when one actually computes all possible diagrams, one generates an infinite tower of different N-fermion diagrams (see figures 2d, 2e ...). The only way to include three-and higher-body forces in a controlled fashion is to identify the induced N-fermion couplings, depending on the auxiliary field interaction coefficients, and, at best, find a convergent and controllable expansion (2)

The Formalism
In order to stay as general as possible, we do not limit ourselves to quadratic interactions between the auxiliary field and the fermion bilinears. The most general form with one type of auxiliary-field self-interaction is given by the following integral whereρ = fψ f ψ f is the fermion density operator, f runs over the different fermion species at a given site and P N (φ) is the normalized probability distribution. This integral is well defined for any coefficients c i ∈ C and N ∈ N. It is sufficient to consider auxiliary fields at just one point in spacetime, because the the density operators at different sites commute and the auxiliary fields get no kinetic term.
The argument of the exponential in (3) describes all possibly allowed interactions between the one-fermion operatorρ and auxiliary fields, with j = 1 . . . 2N − 1 open auxiliary field lines of strength c j . Accordingly, the integral is normalized in the case of c j = 0 ∀ j ⇒ Z c,N = 1, which corresponds to the absence of auxiliary fields. For any c j 0, one effectively induces N-body forces after integrating out the auxiliary fields. Hence, to find a controlled description of N-body fermion forces λ (k) in terms of auxiliary field coefficients c j , one has to match both expressions order by order in the operator expansion of the exponential The expansion of equation (3) involves integration over the probability distribution times another polynomial with arbitrary powers of φ Note that odd powers in φ vanish by symmetry. Using Faà di Bruno's formula [10] ∂ we simultaneously expand equation (3) in φ andρ and equation (2) inρ only Here we have defined the set where m denotes the sum over all components of m. The set M (M) is the same as defined in equation (6). Because one can easily verify 1 that the Z (M) λ is linear in λ (M) , one can recursively determine the dependence of λ (M) on the previous coefficients λ (k) with k < M and the expansion in the c coefficients by requiring We want to highlight that all of the obtained results are smooth in the limit of N → ∞, which would allow one to incorporate an infinite amount of freely tuneable c j -coefficients. In this limit, the probability distribution for the auxiliary fields becomes uniform: P ∞ (φ) = 1/2 for φ ∈ (−1, 1) and zero otherwise. However, for every different order N, the dependence of the induced forces λ (k) on the auxiliary coefficients c j slightly changes. We present an example with N = 2 and N = ∞ in figure 3.
Because the analytic form of the induced forces that depend on the auxiliary field interaction parameters can be tuned in a completely general fashion, we provide a Mathematica notebook in the Supplementary Material of Ref. [9] which can be used to compute the λ (M) for the problem of choice.
Furthermore, we emphasize that the infinite tower of induced interactions can be systematically controlled. It can be inductively shown 2 that the coefficient λ (k) ∝ c i 1 c i 2 · · · c i k . Thus, if the magnitude of c j coefficients is smaller than one, higher forces become increasingly smaller.

Numerical Results
In order to demonstrate the applicability of the proposed method for including N-body forces in lattice algorithms, we set up a toy problem which can be solved analytically and compare numerically sampled results to their analytical counterparts. As a test system, we choose a two-site model (sites i = 0 and 1) and place up to 3 different fermion species f at each site. The Hamiltonian which describes this system is given bŷ Here,â f,i (â   different values in order to explore the capabilities of the proposed transformation. We analytically compute the largest eigenvalue of the transfer matrixT (τ) ≡ : exp(−τĤ) : (see [5]), and compare the result to the large time projection of the transfer matrix acting on an initial target state Ψ T We compute the large time projected matrix element using the proposed prescription for auxiliary fields. In other words we express the Hamiltonian of equation (10) in terms of effective couplings c j and fermion-bilinear auxiliary field interactions 3 . The coefficients c j are matched to the effective twoand three-body forces by the relations listed in figure 3 and the remaining integral is computed by sampling the auxiliary field according to distribution P N (φ) (same distribution for both lattice sites). Finally the function K is obtained by sandwiching the auxiliary-field-transformed transfer matrix with an (arbitrary) initial wave function, which we have chosen as   : E(τ)/κ for different systems. All systems were studied with 2 × 10 8 measurements. Some systems show a signal-to-noise problem and others yield an easy extraction of a constant plateau in the long-time-limit. The left panels all show two-fermion systems and the corresponding right panels show the three-fermion systems with the same parameters. In the top (middle) [bottom] two panels we show a system with repulsive (absent) [attractive] two-body forces. Red (blue) [green] points correspond to repulsive (absent) [attractive] three body forces. The blue points were sampled according to the HS distribution P 1 and the other points according to P ∞ . The data in the middle two panels were generated with P ∞ and coefficients c j = 0 for j > 3 tuned to exactly cancel the two-body force, and we show dashed lines for the corresponding non-interacting energies. This figure is taken from [9]. and generalized transformation following the uniform distribution (N = ∞), for attractive, absent and repulsive two-body forces λ All of the results, in the large-time limit, agree well with the analytical results within uncertainties. However, computations for specific values of the induced two-and three-body forces suffer from relatively large statistical uncertainties. This is expected and can be explained by the Hamburger truncated momentum problem [11,12]. In forthcoming work [13] we intend to analyse such fluctuations in terms of signal-to-noise behaviour and their relation to the sign problem in more detail.

Discussion
We have presented a general prescription which generalizes the Hubbard-Stratonovich transformation such that any N-body contact interactions for fermions (and similarly for bosons) can be cast to a set of linearized fermion (boson) interactions with an auxiliary field. Depending on the distribution of auxiliary fields, ranging from gaussian to uniform, more possible interactions with the auxiliary field can be included. Each of the linear interactions with the auxiliary field comes with an a priori freely tunable interaction coefficient that can be analytically determined to reproduce the original forces in a controled fashion after integrating out the auxiliary fields. If one chooses a proper sampling distribution and as a consequence obtains a sufficient amount of auxiliary interactions, the set of different interaction coefficients which reproduce the original forces is infinite.
We verify these predictions for a simple two-site model and confirm that statistical simulations, using the presented prescription, converge against their analytical result within uncertainties. This non-trivial cross-check is demonstrated for a set of attractive, absent and repulsive two-and threebody interactions. We find that certain combinations of interaction coefficients are more favourable than others in terms of convergence of statistical observables: e.g. one should prefer real over complex coefficients. Furthermore, for certain effective forces, it is not possible to find sets of interaction coefficients which do not suffer from large statistical fluctuations.
For upcoming work, we intend to further investigate the signal-to-noise behaviour for different sets of interaction coefficients. Also, it might be interesting to generalize this transformation even further to include non-central and non-local interactions or understand the renormalisation induced by the effective forces on the level of interaction coefficients.
Last but not least, we expect that our findings impact a variety of physical systems for which the inclusion of many-body forces is relevant. Examples include the non-perturbative addition of threebody forces in nuclear lattice effective field theory, the study of systems near the Efimov threshold and studies of halo nuclei.