UMVUE for $g(\theta)=\theta^2$ of Poisson random variables Let $X_1,...X_n$i.i.d.~$Pois(\theta)$ with unknown $\theta>0$, I want to find the UMVUE for $g(\theta)=\theta^2$.
I know that $T(x)=\Sigma_{i=1}^{n}x_i$ is complete and sufficient for $\theta>0$.
$\mathbf{E}_\theta[T(x)]=\mathbf{E}_\theta[\Sigma_{i=1}^{n}x_i]=\Sigma_{i=1}^{n}\mathbf{E}_\theta[x_i]=n\theta$
$\mathbf{E}_\theta[T(x)^2]=\mathbf{Var}[T(x)]+\mathbf{E}_\theta[T(x)]^2=n\theta+n^2\theta^2$
Hence $\frac{\mathbf{E}_\theta[T(x)^2]-\mathbf{E}_\theta[T(x)]}{n^2}=\theta^2=g(\theta)$
So is $\frac{T(x)^2-T(x)}{n^2}$ the UMVUE for $g(\theta)=\theta^2$?

Yes, that's correct.  To make it rigorous, you can apply Lehmann–Scheffé Theorem with $Y = T(X)$ and $\varphi(Y) = \frac{1}{n^2}(Y^2 - Y)$.

After a closer look at the Wikipedia link above, it seems to confuse $\theta$ with $g(\theta)$.  To avoid such ambiguity, I am quoting Theorem 1.1 in Theory of Point Estimation by E. Lehmann and G. Casella below.

Let $X$ be distributed according to a distribution in
$\mathcal{P} = \{P_\theta, \theta \in \Omega\}$, and suppose $T$ is a complete sufficient statistic for $\mathcal{P}$.
(a) For every $U$-estimable function $g(\theta)$, there exists an
unbiased estimator that uniformly minimizes the risk for any loss
function $L(\theta, d)$ which is convex in its second argument;
therefore, this estimator in particular is UMVU.
(b) The UMVU estimator of (a) is the unique unbiased estimator which
is a function of $T$; it is the unique unbiased estimator with minimum
risk, provided its risk is finite and $L$ is strictly convex in $d$.

The author also gave the following procedure for determining the UMVUE, which you can directly harness for your problem:

If $T$ is a complete sufficient statistic, the UMVU estimator of any $U$-estimable function $g(\theta)$ is uniquely determined by the set of equations
\begin{align}
E_\theta\delta(T) = g(\theta) \quad \text{for all } \theta \in \Omega.
\end{align}

