I posted this question previously to MathSE. However, I have still not solved it, so lets try to ask it here. When doing some calculations with spin-foam models for 3d quantum gravity for some research project, I quite generically stumble over integrals of the following form

$$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(g^{n})$$

where $n\in\mathbb{N}_{>0}$, $\mathrm{d}g$ denotes the normalized Haar measure and where $\delta(g)$ denotes the $\mathrm{SU}(2)$-delta function, i.e. the distribution defined via

$$\langle\delta,f\rangle:=\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(gh^{-1})f(g)=f(h).$$

Now I was told that it is possible to calculate these integrals explicitely and what one gets is formally something like the sum over the $n$-roots of unity, which is also what I would expect, since $\delta(g^{n})$ in the end just says that $g^{n}=1$. So, my approach to this question was the following: It is well-known, using the Theorem of Peter-Weyl, that the delta-function can formally be written as the sum

$$\delta(g)=\sum_{j\in\mathbb{N}/2}(2j+1)\chi^{j}(g),$$

where $\chi^{j}$ are the characters of the spin-$j$-representation, i.e. the unique (up to unitary equivalence) irreducible unitary representation of dimension $(2j+1)$ of $\mathrm{SU}(2)$. Now, as a next step, I parametrize the $\mathrm{SU}(2)$-group elements as $$g=e^{i\varphi\vec{n}\cdot\vec{\sigma}}$$ where $\vec{n}$ is a unit vector and $\varphi\in [0,2\pi]$. In this parametrization, the Haar measure is given by $$\mathrm{d}g=\frac{1}{\pi}\mathrm{sin}(\varphi)^{2}\mathrm{d}\varphi\mathrm{d}^{3}\vec{n}.$$ Since $\chi^{j}$ are class functions, we can always rotate to the $z$-direction, which means that we have that $$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\chi^{j}(g)=\frac{1}{\pi}\int_{0}^{2\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\,\chi^{j}(e^{i\varphi\sigma_{z}})=\frac{2}{\pi}\int_{0}^{\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\,\chi^{j}(e^{i\varphi\sigma_{z}}).$$ Hence, let us look at integrals of the type $$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\chi^{j}(g^{n})=\frac{2}{\pi}\int_{0}^{\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\,\chi^{j}(e^{in\varphi\sigma_{z}})=\frac{2}{\pi}\int_{0}^{\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\,\frac{\mathrm{sin}((2j+1)n\varphi)}{\mathrm{sin}(n\varphi)}$$ where in the last step we used that the characters are explicitly given by $$\chi^{j}(g(\varphi))=\frac{\mathrm{sin}((2j+1)\varphi)}{\mathrm{sin}(\varphi)}.$$

Now, using Mathematica, I got the following results:

Lets consider the following integral for all $j\in\mathbb{N}/2$ and for all $n\in\mathbb{N}$ $$\mathcal{I}_{j,n}=\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\chi^{j}(g^{n})=\frac{2}{\pi}\int_{0}^{\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\,\frac{\mathrm{sin}((2j+1)n\varphi)}{\mathrm{sin}(n\varphi)}$$ For $n=1$, the integral is given by $$\mathcal{I}_{j,1}=\delta_{j0}=\begin{cases}1 &\text{if $j=0$}\\0 > &\text{if $j\neq 0$}\end{cases}$$ and for $n=2$ by $$\mathcal{I}_{j,2}=(-1)^{2j}=\begin{cases}1 &\text{if > $j\in\mathbb{N}$}\\-1 &\text{if > $j\in\mathbb{N}_{\mathrm{odd}}/2$}\end{cases}$$ and for $n=3$ by $$\mathcal{I}_{j,3}=\frac{1+(-1)^{2j}}{2}=\begin{cases}1 &\text{if > $j\in\mathbb{N}$}\\0 &\text{if > $j\in\mathbb{N}_{\mathrm{odd}}/2$}\end{cases}$$

We could continue with higher $n$'s, but let us stick to the cases $n=1,2,3$.

Now, we are interested in integrals of the following type:

$$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(g^{n})=\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\sum_{j\in\mathbb{N}_{0}/2}(2j+1)\chi^{j}(g^{n})=\frac{2}{\pi}\int_{0}^{\pi}\,\mathrm{d}\varphi\,\mathrm{sin}(\varphi)^{2}\sum_{j\in\mathbb{N}_{0}/2}(2j+1)\chi^{j}(e^{in\varphi\sigma_{z}})=\sum_{\mathbb{N}_{0}/2}(2j+1)\mathcal{I}_{j,n}$$

As a quick consistency check, let us see what happens for the case $n=1$:

$$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(g)=\sum_{\mathbb{N}_{0}/2}(2j+1)\mathcal{I}_{j,1}=\sum_{\mathbb{N}_{0}/2}(2j+1)\delta_{0j}=1$$

as it should. However, for $n=2,3$, it seems that we get an undefined or divergent series, as can be seen in the following calculation:

$$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(g^2)=\sum_{\mathbb{N}_{0}/2}(2j+1)\mathcal{I}_{j,2}=\sum_{\mathbb{N}_{0}/2}(2j+1)(-1)^{2j}$$ $$\int_{\mathrm{SU}(2)}\,\mathrm{d}g\,\delta(g^3)=\sum_{\mathbb{N}_{0}/2}(2j+1)\mathcal{I}_{j,3}=\sum_{\mathbb{N}_{0}/2}(2j+1)\frac{1+(-1)^{2j}}{2}$$

Of course, since the delta-function should more precisely be viewed as a distribution, we should also add a test function. However, decomposing the test function also using the Theorem of Peter-Weyl, this should not make any difference.

Does anyone know where my error is? Or maybe some alternative method to calculate the integrals?