Trace of Symmetrized Product of Kronecker Deltas Define the array
$$
I_{\mu_1\mu_2\cdots\mu_{2n}} = \delta_{(\mu_1\mu_2}\cdots\delta_{\mu_{2n-1}\mu_{2n})}
$$
where the $(\cdots)$ indicates we should symmetrize over the indices (and divide by the order of the symmetric group, example below), and the indices run from $1$ to $d$. My goal is to compute the "trace" of this array, meaning the trace over each pair of indices:
$$
S_n = \sum_{\mu_1,\mu_3,\cdots\mu_{2n-1}=1}^d I_{\mu_1\mu_1\mu_3\mu_3\cdots\mu_{2n-1}\mu_{2n-1}}.
$$
As an example,
$$
I_{\mu_1\mu_2} = \frac{1}{2!}(\delta_{\mu_1\mu_2} + \delta_{\mu_2\mu_1}) = \delta_{\mu_1\mu_2}
$$
so $S_1 = d$. As a less trivial example,
$$
I_{\mu_1\mu_2\mu_3\mu_4} = \frac{2^3}{4!}(\delta_{\mu_1\mu_2}\delta_{\mu_3\mu_4} + \delta_{\mu_1\mu_3}\delta_{\mu_2\mu_4} + \delta_{\mu_1\mu_4}\delta_{\mu_2\mu_3})
$$
and $S_2 = \frac{1}{3}d(d + 2)$.
It's clear that $S_n$ is always a polynomial in $d$, and the identity permutation term in the symmetrization gives the highest order, $d^{n}$. With this it's possible to compute the trace by brute force on the computer and fit the polynomial coefficients to find
$$
S_3 = \frac{1}{15}d(d+2)(d + 4),\ \ \ \ S_4 = \frac{1}{105}d(d+2)(d+4)(d+6),\\ S_5 = \frac1{945}d(d+2)(d+4)(d+6)(d+8).
$$
This would seem to suggest that the correct answer is
$$
S_n = \frac{1}{(2n-1)!!}\prod_{k = 1}^{n-1}(d + 2k)
$$
where $x!! = x(x-2)(x-4)\cdots$ is the double factorial. I haven't been able to figure out how to prove that this is the correct expression for all $n$ and would appreciate any help on that front.
 A: The unsymmetrized product is invariant under the subgroup of the symmetric group $S_{2n}$ that permutes the deltas and/or swaps their indices but doesn’t change the index pairs. This has order $2^nn!$, so the result is $\frac{2^nn!}{2n!}=\frac1{(2n-1)!!}$ times a sum over the $(2n-1)!!$ partitions of $2n$ indices into $n$ pairs.
Denote the polynomial resulting from that sum by $p_n(d)$. Let’s construct $p_{n+1}(d)$ from $p_n(d)$. We need to add two more indices $2n+1$ and $2n+2$ to each partition of the first $2n$ indices into $n$ pairs. We can add them as a pair – if so, this is a pair whose indices are contracted in the trace, so that yields an additional factor of $d$, for a contribution $p_n(d)d$ to $p_{n+1}(d)$. Or we can split a pair and pair its indices with the new indices. We can do this in $2n$ ways  for each partition of the first $2n$ indices, since there are $n$ choices which pair to split and then $2$ choices which way around to form the two new pairs. The resulting partition makes the same contribution to the sum as the smaller one did, since the trace forces the two new indices to be equal, and then the deltas force the two old indices they were paired with to be equal, just like before. So that’s a contribution $2np_n(d)$ to $p_{n+1}(d)$. Overall, we have $p_{n+1}(d)=(d+2n)p_n(d)$, and then your formula follows easily by induction.
