Is using an approximation technique to estimate the parameter under a known posterior distribution necessary? Let $X_i \sim \text{IID Poisson}(\theta)$ for $i=1,2,...,n.$ The MLE of $\theta$ is written as
$$\hat{\theta}_\text{MLE}=\bar{x}.$$
Suppose $\theta$ has a probability density function (pdf) $p(\theta) \propto \theta^{-1/2}$, which is Jeffrey's prior. The resulting posterior distribution is $\theta|\mathbf{x} \sim \text{Gamma}(\sum_{i=1}^{n} x_i +\frac{1}{2},\frac{1}{n})$.  Under quadratic loss, it can be shown that the (minimum risk) Bayes estimator is
$$\hat{\theta}_\text{Bayes} = \bar{x}+\frac{1}{2n}.$$

Now, knowing the Bayes estimator in this scenario, if I were to compare the performances of the MLE and the Bayes estimator in terms of bias, can I simply simulate data from Poisson with parameter $\theta_0$ and estimate using $\hat{\theta}_\text{MLE}$ and $\hat{\theta}_\text{Bayes}$ directly? Then, repeat for say, $1000$ iterations and check the means of the estimated values of $\hat{\theta}$?  Alternatively, should I use some sort of approximation techniques such as MCMC?

You can certainly use simulation to examine the performance of estimators, but the difficulty is that you would need to choose one or more values of $\theta_0$, and the performance of the estimators might (in principle) depend on this choice.  In the present case you can get explicit distributional results for the estimators without simulation, so that is preferable.    Specifically, you have the distributions:
$$\hat{\theta}_\text{MLE} \sim \frac{1}{n} \text{Pois}(n \theta)
\quad \quad \quad \quad \quad
\hat{\theta}_\text{Bayes} \sim \frac{1}{n} \bigg[ \text{Pois}(n \theta) + \frac{1}{2} \bigg].$$
Consequently, the biases of the two estimators are:
$$\text{Bias}(\hat{\theta}_\text{MLE}, \theta) = 0
\quad \quad \quad \quad \quad
\text{Bias}(\hat{\theta}_\text{Bayes}, \theta) = \frac{1}{2n}.$$
