I assume the G(s) which you give is the Laplace transform of the Green's function (please indicate next time what your symbols mean otherwise, a question can be hard to understand).
A negative damping rate is not damping but driving. Physically, it results from an external system which constantly pumps energy into your oscillator (I assume your 2nd order ode is an oscillator). In your case with $G(s) =\frac{5}{s^2-6s+16}$, the real-time Green's function is $G(t) = \frac{5}{\sqrt{7}} \sin(\sqrt{7} t) e^{3t}$, i.e., an oscillatory solution whose amplitude is exponentially growing.
By the way, why do you get that the frequency is 4 (I have $\sqrt{7}$). Do you use a different convention for the Laplace transform. My convention is $G(s) = \int_0^\infty \! dt \, G(t) e^{-s t}$.
The answer to your first question, there is no zero in the transfer function when you set $\dot{r}=0$, and this is mostly certainly due to some confusion.
When you differentiate the reference and consider the term $\dot{r}(t)$ in the control input, the corresponds in the Laplace domain to a term $sR(s)$ where $R(s)$ is the Laplace transform of the reference. As a result, you will introduce a zero in the dynamics. This zero will not be there if you do not consider that derivative term.
However, it will be there if you consider an ODE of the form
$$\ddot{x}+\left[ 2-p \right]\dot{x}+\frac{\left[ 2-p \right]^{2}}{4}x=2\dot{u}+\frac{\left[ 2-p \right]^{2}}{4}u.$$
If you consider a step input, then you may have some overshoot which is due to the presence of the zero. I think that there is a confusion here between the two models.
For the second question, I am assuming that we are working with the transfer function with the zero, which can be rewritten as
$$G(s)=\dfrac{2s+\lambda^2}{(s+\lambda)^2}$$
where $\lambda=(p-2)/2$ and I am assuming that $p$ is such that $\lambda<0$.
It can be proven analytically for such a simple system that the zero may yield an overshoot but we can first get some intuition about that.
The zero is given here by $-\lambda^2/2$ and it is known that zeros introduce some amplification (check the Bode diagram). Therefore, if the zero acts on the lower frequency range before the poles attenuate the response, then we will have an overshoot. That is this will be the case if $\lambda^2/2<<\lambda$ or, in other words, if $\lambda<<2$. In fact, since the system is critically damped, this will happen if and only if $\lambda<2$.
Conversely, if the zero kicks in after the poles have started attenuating the response, then the effect of the zero will be negligible and there will be no overshoot. This will be the case if $\lambda>2$ and, similarly, since the system is critically damped, this will happen if and only if $\lambda\ge2$. Note that in the case where $\lambda=2$, we have a pole-zero cancellation and no overshoot since the resulting system is of first-order.
Another, more analytic, approach would be to notice that due the linearity of the system, the response of the system is the sum of the response to the transfer function
$$G_d(s)=\dfrac{2s}{(s+\lambda)^2},$$
and that of
$$G_p(s)=\dfrac{\lambda^2}{(s+\lambda)^2}.$$
We now that the step response of the second one is critically damped. Therefore, if there is an overshoot, it is necessarily due to the presence of the zero.
We have that
$$y_d(t)=2te^{-\lambda t}$$
and
$$y_p(t)=1 - \lambda te^{-\lambda t} - e^{-\lambda t}.$$
Note that $y_d$ increases first and decreases then to zero. In this regard, the only term that can cause the overshoot is the term that depends on $te^{-\lambda t}$ which is given for $y_d+y_p$ by
$$(2-\lambda)te^{-\lambda t}.$$
Since the system is critically damped, then we need this term to be nonpositive, otherwise it will create an overshoot. So, we need that $2-\lambda\le0$.
Best Answer
OK, there's a lot to unpack here, so let's take this one step at a time.
ODE solution
The ODE $$ Q\ddot{x}+2\zeta\omega_n\dot{x}+\omega_n^2x=\omega_n^2r $$ has the characteristic equation $$ Q\lambda^2 + 2\zeta\omega_n\lambda + \omega_n^2 = 0 $$ and eigenvalues $$ \lambda \in \left\{{\omega_n \over Q}\left(-\zeta + \sqrt{\zeta^2 - Q}\right), {\omega_n \over Q}\left(-\zeta -\sqrt{\zeta^2 - Q} \right) \right\} $$ so it can only be critically damped if $\zeta^2 = Q$. For the sake of argument, let's assume it is so that $$ \lambda = -{\omega_n\zeta \over Q} = -{\omega_n \over \zeta} $$ and $x$ has the homogeneous solution $$ x_h = (k_1 + k_2t)e^{-t\omega_n/\zeta}. $$ For the particular solution, since the input $\omega_n^2 r$ is constant, so is the steady-state response, therefore $x_p=C$. The complete solution is $$ x = x_p + x_h = C + (k_1 + k_2t)e^{-t\omega_n/\zeta} $$ with first derivative $$ \dot{x} = \left[k_2 - {\omega_n \over \zeta}(k_1 + k_2t) \right]e^{-t\omega_n/\zeta} $$ In control problems, it's usual to assume that the system has no energy stored at $t=0$ unless you have information indicating otherwise, so $x(0)=\dot{x}(0)=0$. From the initial conditions, we have $k_1 = -C$ and $$ k_2 = {\omega_n \over \zeta}k_1 = -{\omega_n \over \zeta}C. $$ To find $C$, we notice that $x = C$ and $\dot{x}=\ddot{x}=0$ for $t\to\infty$, so, if we were to plug our candidate solution and its derivatives back into the ODE and evaluate it for $t\to\infty$ we'd get $\omega_n^2C = \omega_n^2r$, hence $$ C = r\quad k_1 = -r\quad k_2 = -{\omega_n \over \zeta}r $$ which gives the final solution $$ x = r \left[ 1 - \left(1 + {\omega_n \over \zeta}t\right)e^{-t\omega_n/\zeta}\right]. $$
Having said all that, I believe you may have gotten your coefficients mixed up. If you want $\zeta$ to represent the damping ratio, your ODE has to be expressed in the form $\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2x = f(x)$.
System control
Now, for the control problem. You said you're trying to control a 1st-order unstable plant with a PID controller. Your controlled system then has a closed-loop transfer function of the form (assuming unity feedback) $$ H(s) = {(s^2K_d + sK_p + K_i)/p \over s^2(1 + K_d/p) + s[(K_p-l)/p] + K_i/p}, $$ where $$ \omega_n^2 = {K_i/p \over 1 +K_d/p} = {K_i \over p + K_d} $$ and $$ 2\zeta\omega_n = {(K_p - l)/p \over 1 + K_d/p} = {K_p - l \over p + K_d}. $$ If all you care about is having a critically damped system, you can find, after a few manipulations, the constraint $$ 2K_i^2 = (p+K_d)(K_p-l). $$ If you also need a specific natural frequency $\omega_n$, it can provide you an extra constraint that helps narrow down the possible values for $K_p$, $K_i$ and $K_d$. Of course, you have to make sure the closed-loop poles are stable.