You do not need delta method here.
$$
\sqrt{n}\left(S_n^2-\sigma^2\right) = \sqrt{n}\left(\frac1nS_n^2+\frac{n-1}n S_n^2 - \sigma^2 \right) = \frac{S_n^2}{\sqrt{n}} + \sqrt{n}\left(\frac {n-1}nS_n^2-\sigma^2\right).
$$
Since $S_n^2\xrightarrow{p}\sigma^2$ and $\frac1{\sqrt{n}}\to 0$, Slutsky's theorem implies the first summand tends to zero in distribution, and whence in probability. Again apply Slutsky's theorem and get the sum converges in distribution to $N(0,\text{Var}(X_1-\mu)^2)$.
Below is a back-of-envelop computation. Someone else should please check that I haven't screwed anything. Also, I'm not certain we need to impose tail bounds on $X$ (e.g sub-Gaussianity) as I've done below...
So lat $X$ be a random vector in $\mathbb R^p$ (the scalar case corresponds to $p=1$) with mean $E[X] = \mu \in \mathbb R^p$ and covariance matrix $\Sigma$ (a $p$-by-$p$ psd matrix). Let $\hat{\mu}_N := (1/N)\sum_{i=1}^NX_i$, where $X_1,\ldots,X_N$ are $N$ iid copies of $X$.
The triangle inequality, one computes
$$
\begin{split}
|E[f(X)] - f(\hat{\mu}_N)| &= |E[f(X)] - f(X) + f(X) - f(\mu) + f(\mu) - f(\hat{\mu}_N)|\\
&\le |E[f(X)] - f(X)| + |f(X) - f(\mu)| + |f(\mu)-f(\hat{\mu}_N)|
\end{split}
\tag{*}
$$
Now,
Suppose $X$ is $\sigma$-subGaussian and $f$ is $L$-Lipschitz.
Then, with probability $\ge 1 - \delta$, it holds that
$$
\begin{split}
(1/L)|f(\mu)-f(\hat{\mu}_N)| &\le \|\mu-\hat{\mu}_n\| \le \sqrt{\frac{\text{Tr}(\Sigma)}{N}} + \sqrt{\frac{2\lambda_\max(\Sigma)\log(1/\delta)}{N}}\\
&= \mathcal O(\sqrt{d\log(1/\delta)/N}).
\end{split}
$$
Thus by (*), asymptotically in $N$, we have
$$
\begin{split}
P(|E[f(X)] - f(\hat{\mu}_N)| > \epsilon) &\le P(|E[f(X)] - f(X)| > \frac{\epsilon}{3}) + P(|f(X) - f(\mu)| > \frac{\epsilon}{3})\\
&\quad\quad + P(|f(\mu)-f(\hat{\mu}_N)| > \frac{\epsilon}{3})\\
&\le A_1e^{-C_1\epsilon^2} + A_2e^{-C_2\epsilon^2} + \mathcal O(1/\sqrt{N}),
\end{split}
$$
for some positive absolute constants $A_1,A_2,A,C_1,C_2,C$.
The first term in the above bound is by concentration of Lipschitz functions of subGaussian vectors and and the second is because $|f(X) - f(\mu)| \le L\|X-\mu\|$ and using (again!) the subGaussianity of $X$.
Best Answer
You can get a bound from Markov's inequality, but that's about it. In particular, there is no hope of getting such a strong bound (which is called sub-Gaussian concentration) without much stronger assumptions than just two moments. To see why, suppose that X is a random variable for which $$\mathbb{P}(|X| > t) \leq e^{-ct^2}$$ holds for all $t>0$. Then, using the usual trick for computing moment estimates from tail estimates \begin{align*} \mathbb{E}[e^{\frac{c}{2}X^2}] &= 1 + \int_0^\infty c t e^{\frac{c}{2}t^2}\mathbb{P}(|X|>t)dt \\ &\leq 1+ \int_0^\infty c t e^{\frac{c}{2}t^2}e^{-ct^2} dt = 1+ \int_0^\infty c te^{-\frac{c}{2}t^2} dt <\infty. \end{align*}
Now, let $c>0$ be any number and suppose that the $X_i$ are i.i.d. Exponential(1) random variables. Note that these random variables are much nicer than typical random variables with 2 moments -- they have finite exponential moments for powers less than 1. Even with these nice random variables, this fails already when $n=2$: \begin{align*} \mathbb{E}[e^{\frac{c}{2} (S_2 - 1)^2}] &= \int_{0}^\infty \int_0^\infty e^{\frac{c}{2}\left[\frac{1}{2}(x^2+y^2) - (\frac{1}{2}(x+y))^2 -1 \right]^2} e^{-(x+y)}dxdy = \infty. \end{align*} The same computation works for all $n \geq 2$.
In general, concentration bounds are basically the same thing as moment bounds. You can't hope for much better concentration than you put in with your hypotheses. So what does the hypothesis of having two finite moments buy you in this case? By Markov's inequality, there is a bound of the form
$$ \mathbb{P}(|S_n - Var(X_1)|>t) \leq \frac{\mathbb{E}[|S_n - Var(X_1)|]}{t}. $$
I'll sketch how to show that you cannot significantly improve on this. The same argument as above shows that if there exists constants $c,C,\alpha>0$ so that for all $t>c$ $$ \mathbb{P}(|X|>t) \leq C t^{-\alpha} \implies \mathbb{E}[|X|^\beta] <\infty $$ whenever $\beta < \alpha$. So now let's consider for $\epsilon>0$ i.i.d. random variables with density \begin{align*} f_{X_i}(t) = \begin{cases} \frac{1}{2+\epsilon}t^{-(3+\epsilon)} & t>1 \\ 0 & t \leq 1 \end{cases} \end{align*} Call $\sigma^2 = Var(X_1)$. We can compute that \begin{align*} \mathbb{E}[|S_2-\sigma^2|^{1+\epsilon}] &= \int_1^\infty \int_1^\infty \left|\frac{1}{2}(x^2+y^2) - (\frac{x}{2}+\frac{y}{2})^2-\sigma^2\right|^{1+\epsilon}\frac{1}{x^{2+\epsilon}y^{2+\epsilon}} \frac{1}{(1+\epsilon)^2}dxdy \\ &= \int_1^\infty \int_1^\infty \left|\frac{(x-y)^2}{4} -\sigma^2\right|^{1+\epsilon}\frac{1}{x^{3+\epsilon}y^{3+\epsilon}} \frac{1}{(2+\epsilon)^2}dxdy =\infty. \end{align*} To see that this is infinite, note that the inner integral can be seen to be infinite by limit comparison with $$ \frac{x^{2(1+\epsilon)}}{x^{3+\epsilon}} = x^{-(1-\epsilon)}. $$ For $\alpha>1$, this rules out any bound of the form \begin{align*} \mathbb{P}(|S_n-Var(X_1)| > t) \leq \frac{C}{t^{\alpha}} \end{align*} holding for all sufficiently large $t$ under these hypotheses. More complicated examples can be constructed to rule out weaker bounds as well.
Let me close by commenting on a common point of confusion. Convergence in distribution tells you essentially nothing about the regularity of the pre-limit objects. For example, let $X$ and $Y$ be independent N(0,1) and standard Cauchy random variables respectively and set $Z_n = (1-e^{-n})X + e^{-n} Y.$ $Z_n$ converges exponentially fast to a $N(0,1)$ random variable surely (never mind in distribution), but it doesn't have a finite mean.