For a concentration bound, we can decompose $Y$ by introducing random variables $Y_i = z_i^2 x_i$, where $x_i$ denotes the $i$-th coordinate of $X \in X_k$. Then $Y = \sum_{i=1}^n Y_i$. Assuming vectors $X \in X_k$ are chosen uniformly at random, the $Y_i$ become independent for fixed $X$ due to the independence of the $z_i$.  To apply [Bernstein's inequality][1], we must establish:

 1. Zero-Mean: Since $z_i \sim \mathcal{N}(0, \sigma^2)$, it follows that $\mathbb{E}[Y_i] = \mathbb{E}[z_i^2x_i] = 0$.
 2. Variance: Assuming each  $x_i$ is chosen with probability $P(x_i=1)=k/n$,  we have: 
     $$\text{Var}[Y_i]  = \mathbb{E}[\text{Var}[Y_i|x_i]]= \frac{k}{n} \cdot 2\sigma^4 = \frac{2k\sigma^4}{n}$$ 
 3. Since  a squared Gaussian is non-negative and has a  well-defined bound  determined by $\sigma^2$,  we can say that  $|Y_i| \leq \sigma^2$ almost surely.  

Applying Bernstein's inequality gives:
$$\mathbb{P}\left(\left|\sum_{i=1}^n Y_i - \mathbb{E}\left[\sum_{i=1}^n Y_i\right]\right| > t\right) \leq 2\exp\left(\frac{-t^2/2}{\sum_{i=1}^n \text{Var}[Y_i] + ct/3}\right)$$

Substituting our established mean, variance, and almost-sure bound, we get the concentration inequality for $Y$:
$$\mathbb{P}(|Y - \mathbb{E}[Y]| > t) \leq 2 \exp\left(-\frac{t^2/2}{n (n/k) \cdot 2\sigma^4 + t\sigma^2/3}\right)$$

A basic bound on the expectation,  $\mathbb{E}[Y] \leq  \binom{n}{k} \cdot k \sigma^2$,  is  obtained using  the [union bound][2]. The union bound helps us obtain an initial bound on  $\mathbb{E}[Y]$. Since  there are $|X_k| = \binom{n}{k}$ vectors in $X_k$,  and the squared projection onto any one can't exceed $k\sigma^2$,  we get:
$$\mathbb{E}[Y] \leq  \binom{n}{k} \cdot k \sigma^2$$

Simulations can be used to validate the concentration inequality and explore influences on bound tightness. I ran some simulations and there seemed to be support for the concentration inequality across a range of $n$ and $k$ values  The simulations indicate both the average value of $Y$ and the magnitude of the concentration bound generally increase with greater dimensionality $n$ or higher sparsity parameter $k$.

Here is some python code used for the simulation:

    import numpy as np
    
    def gen_Y(n, k, sigma, trials):
        results = []
        for _ in range(trials):
            Z = np.random.randn(n) * sigma
            X = np.random.choice([0, 1], size=n, p=[1 - k/n, k/n])
            Y = np.max(np.abs(Z @ X)**2)
            results.append(Y)
        return results
    
    def theoretical_bound(n, k, sigma, confidence_level=0.95):
        exponent = - (confidence_level * n * (n/k) * 2 * sigma**4) / 2 
        bound = np.sqrt(exponent) / (sigma**2 / 3)
        return bound
    
    if __name__ == "__main__":
        n_values = [10, 50, 100]
        k_values = [5, 20, 35] 
        sigma = 1 
        trials = 10000 
        confidence_level = 0.95  
    
        for n in n_values:
            for k in k_values:
                Y_vals = gen_Y(n, k, sigma, trials)  
                avg_Y = np.mean(Y_vals)
                tail_bound = theoretical_bound(n, k, sigma, confidence_level)
                tail_count = np.sum(np.abs(Y_vals - avg_Y) > tail_bound) 
                tail_prob = tail_count / trials 
                print(f"n: {n}, k: {k}, Avg Y: {avg_Y:.3f}, Tail Bound: {tail_bound:.3f}, Tail Count: {tail_count}, Empirical Tail Prob: {tail_prob:.4f}") 


  [1]: https://en.wikipedia.org/wiki/Bernstein_inequalities_(probability_theory)
  [2]: https://en.wikipedia.org/wiki/Boole%27s_inequality