Hubbard-Stratonovich (HS) transform (or similar) for higher order-interactions

Question

I have a question about a generalization of the Hubbard-Stratonovich (HS) transformation to decouple two-body interactions. When dealing with Hamiltonians of the kind

\begin{equation} H = -\sum_{a}h_a \sigma_a - \sum_{a,b} \beta_{ab}\sigma_a \sigma_b = -\langle h,\sigma\rangle - \langle \sigma | \beta | \sigma \rangle \end{equation}

one can to decompose the matrix $\beta$ as $\beta = \sum_{k}\lambda_k |v_k\rangle \langle v_k|$. Then to compute the partition function one can do

\begin{equation} Z = \sum_{\{\sigma_a\}} e^{-H} = \sum_{\{\sigma_a\}} e^{\langle h,\sigma\rangle + \sum_{k}\lambda_k \langle \sigma |v_k\rangle \langle v_k | \sigma \rangle } = \sum_{\{\sigma_a\}} e^{\langle h,\sigma\rangle} \prod_{k}e^{\lambda_k \langle \sigma |v_k\rangle^2} \end{equation}

and then employ HS transform

\begin{equation} e^{a^2} = \int \,dx \, e^{-\frac{x^2}{2} + \sqrt{2} a x} \end{equation}

in the second term to decouple the square, bring in the summation and sum over the possible configurations.

My question is: what about more than binary interactions? What if the Hamiltonian is in the form (still linear)

\begin{equation} H = -\sum_{a}h_a \sigma_a - \sum_{a,b} \beta_{ab}\sigma_a \sigma_b + \sum_{a,b,c}\gamma_{abc} \sigma_a \sigma_b \sigma_c + \sum_{a,b,c,d} \delta_{abcd}\sigma_a \sigma_b \sigma_c \sigma_d + \ldots \end{equation}

Any help or insight would be very appreciated

Answer 1

You can generalise higher order interactions, but the calculations are harder. The HS transform is essentially a Fourier transform, which generalises well (for even $p$): $$ e^{-x^p/p} = \int e^{ikx}f_p(k)\frac{dk}{2\pi} \\ f_p(k) = \int e^{-ikx-x^p/p} dx $$ For example, if your higher order interactions can be written as: $$ H_{p,v} = \left(\sum_iv_i\sigma_i\right)^p $$ then you can use the previous Fourier transform to decouple the spins. There are two complications compared to HS.

Firstly $f_p$ is not as simple as in the quadratic case. It will typically involve hypergeometric functions. The only simple result is that at large $k$, you recover the Legendre transform of $x^n$ i.e. $$ f_p(k)\asymp e^{-k^q/q}\\ \frac1p+\frac1q=1 $$

A second hurdle is that you do not have the spectral theorem for higher order tensors. A reduction of the general interacting Hamiltonian in terms of sums of the form of $H_{p,v}$ is not guaranteed. One way around is to compute the multidimensional Fourier transform: $$ \exp\left(-\sum J_{i_1...i_p}\sigma_{i_1}...\sigma_{i_p}\right) = \int D\phi e^{i\phi\cdot \sigma}f_p(\phi) \\ f_p(\phi) = \int D\sigma \exp\left(-i\sigma\cdot \phi-\sum J_{i_1...i_p}\sigma_{i_1}...\sigma_{i_p}\right) \\ $$ but unlike the quadratic case, it is not obvious that this multidimensional integral can be factored into 1D integrals.