Sampling General N-Body Interactions with Auxiliary Fields

TL;DR

This paper introduces a versatile auxiliary field transformation that enables non-perturbative sampling of complex N-body contact interactions in lattice quantum Monte Carlo simulations, with controllable interaction strengths.

Contribution

The authors develop a novel auxiliary field transformation that generates all N-body contact terms and allows for their controlled, non-perturbative sampling in quantum Monte Carlo methods.

Findings

Reproduces exact solutions for a two-site quantum problem

Provides a controllable way to generate N-body interactions

Enables non-perturbative sampling of higher-order contact interactions

Abstract

We present a general auxiliary field transformation which generates effective interactions containing all possible N-body contact terms. The strength of the induced terms can analytically be described in terms of general coefficients associated with the transformation and thus are controllable. This transformation provides a novel way for sampling 3- and 4-body (and higher) contact interactions non-perturbatively in lattice quantum monte-carlo simulations. We show that our method reproduces the exact solution for a two-site quantum mechanical problem.

Tables

Table 1. Table 1: We show the results of the matching for N = 1 𝑁 1 N=1 , N = 2 𝑁 2 N=2 , and N = ∞ 𝑁 N=\infty with all c j > 3 = 0 subscript 𝑐 𝑗 3 0 c_{j>3}=0 . The first column is the Hubbard-Stratonovich case, which produces only a two-body interaction. We repeatedly used Γ ​ ( x ) = ( x − 1 ) ​ Γ ​ ( x − 1 ) Γ 𝑥 𝑥 1 Γ 𝑥 1 \Gamma(x)=(x-1)\Gamma(x-1) to simplify many N = 2 𝑁 2 N=2 coefficients, and use the shorthand γ 3 , 1 = Γ ​ ( 3 / 4 ) / Γ ​ ( 1 / 4 ) subscript 𝛾 3 1 Γ 3 4 Γ 1 4 \gamma_{3,1}=\Gamma(3/4)/\Gamma(1/4) . We provide a Mathematica notebook useful for generating the λ ( M ) superscript 𝜆 𝑀 \lambda^{(M)} for a given N 𝑁 N in the Supplementary Material.

	$N=1$	$N=2$	$N=\infty$, $c_{j>3}=0$
${\lambda }^{\left(1\right)}$	$0$	${\gamma }_{3,1}{c}_2$	$\frac{1}{3}{c}_2$
${\lambda }^{\left(2\right)}$	$-\frac{1}{4}{c}_1^2$	$-\left(\frac{1}{8}-\frac{{\gamma }_{3,1}^2}{2}\right)c_2^2-\frac{1}{4}{c}_1{c}_3-\frac{1}{2}{\gamma }_{3,1}{c}_1^2-\frac{3}{8}{\gamma }_{3,1}{c}_3^2$	$-\frac{1}{6}{c}_1^2-\frac{2}{45}{c}_2^2-\frac{1}{5}{c}_1{c}_3-\frac{1}{14}{c}_3^2$
${\lambda }^{\left(3\right)}$	$0$	$\left(\frac{1}{8}-\frac{{\gamma }_{3,1}^2}{2}\right)c_1^2{c}_2+\left(\frac{5}{32}-\frac{3{\gamma }_{3,1}^2}{8}\right)c_2{c}_3^2$	$\frac{2}{45}{c}_1^2{c}_2+\frac{8}{2835}{c}_2^3+\frac{8}{105}{c}_1{c}_2{c}_3+\frac{2}{63}{c}_2{c}_3^2$
		$+\frac{1}{2}{\gamma }_{3,1}{c}_1{c}_2{c}_3+\frac{1}{3}{\gamma }_{3,1}^3{c}_2^3$
${\lambda }^{\left(4\right)}$	$0$	$\left(\frac{{\gamma }_{3,1}^2}{8}-\frac{1}{96}\right)c_1^4-\frac{{\gamma }_{3,1}^3}{2}{c}_1^2{c}_2^2-\left(\frac{3}{64}-\frac{3{\gamma }_{3,1}^2}{16}\right)c_1^2{c}_3^2$	$\frac{1}{180}{c}_1^4+\frac{4}{14175}{c}_2^4+\frac{1}{105}{c}_1^3{c}_3-\frac{52}{10395}{c}_2^2{c}_3^2-\frac{5}{7644}{c}_3^4$
		$-\frac{{\gamma }_{3,1}}{8}{c}_1{c}_3^3-\left(\frac{1}{8}-\frac{{\gamma }_{3,1}^2}{2}\right)c_1{c}_2^2{c}_3-\frac{{\gamma }_{3,1}+3{\gamma }_{3,1}^3}{8}{c}_2^2{c}_3^2$	$-\frac{4}{945}{c}_1^2{c}_2^2+\frac{13}{3150}{c}_1^2{c}_3^2-\frac{16}{1575}{c}_1{c}_2^2{c}_3-\frac{1}{1155}{c}_1{c}_3^3$
		$-\left(\frac{1}{192}-\frac{{\gamma }_{3,1}^4}{4}\right)c_2^4-\left(\frac{15}{512}-\frac{9{\gamma }_{3,1}^2}{128}\right)c_3^4$

Equations

$${\mathcal{Z}}_{c,N}\equiv \underset{-\infty }{\overset{\infty }{\int }}\text{d}\varphi P_N(\varphi )\exp \left\{-\sum _{j=1}^{2N-1}{c}_j{\varphi }^j\hat{\rho }\right\},$$

$$P_N(\varphi )=\frac{N}{\Gamma \left(\frac{1}{2N}\right)}{e}^{-{\varphi }^{2N}}.$$

$$\underset{-\infty }{\overset{\infty }{\int }}\text{d}\varphi e^{-{\varphi }^{2N}}{\varphi }^{2k}=\frac{\Gamma \left(\frac{1+2k}{2N}\right)}{N}\stackrel{N\to \infty }{⟶}\frac{2}{1+2k},\forall k\in {\mathbb{N}}_0,N\in \mathbb{N}.$$

$${\mathcal{Z}}_{\lambda }$$

$${\mathcal{Z}}_{c,N}$$

$${\mathcal{Z}}_{c,N}^{(M)}$$

$${\mathcal{M}}_{NM}^{(2k)}=$$

$$\wedge \left(\sum _{j=1}^{2k}{m}_j=M\right)\wedge \left(j\ge 2N\Rightarrow m_j=0\right)\},$$

$${\mathcal{Z}}_{\lambda }=\sum _{M=0}^{\infty }{\mathcal{Z}}_{\lambda }^{(M)}{\hat{\rho }}^M,{\mathcal{Z}}_{\lambda }^{(M)}=\sum _{\vec{m}\in {\mathcal{M}}^{(M)}}\prod _{k=1}^M\left[\frac{{\left(-{\lambda }^{(k)}\right)}^{m_k}}{m_k!}\right]$$

$${\mathcal{M}}^{(M)}=\left\{\vec{m}\in \text{mathbbm}{N}_0^M|\sum _{k=1}^{M}k{m}_k=M\right\}.$$

$${\mathcal{Z}}_{\lambda }^{(M)}$$

$$\hat{H}=\kappa \left(-3\sum _{f=1}^3\sum _{⟨i,j⟩}{\hat{a}}_{f,i}{\hat{a}}_{f,j}^{†}+{\lambda }^{(2)}\sum _{i\in \{0,1\}}{\hat{\rho }}_i^2+{\lambda }^{(3)}\sum _{i\in \{0,1\}}{\hat{\rho }}_i^3\right),$$

$$E=-\underset{\tau \to \infty }{\lim }{\partial }_{\tau }\log Z[\tau ,{\Psi }_T],$$

$$Z[\tau ,{\Psi }_T]$$

$$=\int \left(\prod _{x}d{\varphi }_x{P}_N({\varphi }_x)\right)\mathcal{K}[\tau ,\vec{\varphi },c_j,{\Psi }_T],$$

$$N$$

$$\mathcal{K}[\tau ,\vec{\varphi },c_j,{\Psi }_T]={\left(K^{-1}[\tau ;\vec{\varphi },c_j]\right)}^{\alpha }$$

$$K^{-1}[\tau ;\vec{\varphi },c_j]$$

$$-K^{-1}[1,\tau ;0,0;\vec{\varphi },c_j]+K^{-1}[1,\tau ;1,0;\vec{\varphi },c_j]),$$

$$K[i,t;i^{\prime },t^{\prime };\vec{\varphi },c_j]\equiv {\delta }_{i,i^{\prime }}{\delta }_{t,t^{\prime }}+a_{\tau }{\delta }_{⟨i,i^{\prime }⟩}{\delta }_{t,t^{\prime }+1}$$

$$+{\delta }_{i,i^{\prime }}{\delta }_{t,t^{\prime }+1}\left(c_1{\varphi }_{i,t}+c_2{\varphi }_{i,t}^2+c_3{\varphi }_{i,t}^3-1\right).$$

$$E(\tau )\equiv -\frac{1}{a_{\tau }}\log \left(\frac{Z[\tau +a_{\tau },{\Psi }_T]}{Z[\tau ,{\Psi }_T]}\right)$$

Full text

The introduction of auxiliary fields, e.g. via the Hubbard-Stratonovitch (HS) transformation Stratonovich (1957); Hubbard (1959), as a means of linearizing the interaction term of the Hamiltonian in terms of its density operator is common practice in many areas of theoretical physics, and is particularly prevalent in condensed matter and nuclear physics Gubernatis et al. (2016); Kleinert (2011). Upon linearization of the theory, the problem becomes that of many particles undergoing one-body interactions with the fluctuating background auxiliary field. In occupation number formalism, the problem becomes very similar to that of non-interacting (amongst themselves) particles, which in some cases can be evaluated via steepest-descent methods (e.g. 1-D Ising model). At the very least, the transformation greatly facilitates numerical treatment of the many-body problem Lee (2009); Epelbaum et al. (2010, 2009); Carlson et al. (2015). In path integral formalisms that involve Grassmann fields, the transformation is essential as it reduces the Lagrangian to terms bilinear in fermionic fields, which can subsequently be integrated out exactly via Grassmann gaussian integration (e.g. bosonization of fermionic theories Kopietz (2008)). The HS transformation is ideal for theories that initially have 2-body interactions (i.e. terms quadratic in the density operator). In principle many-body interactions can be linearized through recursive application of the HS transformation, but at the cost of introducing numerous auxiliary fields. A transformation that naturally includes 3-body interactions, for example, would be beneficial for studying ultra-cold gases of polar molecules, where 3-body (and higher) interactions can be tuned to become dominant Büchler et al. (2007); Bonnes et al. (2010).

In this Letter we detail a generalization of the HS transformation that includes two-body and all possible higher contact interactions that involves only one auxiliary field. The standard HS transformation can be viewed as a particular limit of this general transformation. Our transformation, in principle, induces all possible $n$-body interactions $\lambda^{(n)}\hat{\rho}^{n}$ where $\hat{\rho}$ is the density operator. The coefficients $\lambda^{(n)}$ are analytically determined and controllable through a set of accompanying coefficients $c_{j}$ that control the coupling of the density operator with the $j^{\text{th}}$ power of the auxiliary field in the linearized theory. The numerical implementation of this transformation is trivial as the sampling of the fields can be obtained from known distributions. For 2-flavor fermionic systems (e.g. nucleons), this transformation allows for complete control of contact interactions up to 4-body in nature.

Formalism. Consider the following integral which is intended to correspond to the interacting part of a partition function at a single space-time point

[TABLE]

where $P_{N}(\phi)$ is the normalized probability distribution

[TABLE]

Here, $\phi$ is an auxiliary field which couples to the fermionic density operator $\hat{\rho}=\sum_{f}{\bar{\psi}}_{f}\psi_{f}$, where $f$ runs over the different fermion species at a given site. Because the density operators at different spacetime points commute, we present only the derivation for a single point.

The argument of the exponential in (1) describes interaction vertices with an incoming and outgoing fermionic field (the density operator) and from one up to $2N-1$ auxiliary fields associated with interaction strength $c_{j}$. The integral is normalized such that $\mathcal{Z}_{c,N}\to 1$ for $c_{j}\to 0$. The result of the integration over the largest exponent of $\phi$ times another polynomial in $\phi$ is given by

[TABLE]

It is sufficient to only consider polynomials in even powers of $\phi$ because of the symmetry of the integration – odd powers vanish. Also, in order for the integral to converge, the leading exponent is not allowed to be odd in $\phi$.

In this work we identify the integral $\mathcal{Z}_{c,N}$ with an effective action consisting of general $2k$-fermionic field vertices

[TABLE]

We expand $\mathcal{Z}_{c,N}$ and $\mathcal{Z}_{\lambda}$ in powers of the density operator $\hat{\rho}$ and systematically match order-by-order, relating the auxiliary field interactions to the induced many-body forces. Using Faà di Bruno’s formula Faà di Bruno (1855) we find that

[TABLE]

where the sum runs over the set

[TABLE]

and

[TABLE]

where the sum in $\mathcal{Z}_{\lambda}^{(M)}$ is over the set

[TABLE]

As the highest coefficient $\lambda^{(k)}$ (in terms of $k$) in (8) is given by $\lambda^{(M)}$ and is also linear in exactly this coefficient (because if $m_{M}=1$ then $m_{j\neq M}=0$), one can recursively determine all coefficients $\lambda^{(M)}$ for $M>0$ by

[TABLE]

Note that one can prove by induction that each coefficient $\lambda^{(M)}$ is proportional to a sum where each term is a products of $M$ fermion auxiliary-field coefficients: $\lambda^{(M)}\propto c_{j_{1}}\cdots c_{j_{M}}$ Körber et al. (2017). In Tab. 1 we show the coefficients of the induced forces up to the order of four-body forces ($M=4$) for three different choices of $N$.

Numerical Results.

To demonstrate the efficacy of our transformation, we consider a system of 3 different fermion species $f$ interacting on a two-site model (with sites $\{0,1\}$) with the following Hamiltonian,

[TABLE]

where $\hat{a}^{\dagger}_{f,i}$ ($\hat{a}_{f,i}$) is a creation (annihilation) operator for a fermion of species $f$ at site $i$ and $\hat{\rho}_{i}=\sum_{f}\hat{a}^{\dagger}_{f,i}\hat{a}_{f,i}$ is the number operator at site $i$. Here the sum over $\langle i,j\rangle\in\{(0,1),(1,0)\}$ represents nearest neighbor hoppings between the two sites, and $\kappa$ is a dimensionful parameter that sets the dynamical scale of our problem. As is evident in (11), $\lambda^{(2)}$ and $\lambda^{(3)}$ represent the size of two- and three-body interactions (relative to $\kappa$), respectively, and are dimensionless. A positive (negative) value of $\lambda^{(i)}$ indicates a repulsive (attractive) interaction. We quote results in units of $\kappa$. The two-site model can be diagonalized and its spectrum directly determined from its eigenvalues. We use the transfer matrix formalism Lee (2009) to obtain this spectrum 111We compare the spectrum obtained from the transfer matrix (instead of direct diagonalization of the Hamiltonian) since it incorporates lattice discretization effects. This enables a direct comparison with our lattice projection calculations.. With the spectrum in hand, we compare these results with those obtained from a stochastic lattice projection calculation Lee (2009) where the introduction of auxiliary fields via our transformation is needed. A detailed description of these calculations is in preparation Körber et al. (2017). In what follows, we give a succinct description.

The projection method extracts the lowest energy level $E$ in the spectrum of the system via

[TABLE]

where $Z$ is given by

[TABLE]

$\Psi_{T}$ is an initial trial wavefunction, and $\vec{\phi}$ indicates the collection of all $\phi$ over $x$, all space (both sites) and time from 0 to $\tau$.

In the last step of the (13) we have applied our transformation with its associated probability distribution $P_{N}(\phi)$ at each spacetime point and introduced $\mathcal{K}$ which is a functional of $\phi$ and depends explicitly on the coefficients $c_{j}$, the trial wavefunction $\Psi_{T}$, and the time separation $\tau$. The form of $P_{N}(\phi)$ depends on the order $N$ of the transformation that is applied. For the work presented here, we consider the two extreme cases

[TABLE]

The case $N=1$ corresponds to the gaussian distribution of the original Hubbard-Stratonovich transformation (and thus only $\lambda^{(2)}$ is nonzero), whereas for $N=\infty$ we have uniform sampling and in principle all allowed many-body contact interactions. In the $N=\infty$ case, there are in principle an infinite number of $c_{j}$—we set them all to zero for $j>3$.

As our trial wavefunction for multifermion systems $|\Psi_{T}\rangle$ we pick a direct product state of single-particle ground state wavefunctions, $|\psi_{f}\rangle=\frac{1}{\sqrt{2}}[a^{\dagger}_{f,0}-a^{\dagger}_{f,1}]|0\rangle$ for each flavor $f$. The direct product structure allows us to write the functional $\mathcal{K}$ in (13) as

[TABLE]

where $\alpha$ is the number of fermion types and

[TABLE]

After introducing a time discretization $\tau=a_{\tau}t$, the $K$ matrix is given by

[TABLE]

Here we have used the fact that we let $c_{j}=0$ for $j>3$. We attach the indices $i$ and $t$ on the auxiliary field to indicate that there exists auxiliary fields for each spatial and temporal point.

Because we work with discretized time, we analyze the discretized version of eq. (12),

[TABLE]

and search for constant plateaus at long times to numerically determine the $\tau\rightarrow\infty$ limit.

In Fig. 1 we show the numerical results of $E(\tau)/\kappa$ for combinations of $\lambda^{(M)}$ where $\lambda^{(2)}\in\{\pm 1.74,0.0\}$ and $\lambda^{(3)}\in\{\pm 0.245,0.0\}$ for two and three-fermion systems. For calculations with $\lambda^{(3)}=0.0$, $\phi$ was sampled according to the HS distribution $P_{1}$ and all other calculations sampled according to $P_{\infty}$. Note that the coefficients $c_{j}$ can be complex. The solid lines correspond to the exact answers from diagonalization of the transfer matrix, and the dashed lines correspond to non-interacting energies. We note that for all two-fermion calculations, both $N=1$ and $N=\infty$ calculations agree with each other, regardless of the value of $\lambda^{(3)}$, since the 3-body interaction plays no role. Furthermore, with $\lambda^{(2)}=0.0$ but nonzero $\lambda^{(3)}$, the two-fermion system reproduces the non-interacting result. In the three-fermion system, the effects of $\lambda^{(3)}$ are apparent and agree well with exact diagonalization. A quantitative analysis will be given in Körber et al. (2017), as well as an analysis of signal-to-noise behavior for the different systems.

Discussion. The Hubbard-Stratonovitch transformation has been a useful tool for making analytic and numerical progress in problems of physical interest, and is a special case of the transformation described in this work. Our transformation, which uses a self-interacting auxiliary field, allows for the direct inclusion of controllable many-body forces into numerical calculations.

Choosing a degree of auxiliary field self-interaction fixes the sampling distribution and limits the types of fermion/auxiliary-field interaction vertices. These interactions in turn generate a slew of $n$-body forces. Sampling a gaussian distribution for the auxiliary field recovers the original HS transformation, while uniform sampling in principle allows for the independent control of all possible $n$-body forces, while intermediate distributions yield correlated forces.

We have demonstrated that on a two-site model, different sampling techniques reproduce the exact results in a variety of cases, correctly handling all combinations of attractive, repulsive, and absent two-body forces with attractive and repulsive three-body forces. While even the agreement of the one-body energies found with the different methods is nontrivial, the agreement of the two-body systems between methods and agreement of three-body systems with exact results indicate a reliable understanding and definite control over many-fermion interactions.

One drawback of our method is that some combinations of strengths of many-body forces are only achievable with complex $c_{j}$, opening the possibility of introducing a numerical sign problem Chen et al. (2004). Because the $c$ coefficients appear nonlinearly in the higher-body $\lambda$ coefficients, there is not a unique set of $c$s that yield a particular set of $\lambda$s.

It would be interesting to consider sampling by more generic functions such as $P_{N}(\phi)=\exp(-\sum_{n}\alpha^{(n)}\phi^{n})$ with multiple tuneable $\phi$ self-interactions. For field-theoretic applications, it may prove fruitful to understand how the renormalization of different $\lambda^{(M)}$ control renormalization of the parameters $c_{j}$.

We expect our transformation to be useful in a variety of physical systems where many-body interactions are important. Numerical applications may include density functional theory approaches to nuclear physics Drut et al. (2010), the nonperturbative inclusion of multinucleon forces into NLEFT which might unlock precision characterizations of halo nuclei, the study of systems near the Efimov threshold such as cold atoms (see Refs. Braaten and Hammer (2007); Naidon and Endo (2017) and references therein), systems at high density, and any other system where contact interactions need stochastic implementation.

Appendix A Acknowledgements

We thank A. Rokash, D. Lee and A. Cherman for insightful discussions. This work was done in part through financial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110).