Settling time is usually defined as when when the step response is within 2% of the steady state response. Even though the difference between the step response and the steady state response gets bigger for all times when you use a bigger gain $K$, the percentage difference remains the same. That is why $K$ does not affect settling time.
The poles of your system are,
$$
p = \omega_n \left(-\zeta \pm \sqrt{\zeta^2 - 1}\right).
$$
The transient response of such a system will look like,
$$
y_t(t) = C_1\, e^{\omega_n \left(-\zeta + \sqrt{\zeta^2 - 1}\right) t} + C_2\, e^{\omega_n \left(-\zeta - \sqrt{\zeta^2 - 1}\right) t},
$$
where $C_1$ and $C_2$ depend on the initial conditions.
For $\zeta<1$ the poles will be each others complex conjugate and the transient can also be written as,
$$
y_t(t) = e^{-\omega_n\, \zeta\, t} \left(C_3 \cos\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right) + C_4 \sin\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right)\right).
$$
For a step response the total response can be written as the sum of the steady state and the transient response ($y(t)=y_{ss}(t)+y_t(t)$), with the steady state equal to $K$ (which can be found by setting $s$ to zero). For a step response it is assumed that the system is initial at rest, so position and velocity are zero for all $t\leq0$. Applying these constraints yields $C_3=-K$ and $C_4=0$,
$$
y(t) = K \left(1 - e^{-\omega_n\, \zeta\, t} \cos\left(\omega_n\,\sqrt{1 - \zeta^2}\, t\right)\right).
$$
Using an upper bound of one for the cosine, then the settling time can be approximated by only considering the exponential,
$$
\left|\frac{y(t)-y_{ss}(t)}{y_{ss}(t)}\right| \approx e^{-\omega_n\, \zeta\, t}.
$$
It can be noted that if you work out the left hand side of this equation yourself you will see that all $K$'s cancel.
The settling time for this approximation can be found by solving when $e^{-\omega_n\, \zeta\, t}=0.02$. This has the solution,
$$
T_{settling} = \frac{-\ln(0.02)}{\zeta\,\omega_n},
$$
$-\ln(0.02)$ is indeed close to $4$.
Solving for the constants for two real poles (when $\zeta\geq 1$) you get,
$$
C_1 = K \frac{\sqrt{\zeta^2 - 1} + \zeta}{2 \sqrt{\zeta^2 - 1}},
$$
$$
C_2 = K \frac{\sqrt{\zeta^2 - 1} - \zeta}{2 \sqrt{\zeta^2 - 1}}.
$$
For $\zeta$ equal or slightly bigger than 1 you will have to solve the following equation for $t$ in order to get the settling time,
$$
\left|e^{-\omega_n\, \zeta\, t} \left(\cosh\left(\omega_n \sqrt{\zeta^2 - 1}\, t\right) + \frac{\zeta}{\sqrt{\zeta^2 - 1}} \sinh\left(\omega_n \sqrt{\zeta^2 - 1}\, t\right)\right)\right| = 0.02.
$$
When $\zeta\gg1$ then $\zeta\, (\zeta^2 - 1)^{-1/2} \approx 1$ and $\cosh(t)+\sinh(t)=e^t$ can be used. Solving the above equation then has the solution,
$$
T_{settling} = \frac{-\ln(0.02)}{\left(\zeta - \sqrt{\zeta^2 - 1}\right) \omega_n}.
$$
This can also be seen as taking the biggest (least negative) pole, instead of the real part of the complex conjugate poles. Because when $\zeta$ is large the difference between the two real poles is large and only the "slowest" pole will dominate the transient response for the settling time.
Best Answer
The peak overshoot $\text{PO}$ (in percent) of a step response for stable systems can be calculated as the difference between the supremum of the time response $x_{\sup}$ (similar to the maximal value) and the terminal value $x_{\infty}$ of the step response for $t \to \infty$ divided by the terminal value. Hence
$$\text{PO}=100 \%\cdot \dfrac{x_{\sup}-x_{\infty}}{x_{\infty}}.$$
Let us consider the $n^{\text{th}}$ order system with a Heaviside step function $H(t)$ as input and zero initial conditions.
$$a_nx^{(n)}(t)+\ldots+a_1x^{(1)}(t)+a_0x(t)=H(t) \qquad x(0)=x^{(1)}(0)=\ldots=x^{(n-1)}(0)=0$$
In order to determine the solution of the previous ODE we need to determine the homogeneous $x_{\text{h}}(t)$ and the particular solution $x_{\text{p}}(t)$. The general solution is then given as the linear superposition of the homogeneous and particular solution:
$$x(t)=x_{\text{h}}(t)+x_{\text{p}}(t).$$
While the particular solution is easy to determine as $x_{\text{p}}(t)=1/a_0$ (note that $a_0>0$ is necessary for stability by the Stodola criterion) it is much more difficult to determine the homogeneous solution. In general, it is only possible to determine the exact solution for systems with order $n<5$, because we need to determine the roots of the characteristic polynomial of degree $n$ when determining the solution of the system (see Euler's exponential ansatz). It is a well-known fact that in general, it is not possible to determine the roots (using only radicals) for polynomials with degree higher than $4$.
Hence, if we cannot determine the roots of a general polynomial of degree $n\geq 5$ than we can not expect a closed form solution for the $\text{PO}$ for systems of order $n\geq 5$, because we cannot determine the exact form of $x(t)$.
For systems with an order $n\leq 4$ you might be able to determine the $\text{PO}$ with a closed formula but I expect it to be quite messy (if possible). You will need to solve $\dot{x}(t)=0$ (Fermat criterion for an extremal value of function aka first derivative equal to zero).