I'm interested in finding tail bound for $sum_{i=1}^k X_i^2 Y_i^2$, where $X_i$ and $Y_i$ are independent standard normal random variables.
It should be roughly as tight as the standard Chernoff bound, something like $e^{-Omega(ksqrt t)}$ would be nice.

My first instinct was to look at the mgf. $exp(t X^2 Y^2)$, but naturally it doesn't exist. I looked at $exp(i t X^2 Y^2) = e^{i/(8t)} K_0(i/(8t))/sqrt{pi i t}$, but I don't know how to get a tail bound using the characteristic function. I also considered moment bounds. We have $E(X^2 Y^2)^k=2^{2k}Gamma(k+1/2)^2/pile2(2k/e)^{2k}$, but to get a tail bound I need $E(sum_i X_i^2 Y_i^2)^k$, which is of course a lot harder to estimate. I also considered the Cauchy Schwarz bound: $E(sum_i X_i^2 Y_i^2)^kle E(sum_i X_i^4)^{k/2}(sum_i Y_i^4)^{k/2}=E(sum_i X_i^4)^{k}$, but even those moments seem pretty involved.

We have that $Pr[sum_{i=1}^kX_iY_ige tk]le expleft(frac{-t^2k}{2+t}right)$ by Chernoff bounds, which is close to gaussian at least for small $t$. Perhaps we might also expect that $sum X^2_iY^2_i$ is close to Chi-Squared for small $t$?

Does anyone know if there is a standard bound for this sum? Or if one of my approaches might be workable?

I found that a good way to handle this sum is using Bernstein's inequality:

Let $X_1, ldots, X_n$ be independent zero-mean random variables. Suppose that $|X_i|leq M$ almost surely, for all $i$. Then, for all positive
$t$,

$$Pr left (sum_{i=1}^n X_i ge t right ) leq exp left(
-frac{t^2/2}{sum mathbb{E} X^2+ Mt/3} right).$$

We have
$$
Pr[X^2Y^2ge t]
le E[(XY)^{2p}]t^{-p}
le 2(2p/e)^{2p}t^{-p}
le 2 e^{-sqrt t},
$$
taking $p=sqrt{t}/2$. Hence we can take a union bound over all $XY$ to get:
$$
begin{align}
Prleft(
sum_k X_k^2Y_k^2ge (1+epsilon)k
right)
&le
exp left(
-frac{kepsilon^2/2}{9+ Mepsilon/3} right)
+2k exp(-sqrt M)le2delta
end{align}
$$
When taking $M=log^2((2 k)/delta)$ and
$$begin{align}t
&=
frac{18 log left(frac{1}{delta}right)}{epsilon^2}+frac{2 log left(frac{1}{delta}right) log^2 left(frac{2 k}{d}right)}{3 epsilon}
=O(epsilon^{-2}log1/delta+epsilon^{-1}log^32k/delta).
end{align}$$

This matches exactly what we would expect from the central limit theory in the first term, and nearly bound from a "single large term" in the second term.

The only remaining question is whether we can get rid of the third power, to just have $epsilon^{-1}log^2(k/delta)$ in the second term. I don't know how to do that though.

Update

The previous approach lost a factor $log1/delta$ in $t$.
This can be avoided by the following neat trick: Instead of using Bernstein's inequality, just use $(XY)^2le M^{1/2}|XY|le sqrt M (X^2+Y^2)/2$.

That gives us
$$begin{align}
Pr[sum_iX_i^2Y_i^2 ge (1+epsilon)k mid X^2Y^2le M]
&le expleft(lambda k X^2Y^2kright)/expleft(lambda(1+epsilon)right)
\&le expleft(lambda k sqrt M (X^2+Y^2)/2right)/expleft(lambda(1+epsilon)kright)
\&= frac{1}{1-lambda ksqrt M}/expleft(lambda(1+epsilon)kright)
\&=frac{(1+epsilon)k}{sqrt M}exp(1 - frac{(1+epsilon)k}{sqrt M})
\&=frac{(1+epsilon)k}{log(2k/delta)}exp(1 - frac{(1+epsilon)k}{log(2k/delta)})
%\&=sqrt{(1+epsilon)k}exp(1 - sqrt{(1+epsilon)k}).
end{align}$$

Taking $lambda=((1+epsilon)k-M)/(M (1+epsilon)k)$ and $M=(1+epsilon)k$.
The union bound over from the single elements is then $+2kexp(-sqrt{(1+epsilon)k})$.

We see that taking $kapproxepsilon^{-2}log1/delta + epsilon^{-1}(log1/delta)^2$ now suffices.

The method is nearly as versatile as the Bernstein approach, and it gives the optimal values.