Simply put, the second diagram cannot exist as shown. If an LTI system overshoots at all, there must be an oscillatory component to the decay of the error waveform. Remember, the behavior of the system must be qualitatively identical, regardless of the initial value of the error!
Critically damped is when the system is on the verge of oscillatory decay. There is only one point at which this happens.
On the other hand, minimum settling time to within a certain error tolerance does allow a small amount of underdamping, with oscillation. Perhaps that's the distinction that the book is trying to make.
You need to be careful when using settle time approximations that relate \$\omega_n\$ and \$\zeta\$ to a \$T_s\$.
First, \$T_s\$ is defined as the time where the signal remains within \$\pm x\$% of the final value. Typically the settle limits are \$\pm 2\$% or \$\pm 5\$%, and these limits impact all second order system settle time approximations. Yes, every 2nd order settle time equation is an approximation, and they all differ as a function of \$\zeta\$.
Next, $$\omega_n \approx \frac{5.8335}{T_s},$$ is only a numerical approximation to the normalized unit step response for a critically damped system (i.e., \$\zeta=1\$), which is $$\theta(t)= \left(1-(1+\omega_n t)e^{-\omega_n t}\right)u(t),$$ where \$u(t)\$ is a unit step (aka Heaviside step). The settle time approximation for this is eloquently derived in your "this answer" link using a \$\pm 2\$% window to numerically arrive at \$\omega_n \approx \frac{5.8335}{T_s}\$.
When I look at the graphs of the data you shared, it looks like the 2% settle time for the critically damped trace is closer to 35. This results in \$\omega_n\approx\frac{5.8335}{35}= 0.167\$ rad/sec. A result is closer to the actual \$\omega_n\$ in the article. Better graphs, and clarification on which trace is actually the critically damped one, will improve this estimate.
Sorry, but I do not have free access to your article, so cannot answer your question directly, but most likely the authors knew \$\omega_n\$ directly.
Best Answer
TL;DR: NO, you can't use the underdamped settling time formula to find out the settling time of an overdamped system. And you can't use it for a critically damped system either.
LONG FORM answer follows...
Critically damped case
For the critically damped case (\$\zeta=1\$), the step response is:
$$ v_{out}(t) = H_0 u(t) \lbrack 1 - (1+\omega_0 t) e^{-\omega_0 t} \rbrack $$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$ 0.02 = (1+\omega_0 T_s) e^{-\omega_0 T_s}\\ $$
Solving numerically for \$\omega_0 T_s\$ (by simply using Excel's solver) we obtain:
$$ T_s \approx \frac{5.8335}{\omega_0} $$
Overdamped case
For the overdamped case (\$\zeta>1\$), the step response is:
$$ v_{out}(t) = H_0 u(t) \left[ 1 - \frac{s_2}{s_2-s_1}e^{s_1 t} - \frac{s_1}{s_1-s_2}e^{s_2 t} \right] $$
where \$s_1, s_2\$ are the real roots of the transfer function denominator:
$$ s_1 = -\zeta \omega_0 + \omega_0 \sqrt{\zeta^2-1} \\ s_2 = -\zeta \omega_0 - \omega_0 \sqrt{\zeta^2-1} $$
For convenience we define:
$$ \begin{align} \Delta &= \frac{s_2-s_1}{2} = - \omega_0 \sqrt{\zeta^2-1} \\ \Sigma &= \frac{s_1+s_2}{2} = - \zeta \omega_0 \\ K &= \frac{\Sigma}{\Delta} = \frac{\zeta}{\sqrt{\zeta^2-1}} \end{align} $$
So that:
$$ \begin{align} s_1 &= \Sigma-\Delta \\ s_2 &= \Sigma+\Delta \end{align} $$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$ \begin{align} 0.02 &= \frac{s_2}{s_2-s_1} e^{s_1 T_s} + \frac{s_1}{s_1-s_2} e^{s_2 T_s} = \\ &= \frac{\Sigma + \Delta}{2 \Delta} e^{(\Sigma - \Delta) T_s} - \frac{\Sigma - \Delta}{2 \Delta} e^{(\Sigma + \Delta) T_s} = \\ &= \frac{e^{\Sigma T_s}}{\Delta} \left[ \frac{\Sigma+\Delta}{2} e^{-\Delta T_s} - \frac{\Sigma-\Delta}{2} e^{\Delta T_s} \right] = \\ &= \frac{e^{\Sigma T_s}}{\Delta} \left[ \frac{\Delta}{2} \left( e^{\Delta T_s} + e^{-\Delta T_s} \right) - \frac{\Sigma}{2} \left( e^{\Delta T_s} - e^{-\Delta T_s} \right) \right] = \\ &= \frac{e^{\Sigma T_s}}{\Delta} \left[ \Delta \cosh{(\Delta T_s)} - \Sigma \sinh{(\Delta T_s)} \right] = \\ &= e^{K \Delta T_s} \left[ \cosh{(\Delta T_s)} - K \sinh{(\Delta T_s)} \right] = \\ &= e^{-K |\Delta| T_s} \left[ \cosh{(-|\Delta| T_s)} - K \sinh{(-|\Delta| T_s)} \right] \end{align} $$
And finally:
$$ 0.02 = e^{-K |\Delta| T_s} \left[ \cosh{(|\Delta| T_s)} + K \sinh{(|\Delta| T_s)} \right] \\ $$
Now that we have rewritten the expression in term of \$ |\Delta| T_s\$ and \$K\$ (instead of in terms of \$s_1\$ and \$s_2\$), we can numerically solve for \$ |\Delta| T_s\$, (by simply using Excel's solver) for any arbitrary given \$\zeta>1\$.
Example 1: a moderately overdamped system with \$\zeta = 1.1\$. Thus \$K = \frac{1.1}{1.1^2-1} \approx 2.4\$, and then solving numerically:
$$ T_s \approx \frac{3.172}{|\Delta|} = \frac{3.172}{\omega_0 \sqrt{1.1^2-1}} \approx \frac{6.922}{\omega_0} $$
Example 2: a heavily overdamped system with \$\zeta = 5\$. Thus \$K = \frac{5}{\sqrt{24}} \approx 1.0206\$, and then solving numerically:
$$ T_s \approx \frac{190.21}{|\Delta|} = \frac{190.21}{\omega_0 \sqrt{24}} \approx \frac{38.827}{\omega_0} $$
There is also an approximation for heavily overdamped (\$\zeta \gg 1\$) systems based on the dominant pole:
$$ v_{out}(t) \approx H_0 u(t) \left[ 1 - e^{s_1 t} \right] $$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$ 0.02 \approx e^{s_1 T_s} $$
and:
$$ T_s \approx \frac{\ln(0.02)}{s_1} = \frac{-\ln(0.02)}{\omega_0 (\zeta-\sqrt{\zeta^2-1})} $$
We can compare this approximation with the exact results that we have derived before.
For \$\zeta = 5\$:
$$ T_s \approx \frac{38.725}{\omega_0} $$
An estimation error just about -0.25%. Quite good indeed.
For \$\zeta = 1.1\$:
$$ T_s \approx \frac{6.096}{\omega_0} $$
An estimation error of approx -12%. Not bad taking into account that \$\zeta = 1.1\$ is just marginally above the critically damped case!.
Bonus
We can write a generic settling time expression for \$\zeta>1\$ as follows
$$ T_s = \frac{\psi}{\omega_0} $$
where \$\psi\$ is a coefficient roughly proportional to the damping factor \$\zeta\$.
I've numerically calculated the value of \$\psi\$ for a range of \$1<\zeta<9\$ using the expression previously derived for settling within 2% of the final value,
$$ 0.02 = e^{-K |\Delta| T_s} \left[ \cosh{(|\Delta| T_s)} + K \sinh{(|\Delta| T_s)} \right] $$
Then I've calculated (for comparison purposes) 1) the dominant pole approximation, 2) a 3rd order polynomial regression on my numerically calculated dataset, and 3), 4) the relative error due to these two approximations.
Here is an Excel plot with the results: