Solution (source code)

= Solution

Use the <Fourier transform> convention $\zeta(\mathbf x)=\int d^3k\,e^{i\mathbf k\cdot\mathbf x}\zeta(\mathbf k)/(2\pi)^3$. Statistical homogeneity and isotropy define the <primordial bispectrum> by
$$
\boxed{\langle\zeta(\mathbf k_1)\zeta(\mathbf k_2)\zeta(\mathbf k_3)\rangle_c=(2\pi)^3\delta^{(3)}(\mathbf k_1+\mathbf k_2+\mathbf k_3)B(k_1,k_2,k_3).}
$$
The <Dirac delta function> enforces <momentum conservation>; the three magnitudes describe a closed triangle. Set $\alpha=3f_{\mathrm{NL}}/5$. For the <centered quadratic Gaussian transformation>, its quadratic Fourier term is
$$
\alpha\int\frac{d^3q}{(2\pi)^3}G(\mathbf q)G(\mathbf k-\mathbf q)-\alpha(2\pi)^3\delta^{(3)}(\mathbf k)\langle G^2\rangle.
$$
At first order in $\alpha$, choose the quadratic field at one of the three external positions. For the first position, <Wick contractions> pair the two internal fields with the two remaining external fields in two ways. The self-pairing is precisely removed by the subtracted mean. The surviving contribution is $2\alpha P(k_2)P(k_3)$ times the momentum delta. Adding the other positions gives the tree-level <local-type primordial non-Gaussianity> result
$$
\boxed{B^{\mathrm{loc}}=\frac65f_{\mathrm{NL}}[P(k_1)P(k_2)+P(k_2)P(k_3)+P(k_3)P(k_1)].}
$$
\b[This expression is at leading order in $f_{\mathrm{NL}}$.] The exact quadratic map also has a <loop correction to the local primordial bispectrum>, from three quadratic vertices:
$$
B_{\mathrm{loop}}=8\alpha^3\int\frac{d^3q}{(2\pi)^3}P(q)P(|\mathbf k_1-\mathbf q|)P(|\mathbf k_2+\mathbf q|).
$$
There is no quadratic-in-$\alpha$ bispectrum term, because it would involve five centered Gaussian fields. A regulator may be needed for this higher-order integral. The displayed leading formula uses the Gaussian power $P$, consistently neglecting its order-$f_{\mathrm{NL}}^2$ correction.