is there a way to extend the HubbardStratonovich transformation $$e^{\frac{1}{2}Ks^2}=\left(\frac{K}{2\pi}\right)^{1/2}\int_{–\infty}^\infty e^{–Kx^2+Ksx}dx$$ to the case $e^{\frac{1}{2}Ks^p}$ for $p \in \mathbb{N}$?
2 Answers
@HyyFly: It would be nice to elaborate a bit on your question and in particular give some motivation. What applications do you have in mind for possible generalizations of the HubbardStratonovich transformation? The latter is just the physics parlance for the formula for the Fourier transform of a Gaussian, used in reverse. In any case, if you want to entice people to write useful and thoughtful answers, how you write your question must show some effort.
As to the question itself:
The problem of figuring out when "simple functions" like the exponential of a quadratic has a simple Fourier transform was addressed in the article "When is the Fourier transform of an elementary function elementary?" by Etingof, Kazhdan and Polishchuk in Selecta Math. 2002.
If all you want is a way of breaking quantum field theory vertices with $p$ legs into trivalent vertices, then one can do that by iterating the HubbardStratonovich transformation (aka intermediate field representation). See, for instance, the article "Loop vertex expansion for $\Phi^{2k}$ theory in zero dimension" by Rivasseau and Wang in J. Math. Phys. 2010 and the follow up "Note on the intermediate field representation of $\phi^{2k}$ theory in zero dimension" by Lionni and Rivasseau.

$\begingroup$ For a QFT where you want to split fermion multilinears into a bilinear and auxiliary field there is also arxiv.org/abs/1706.06494 $\endgroup$– evanbCommented Aug 10, 2020 at 21:01
As Abdelmalek Abdesselam mentions you need to iterate such transformations. Let me sketch one way this can work, (in a formal perturbative sense):
Starting from $\exp(s^p)$, you can consider the integral over two new variables $x,y$ $$ \int dx \int dy \exp(s^{p2} y + x(ys^2)) $$ Clearly doing the $x,y$ integrations together will set $y=s^2$ (and there will also be some numerical prefactor you can adjust) so you obtain the original term $\exp(s^p)$.
However the log of the integrand here is now of order $p1$ rather than $p$, so you can do another such transformation involving $s,x,y$ and more new variables to further reduce the total order until you hit $p=3$, where the procedure terminates.
There is a mathematical physics interpretation here, if you are interested: if we view $\exp(s^p)$ as the (exponential of an) action of some physical system, and canonically associate to that an Linfinity algebra, what this procedure is doing is producing a strictification of said algebra: namely an Linfinity algebra quasiisomorphic (read: physically equivalent) to the original one, but all of whose higher (than bilinear) brackets vanish. The end result looks like a ChernSimons type theory.