Where do I make a mistake in the proof of the sample variance for the unbiased estimator? I am trying to prove that the sample variance formula is unbiased.
Firts, let $\mathbb{E}(X_i) = \mu, \text{Var}(X_i)=\sigma^2=\mathbb{E}(X_i^2)-\mathbb{E}(X_i)^2, \bar{X}=\frac{\sum X_i}{n}$.
$\mathbb{E}( \frac{\sum(X_i-\bar{X})^2}{n-1}) =\sigma^2 $ can be transformed as follows.
$\mathbb{E}( \sum(X_i-\bar{X})^2/(n-1))= \frac{\sum\mathbb{E}(X_i-\bar{X})^2}{n-1}\cdot\dots\cdot (1)$
$ \sum\mathbb{E}(X_i-\bar{X})^2=  \sum\mathbb{E}(X_i- \mu +\mu-\bar{X})^2 \cdot\dots\cdot (2)$
$ \mathbb{E}(X_i- \mu +\mu-\bar{X})^2 = \mathbb{E}( \underbrace{(X_i- \mu)^2}_{a} +\underbrace{(\mu-\bar{X})^2}_{b}  +2\underbrace{(X_i- \mu)(\mu-\bar{X})}_{c}) $
$a: \mathbb{E}((X_i- \mu)^2 = \sigma^2 $
$b: \mathbb{E}((\mu-\bar{X})^2 = \frac{\sigma^2}{n} $
$c$ is arranged as follows.
$ \mathbb{E}(X_i- \mu)(\mu-\bar{X}) = \mathbb{E}(X_i\mu-X_i\bar{X} -\mu^2+\mu\bar{X})$
$ \mathbb{E}(X_i\mu-X_i\bar{X} -\mu^2+\mu\bar{X})=\mathbb{E}(-X_i\bar{X} +\mu\bar{X})$
$ \because \mathbb{E}(X_i\mu) = \mu^2$
$ \mathbb{E}(-X_i\bar{X}) =- \mathbb{E}(X_i\sum{X_i}/n)=-\frac{\sum\mathbb{E}(X_i^2)}{n} = -\mathbb{E}(X_i^2)$
$ \mathbb{E}(\mu\bar{X})= \mu\mathbb{E}(\sum{X_i}/n)= \mu\frac{\sum{\mathbb{E}(X_i)}}{n}=\mathbb{E}(X_i)^2$
Therefore,
$ \mathbb{E}(X_i\mu-X_i\bar{X} -\mu^2+\mu\bar{X})= -\mathbb{E}(X_i^2)+\mathbb{E}(X_i)^2=-\sigma^2$
Then, $(2)$ becomes $ \sigma^2+ \frac{\sigma^2}{n} -2\sigma^2=-\sigma^2+ \frac{\sigma^2}{n}=-\frac{n-1}{n}\sigma^2$, and plugging this into (1) leads to $-\sigma^2$.
Which stage is incorrect?

Something that you have not explicitly stated in your question, but which I assume that you are implicitly assuming, is that your $X_i$'s are independent of one another (or at least uncorrelated). In that case, You have made an error in your calculation of $\mathbb E[-X_i\bar X]$. I believe the error is due to a bad notational choice when expanding out $\bar X$. Specifically, the correct way to make this calculation is:
$$\mathbb E[-X_i\bar X] = \mathbb E\left[-X_i \frac{\sum_{j=1}^n X_j}{n}\right] = -\frac1n\mathbb E[X_i^2] - \mathbb E\left[\frac{\sum_{j\neq i} X_iX_j}{n}\right] = - \frac1n \mathbb E[X_i^2] - \frac{(n-1)\mu^2}{n}$$
where in the last line, I use the fact that $X_i$ and $X_j$ are uncorrelated, so $\mathbb E[X_i X_j] = \mathbb E[X_i]\mathbb E[X_j] = \mu^2$. In that case, we have $c = \mu^2 - \frac{1}{n}\mathbb E[X_i^2] - \frac{(n-1)\mu^2}{n} = -\frac1n\left[\mathbb E[X_i^2] - \mu^2\right] = -\frac1n \sigma^2$. Thus, (2) becomes $\sigma^2 + \frac{\sigma^2}{n} - 2 \frac{\sigma^2}{n} = \frac{(n-1)\sigma^2}{n}$ as desired.
