Given independent Gaussian random vectors $a,c \sim \mathcal{N}(0,I_p)$ and $b \sim \mathcal{N}(0,\frac1p I_p)$, I would like to compute the following expectation $$\mathbb{E}_{a,b,c} f(a^T b) f(b^T c) = \int_{\mathbb{R}^p} \int_{\mathbb{R}^p} \int_{\mathbb{R}^p} f(a^T b) f(b^T c) (2\pi)^{-\frac{3p}{2}} e^{-\frac12 (\|a\|^2 + \frac1p \|b\|^2 + \|c\|^2)} da\ db\ dc$$ for arbitrary function $f$. Some existing results (via the so-called Orthogonal polynomials) suggest that this expectation may depend on the following two quantities $$d_1 = \int_{\mathbb{R}} \frac1{\sqrt{2\pi}} f^2(x) e^{-\frac{x^2}2} dx$$ and $$d_2 = \int_{\mathbb{R}} \frac1{\sqrt{2\pi}} f^\prime(x) e^{-\frac{x^2}2} dx$$ for example, if $a=c$, it resorts to compute $\mathbb{E}_{a,b} f^2(a^T b)$. Since $a,b$ have independent normal entries, the variable $a^T b = \sum_{i=1}^p a_i b_i \equiv \sum_{i=1}^p z_i$ with $z_i \sim \mathcal{N}(0,\frac1{\sqrt{p}})$ and therefore $a^T b \sim \mathcal{N}(0,1)$ with is connected to $d_1$.
Any one has some idea?
I don't think your last paragraph is correct. $a^Tb = \frac{1}{\sqrt{p}}a^Ta$ with $a \sim N(0,I_p)$. This gives us:
\begin{equation} a^Tb = \frac{1}{\sqrt{p}}a^Ta = \frac{1}{\sqrt{p}}\sum_{i=1}^pa^2_i = \frac{1}{\sqrt{p}}\sum_{i=1}^pZ_i^2 \sim \frac{1}{\sqrt{p}}\chi^2_p \overset{d}{=}\textrm{Gamma}(\frac{p}{2},\frac{2}{\sqrt{p}}) \end{equation}
By symmetry of the real inner product $a^Tb = b^Tc$.
Specifying $a^Tb = b^Tc = X \sim \textrm{Gamma}(\frac{p}{2},\frac{2}{\sqrt{p}})$. We have that your expectation is:
\begin{equation} \int_{0}^\infty f(x)^2 \frac{(\frac{2}{\sqrt{p}})^{\frac{p}{2}}2x^{\frac{p}{2}-1}}{\Gamma(\frac{p}{2})}\textrm{Exp}(-\frac{2x}{\sqrt{p}})dx \end{equation}
Which is:
\begin{equation} E_X[f(X)^2], \hspace{2mm} X \sim \textrm{Gamma}(\frac{p}{2},\frac{2}{\sqrt{p}}) \end{equation}
This is the expectation of the second moment of $f(X)$ under the given Gamma distribution.