I don't work directly through your derivation, but provide a more general formulation below.
For a more general formulation, let your regression model be $Y = X\beta + \epsilon$, $P_X = X(X^\prime X)^{-1} X^\prime$, and $M_X = I_N - P_X$ ($I_N$ is a $N\times N$ identity matrix). $X$ is $N\times K$ and of full column rank. We assume homoskedasticity and no serial correlation.
We show that $\hat{\sigma}^2$ is unbiased:
$$\begin{align*}
\mathbb{E}\left[\frac{\hat{\epsilon}^\prime \hat{\epsilon}}{N - K}\mid X\right] &= \mathbb{E}\left[\frac{\epsilon^\prime M^\prime M \epsilon}{N - K}\mid X\right] \\
&= \mathbb{E}\left[\frac{\epsilon^\prime M \epsilon}{N - K}\mid X\right] \\
&= \frac{\sum_{i=1}^N{\sum_{j=1}^N{m_{ji}\mathbb{E}[\epsilon_i\epsilon_j\mid X]}}}{N - K} \\
&= \frac{\sum_{i=1}^N{m_{ii}\sigma^2}}{N - K} \\
&= \frac{\sigma^2\mathop{\text{tr}}(M)}{N - K} \\
\end{align*}$$
$$\begin{align*}
\text{tr}(M) &= \text{tr}(I_N - P_X) \\
&= \text{tr}(I_N) - \text{tr}(P_X) \\
&= N - \text{tr}\left(X\left(X^\prime X\right)^{-1}X^\prime\right) \\
&= N - \text{tr}\left(\left(X^\prime X\right)^{-1}X^\prime X\right) \\
&= N - \text{tr}(I_{K}) = N - K \\
\Longrightarrow \mathbb{E}\left[\frac{\hat{\epsilon}^\prime \hat{\epsilon}}{N - K}\mid X\right] &= \frac{\sigma^2 (N-K)}{(N-K)} = \sigma^2.
\end{align*}$$
Update:
Consider the estimator
$$\hat 0 = \bar{X} - cS$$
where $c$ is given in your post. This is is an unbiased estimator of $0$ and will clearly be correlated with the estimator given below (for any value of $a$).
Theorem 6.2.25 from C&B shows how to find complete sufficient statistics for the Exponential family so long as $$\{(w_1(\theta), \cdots w_k(\theta)\}$$ contains an open set in $\mathbb R^k$. Unfortunately this distribution yields $w_1(\theta) = \theta^{-2}$ and $w_2(\theta) = \theta^{-1}$ which does NOT form an open set in $R^2$ (since $w_1(\theta) = w_2(\theta)^2$). It is because of this that the statistic $(\bar{X}, S^2)$ is not complete for $\theta$, and it is for the same reason that we can construct an unbiased estimator of $0$ that will be correlated with any unbiased estimator of $\theta$ that is based on the sufficient statistics.
Another Update:
From here, the argument is constructive. It must be the case that there exists another unbiased estimator $\tilde\theta$ such that $Var(\tilde\theta) < Var(\hat\theta)$ for at least one $\theta \in \Theta$.
Proof: Let suppose that $E(\hat\theta) = \theta$, $E(\hat 0) = 0$ and $Cov(\hat\theta, \hat 0) < 0$ (for some value of $\theta$). Consider a new estimator
$$\tilde\theta = \hat\theta + b\hat0$$
This estimator is clearly unbiased with variance
$$Var(\tilde\theta) = Var(\hat\theta) + b^2Var(\hat0) + 2bCov(\hat\theta,\hat0)$$
Let $M(\theta) = \frac{-2Cov(\hat\theta, \hat0)}{Var(\hat0)}$.
By assumption, there must exist a $\theta_0$ such that $M(\theta_0) > 0$. If we choose $b \in (0, M(\theta_0))$, then $Var(\tilde\theta) < Var(\hat\theta)$ at $\theta_0$. Therefore $\hat\theta$ cannot be the UMVUE. $\quad \square$
In summary: The fact that $\hat\theta$ is correlated with $\hat0$ (for any choice of $a$) implies that we can construct a new estimator which is better than $\hat\theta$ for at least one point $\theta_0$, violating the uniformity of $\hat\theta$ claim for best unbiasedness.
Let's look at your idea of linear combinations more closely.
$$\hat\theta = a \bar X + (1-a)cS$$
As you point out, $\hat\theta$ is a reasonable estimator since it is based on Sufficient (albeit not complete) statistics. Clearly, this estimator is unbiased, so to compute the MSE we need only compute the variance.
\begin{align*}
MSE(\hat\theta) &= a^2 Var(\bar{X}) + (1-a)^2 c^2 Var(S) \\
&= \frac{a^2\theta^2}{n} + (1-a)^2 c^2 \left[E(S^2) - E(S)^2\right] \\
&= \frac{a^2\theta^2}{n} + (1-a)^2 c^2 \left[\theta^2 - \theta^2/c^2\right] \\
&= \theta^2\left[\frac{a^2}{n} + (1-a)^2(c^2 - 1)\right]
\end{align*}
By differentiating, we can find the "optimal $a$" for a given sample size $n$.
$$a_{opt}(n) = \frac{c^2 - 1}{1/n + c^2 - 1}$$
where
$$c^2 = \frac{n-1}{2}\left(\frac{\Gamma((n-1)/2)}{\Gamma(n/2)}\right)^2$$
A plot of this optimal choice of $a$ is given below.
It is somewhat interesting to note that as $n\rightarrow \infty$, we have $a_{opt}\rightarrow \frac{1}{3}$ (confirmed via Wolframalpha).
While there is no guarantee that this is the UMVUE, this estimator is the minimum variance estimator of all unbiased linear combinations of the sufficient statistics.
Best Answer
Say $Q$ is unbiased for $\theta^2$, i.e. $E(Q) = \theta^2$, then because of Jensen's inequality,
$$\sqrt{ E(Q) } = \theta < E \left( \sqrt{Q} \right)$$
So $\sqrt{Q}$ is biased high, i.e. it will overestimate $\theta$ on average.
Note: This is a strict inequality (i.e. $<$ not $\leq$) because $Q$ is not a degenerate random variable and square root is not an affine transformation.