We are given a probability space on $\Omega = [k] = \{1,...,k\}$. Let $p(i)$ be the probability assigned to $i$. These are unknown.
I need, for real $l > 0$, an unbiased estimator for $\sum_{i=1}^k p(i)^l$. That is, a function $f : \mathbb [k]^n \to [0,1]^k$ such that $E[f(X)] = \sum_{i=1}^n p(i)^l$.
In the case $l=1$, what works is letting $Y_i$ be the number of times that $i$ appears in the given sample vector, and letting $f(X) = \frac 1n[Y_1,Y_2,...,Y_n]$. Naturally, since $Y_i$ is binomially distributed, we get our result.
The idea, is to look now for random variables $Z_i : [k]^n \to [0,1]$,probably depending on $Y_i$, such that $E[Z_i] = p(i)^l$. But, $Z_i = Y_i^l$, an obvious candidate, doesn't do the job. Of course, obtaining such $Z_i$ will allow me to use $f(X) = [Z_i]$ as the desired answer.
For $n\geq l$, $$\mathbb E\left[Y_1(Y_1-1)\ldots(Y_1-l+1)\right]$$ $$=\sum_{k=l}^n k(k-1)\ldots(k-l+1) \binom{n}{k} \bigl(p(1)\bigr)^k \bigl(1-p(1)\bigr)^{n-k}$$ $$ =n(n-1)\ldots(n-l+1)\cdot \bigl(p(1)\bigr)^l\cdot\sum_{k-l=0}^{n-l} \binom{n-l}{k-l} \bigl(p(1)\bigr)^{k-l} \bigl(1-p(1)\bigr)^{(n-l)-(k-l)} $$ $$ =n(n-1)\ldots(n-l+1)\cdot \bigl(p(1)\bigr)^l. $$ So, $$\mathbb E\left[\frac{Y_1(Y_1-1)\ldots(Y_1-l+1)}{n(n-1)\ldots(n-l+1)}\right]=\bigl(p(1)\bigr)^l$$