Bounding an expectation involving i.i.d. standard Gaussians and Rademacher I have tried to bound the following quantity, but cannot get the "right" (conjectured) bound:
$$
\phi(\gamma,d,n) = -1+e^{\frac{1}{2}n\gamma^2 d}
 \mathbb{E}_{X}\left[\frac{\mathbb{E}_Z[\prod_{j=1}^n(1+\gamma\langle X^{(j)},Z\rangle)]^2}{\prod_{i=1}^d \cosh(\gamma\sum _{j=1}^nX^{(j)}_i)}\right]\tag{1}$$
where $\gamma\in(0,1]$, $d\gg 1$, $X=(X^{(j)})_{1\leq j\leq n}$ is a collection of i.i.d. standard $d$-dimensional Gaussian r.v.'s, and $Z$ is uniform on $\{-1,1\}^d$ (Rademacher) independent of the $X^{(j)}$'s.

Conjecture. Assuming $\gamma^2 d \leq 1$, as long as  $n\gamma^4 d \ll 1$ we have
  $$\phi(\gamma,d,n)\ll 1\tag{2}$$

I haven't been able to show it, even for $n=2$; there seems to be a delicate balancing act for things to exactly cancel the $e^{\frac{1}{2}n\gamma^2 d}$ factor...

Update: For (1) in the case $n=2$, my numerical experiments (a bit noisier than I had hoped) seem to be consistent with the conjecture:


*

*$\gamma$ fixed, varying $d$


*$d=5$ fixed, varying $\gamma$ (fit made with SciPy, in Python)

Update 2: For $d=1$ and $n\in\{1,2,3\}$, Mathematica could compute explicitly the expectation (though already it takes some time for $n=3$. The behavior on these few points is clearly linear wrt $n$:


 A: Something looks fishy at least with the first conjecture. Perhaps I'm again misinterpreting something, but the argument is as follows:
Consider $n=2$. Write $X^{(1)}=\frac{X+Y}{\sqrt 2}, X^{(2)}=\frac{X-Y}{\sqrt 2}$. Then $X,Y$ are independent standard Gaussians in $\mathbb R^d$. Also
$$
E_Z[(1+\gamma\langle X^{(1)},Z\rangle)(1+\gamma\langle X^{(2)},Z\rangle)]=1+\gamma^2\langle X^{(1)},X^{(2)}\rangle
\\
=1+\frac{\gamma^2}2(\|X\|^2-\|Y\|^2)
$$
Now, since $\|Y\|^2$ is a sum of squares of $d$ independent Gaussians, it deviates from any fixed number by about $\sqrt d$ with constant probability, so the expectation $E_Y[E_Z[(1+\gamma\langle X^{(1)},Z\rangle)(1+\gamma\langle X^{(2)},Z\rangle)]^2]$ is at least $c\gamma^4d$ regardless of $X$ and the whole expression you are interested in is at least 
$$
-1+c\gamma^4d\left(E[\cosh(\gamma \sqrt 2 W)]E[\frac 1{\cosh(\gamma \sqrt 2 W)}]\right)^d
$$ 
but for fixed $\gamma$ and $d\to+\infty$, this is, clearly, exponential in $d$. What am I missing this time? 
A: Partial progress: here is a proof for the analogous "Gaussian case" (i.e., $Z\sim N(\mathbf{0}_d, \mathbf{I}_d)$) and $n=2$. I was hoping to get some inspiration for it to handle the "Rademacher" $Z$ case I am interested in, though it doesn't seem to be helpful after all...

The goal is to bound an analogous quantity as (1), which for $n=2$ is
$$
\psi(\gamma,d,2) = \frac{e^{-\gamma^2 d}}{(1-2\gamma^2)^{d/2}}-2 + e^{\gamma^2d}\mathbb{E}_{XY}\left[ {\frac{\mathbb{E}_{Z}[ (1+\gamma\langle X,{Z}\rangle)(1+\gamma\langle Y,{Z}\rangle) ]^2}{ e^{\frac{\gamma^2}{2}\lVert X+Y \rVert_2^2} } }\right]
$$
(the denominator changes because it comes from $\mathbb{E}_Z[e^{\gamma \langle\sum_{j=1}^n X^{(j)},Z\rangle}]$, and that expression changes between the Gaussian and Rademacher cases. Similarly for why the additive "-1" is now a less nice expression.)
We can expand and compute the numerator (before squaring) as
\begin{align*}
  \mathbb{E}_{Z}\left[(1+\gamma\langle X, Z\rangle)(1+\gamma\langle Y, Z\rangle) \right]
  = \mathbb{E}_{Z}\left[1+\gamma\langle{X+Y,Z}\rangle + \gamma^2\langle X, Z\rangle\langle Y, Z\rangle \right]
  = 1+ \gamma^2\langle X, Y\rangle
\end{align*}
since $Z$ is Gaussian. From there, observing that $V := \frac{\lVert X+Y\rVert_2^2}{2}$ is a $\chi^2(d)$ r.v., we get
\begin{align*}
  e^{-\gamma^2d }&(1+\psi(d ,\gamma,2))\\
  &= \mathbb{E}_{XY}[{ e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} (1+ \gamma^2\langle X, Y\rangle)^2 }]\\
  &= \mathbb{E}_{V}[{ e^{-\gamma^2 V} }] + 2\gamma^2 \mathbb{E}_{XY}[{\langle X, Y\rangle e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] + \gamma^4\mathbb{E}_{XY}[{\langle X, Y\rangle^2e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] \\
  &\leq \mathbb{E}_{V}[{ e^{-\gamma^2 V} }] + 2\gamma^2 \mathbb{E}_{XY}[{\langle X, Y\rangle e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] + \gamma^4\mathbb{E}_{XY}[{\langle X, Y\rangle^2}] \\
  &= \frac{1}{(1+2\gamma^2)^{d /2}} + \gamma^2 \mathbb{E}_{XY}[{(\lVert X+Y\rVert_2^2 - \lVert X\rVert_2^2 - \lVert Y\rVert_2^2) e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] + d \gamma^4\,.
\end{align*}
The second term can be computed explicitly. Indeed, we have
\begin{align*}
  \mathbb{E}_{XY}[{\lVert X+Y\rVert_2^2e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }]
  &= 2 \mathbb{E}_{V}[{V e^{-\gamma^2V} }] = \frac{2d }{(1+2\gamma^2)^{1+d /2}} \\
  \mathbb{E}_{XY}[{\lVert X\rVert_2^2e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] 
  &= \mathbb{E}_{X}[{\lVert X\rVert_2^2 \mathbb{E}_{Y}[{ e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }] }]
  = \mathbb{E}_{X}[{\lVert X\rVert_2^2 \prod_{i=1}^d \mathbb{E}_{Y}{ e^{-\frac{\gamma^2}{2}(X_i+Y_i)^2} } }]\\
  &= {(1+\gamma^2)^{-d/2}}\mathbb{E}_{X}[{\lVert X\rVert_2^2e^{-\frac{\gamma^2}{2(1+\gamma^2)}\lVert X\rVert_2^2} }]
  = \frac{d (1+\gamma^2)}{(1+2\gamma^2)^{1+d /2}}
\end{align*}
from which the second term equals
\begin{align*}
  \mathbb{E}_{XY}[{(\lVert X+Y\rVert_2^2 - \lVert X\rVert_2^2 - \lVert Y\rVert_2^2) e^{-\frac{\gamma^2}{2}\lVert X+Y\rVert_2^2} }]
  &= - \frac{2d \gamma^2}{(1+2\gamma^2)^{1+d /2}}\,.
\end{align*}
Overall, we get
\begin{align*}
  e^{-\gamma^2d }\left(2-\frac{e^{-\gamma^2 d}}{(1-2\gamma^2)^{d/2}}+\psi(d ,\gamma,2)\right)
  &\leq \frac{1}{(1+2\gamma^2)^{d /2}}\left(1- \frac{2d \gamma^4}{1+2\gamma^2}\right) + d \gamma^4
  \leq \frac{1}{(1+2\gamma^2)^{d /2}} + d \gamma^4\,,
\end{align*}
from which
\begin{align*}
  \psi(d ,\gamma,2)
  = O( d \gamma^4 )
\end{align*}
as long as $\gamma^2d = O(1)$.
