Apologies for using screenshots, I just don't have time to remember and look up all of the formatting tips.

Alpha is known. Taking the expected value of the MLE I get :
My first question is what am I doing wrong as I've highlighted in red? My professor is getting a different integral than me:
I asked him and all he mentioned was that it was because summation of xi was in the denominator, which still doesn't make sense how he magically came up with 1/t instead of 1/summation of xi!
My second question. Even if taking his integral and solving for it using integration by parts I keep getting an integral that requires repeat integration by parts! How can I solve this integral?



What you wrote does not make sense because the $X_i$ are random variables. If I asked you to compute for some continuous random variable $Y$ with density $f_Y(y)$ the expectation $\operatorname{E}[Y^2]$, what would you write? You would write $$\operatorname{E}[Y^2] = \int_{y \in \Omega} y^2 f_Y(y) \, dy,$$ where $\Omega$ is the support of $Y$. You would not write $$\operatorname{E}[Y^2] = \int_{y \in \Omega} Y^2 f_Y(y) \, dy.$$
Now, recall that if $X_1, X_2, \ldots, X_n$ are independent and identically distributed Gamma random variables with common shape parameter $\alpha$ and common rate parameter $\lambda$, i.e. $$X_i \sim \operatorname{Gamma}(\alpha, \lambda), \\ f_X (x) = \frac{\lambda^\alpha x^{\alpha - 1} e^{-\lambda x}}{\Gamma(\alpha)}, \quad x > 0,$$ then the sum $$T = \sum_{i=1}^n X_i \sim \operatorname{Gamma}(n\alpha, \lambda)$$ is also gamma distributed, but with shape $n\alpha$ and rate $\lambda$. Thus the density of $T$ is $$f_T(t) = \frac{\lambda^{n\alpha} t^{n\alpha-1} e^{-\lambda t}}{\Gamma(n\alpha)}.$$ The expectation of $1/T$ is then $$\operatorname{E}[1/T] = \int_{t=0}^\infty \frac{1}{t} f_T(t) \, dt = \int_{t=0}^\infty \frac{\lambda^{n \alpha} t^{n\alpha - 2} e^{-\lambda t}}{\Gamma(n \alpha)} \, dt.$$ To evaluate this integral, we employ a little trick: we write the integrand as $$\frac{\lambda^{n \alpha} t^{n\alpha - 2} e^{-\lambda t}}{\Gamma(n \alpha)} = \frac{\lambda}{n\alpha - 1} \frac{\lambda^{n\alpha - 1} t^{(n\alpha-1)-1} e^{-\lambda t}}{\Gamma(n \alpha - 1)},$$ where we have employed the property that $\Gamma(z) = (z-1)\Gamma(z-1)$. Then the second factor $$\frac{\lambda^{n\alpha - 1} t^{(n\alpha-1)-1} e^{-\lambda t}}{\Gamma(n \alpha - 1)}$$ is the density of a gamma-distributed random variable with shape $n\alpha - 1$ and rate $\lambda$, thus integrates to $1$ over its support. It follows that $$\operatorname{E}[1/T] = \frac{\lambda}{n\alpha-1}$$ and $$\operatorname{E}[n\alpha/T] = \frac{n\alpha}{n\alpha-1} \lambda.$$ Thus, this estimator is in fact biased for the rate parameter $\lambda$. The bias is especially large when $\alpha$ is small, something that is easily seen by performing a simulation.
I honestly do not know how your professor could have made such a mistake. It is well-known that the MLE for the rate parameter of a gamma distributed sample is biased. The reciprocal of the sample mean is inverse gamma distributed. When I look at the solution you furnished, it is clear that there is an error: the addition of iid gamma random variables does not change the rate, but the shape, in as much as the sum of iid exponential random variables makes a gamma distribution with shape equal to the number of exponential variables added together; e.g. if $$Y_1, Y_2, \ldots, Y_\alpha \sim \operatorname{Exponential}(\lambda), \\ f_Y (y) = \lambda e^{-\lambda}, \quad y > 0,$$ then $$X = Y_1 + \cdots + Y_\alpha \sim \operatorname{Gamma}(\alpha, \lambda)$$ and in turn, $$X_1 + X_2 + \cdots + X_n$$ is simply the sum of $n\alpha$ iid exponential random variables hence is gamma with shape $n\alpha$, not rate $n\lambda$.