Expectation of exponential of a function of independent Rademacher r.v.'s involving the error function Let $Z,Z'\in\{-1,1\}^n$ be two independent vectors of i.i.d. Rademacher r.v.'s, where $1\leq n \leq d$ are two integers ($d\gg 1$). I am trying to get an upper bound on
$$
\mathbb{E}_{ZZ'}\left[ \exp\left({k\left(e^{\sum_{j=1}^n \mathrm{Erf}(\varepsilon Z_j /\sqrt{d})} -1\right)\left(e^{\sum_{j'=1}^n \mathrm{Erf}(\varepsilon Z'_{j'} /\sqrt{d})} -1\right)}\right) \right] \tag{$\dagger$}
$$
where $t\geq 1$, $\varepsilon \in(0,1]$, $k\geq 1$ is an integer, and $\mathrm{Erf}$ is the error function. In particular, I am interested in the dependence on $n$: since
$$
\mathrm{Erf}(x) \operatorname*{\sim}_{x\to\ 0} \frac{2}{\sqrt{\pi}}x
$$
I expect that, for $n = n(d) \ll \sqrt{d}$, $(\dagger)$ should be roughly
$$
\mathbb{E}_{ZZ'}\left[ \exp\left({\frac{4}{\pi}\frac{\varepsilon^2 k}{d}\sum_{j,j'} Z_jZ'_{j'}}\right) \right] \leq e^{C \frac{\varepsilon^4 k^2}{d^2}\cdot n^2}
$$
(the last bound for $d$ big enough). But I'm unclear if the above qualitative behavior still holds for $n \gg \sqrt{d}$: do we still get that $n^2$ dependence in the exponent?

Note: I am actually interested in the quantity
$$
\mathbb{E}_{ZZ'}\left[ \exp\left({k \left(\prod_{j=1}^n  (1+\mathrm{Erf}(\varepsilon Z_j /\sqrt{d}))-1\right)\left(\prod_{j'=1}^n  (1+\mathrm{Erf}(\varepsilon Z'_{j'} /\sqrt{d}))-1\right)}\right) \right] \tag{$\ddagger$}
$$
(same questions), but $(\dagger)$ is similar and hopefully a bit cleaner to analyze.
 A: (This was supposed to be a comment, but I wrote in the Answer box because I don't have the privilege to comment. Apologies for the confusion.)
Every $Z_i = \pm 1$, so you have $\mathrm{Erf}(Z_i \delta) = Z_i \mathrm{Erf}(\delta) \equiv Z_i \tau$, and the products become $(1 + \tau)^b (1 - \tau)^{n - b}$ where $b$ follows the Binomial distribution $B(n, 1/2)$, and the same for the primed quantities. The expectation is then $$\mathbb{E}_{bb'}\left[\exp\left({k \left((1 + \tau)^b (1 - \tau)^{n-b}-1\right)\left((1 + \tau)^{b'} (1 - \tau)^{n-b'}))-1\right)}\right) \right] .$$
Perhaps this allows you to get the kind of bound you are looking for. You have $0 < \tau < 1$, so you can bracket $(1 + \tau)^b(1 - \tau)^{n-b}$ by $(1-\tau)^n$ and $(1+\tau)^n$. Without having too good an idea about how to proceed, I would naively expect the expectation to behave like $\exp\bigl(k \, (1+\tau)^{2n})\bigr)$ as $n$ increases. A quick Mathematica experiment with $k = 1.23$, $\delta = 0.1$ seems to bear that out, but there could be all kinds of things going wrong here.

(Since you now have only two r.v.s to contend with, you can compute the expectation explicitly to get some additional insight, at least until the arithmetic overflows for large $n$.)
