For the impulse response of a 2nd order system given by
$$ H(s) = \frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2}$$
I was surprised to see there exists a general formula found by Wolfram Alpha for its inverse Laplace transform:
$$ h(t) = \frac{\omega_n}{\sqrt{\zeta^2-1}}\exp(-t\omega_n\zeta )\sinh\Big(t\omega_n \sqrt{\zeta^2-1}\Big)$$
Although this formula is not usually found engineering textbooks, it does seem to cover all cases. For example, in the under-damped case $0 < \zeta < 1$, usual formula can be recovered from the above by the simple facts that $\sqrt{\zeta^2-1} = i \sqrt{1-\zeta^2}$ in this case and $\sinh(ix)=i\sin x$ in general, so by cancelling out the imaginary units one obtains the [engineering] textbook result for the under-damped case as:
$$ h_\text{underdamped}(t) = \frac{\omega_n}{\sqrt{1-\zeta^2}}\exp(-t\omega_n\zeta )\sin\Big(t\omega_n \sqrt{1-\zeta^2}\Big)$$
Now what my question is actually about is whether there exist similarly nice, general formula for the step response of a 2nd system, which means finding the inverse Laplace transform of $Y(s)=\frac{1}{s}H(s)$, i.e. of
$$ Y(s) = \frac{\omega_n^2}{s(s^2+2\zeta\omega_ns+\omega_n^2)}$$
Alas the brains of Wolfram Alpha did not have much insight on this, resulting in a rather ugly looking solution. It is known however that at least for the under-damped case the step response can be written in a fairly nice form:
$$ y_\text{underdamped}(t) = 1 - \frac{1}{\sqrt{1-\zeta^2}}\exp(-t\omega_n\zeta )\sin\Big(t\omega_n \sqrt{1-\zeta^2} + \arccos \zeta\Big)$$
So, finally, my question is: is there a nice general form for the inverse Laplace transform of $Y(s)$ given above? By nice I mean something that doesn't look too far from the under-damped case so that it can be derived from the general case solution with little computation, hopefully similar to the way $h(t)$ relates to $h_\text{underdamped}(t)$.
Noting that the step response is just the integral of the impulse response
$$a(t)=\int_0^th(\tau)d\tau\tag{1}$$
you can enter the impulse response you found in Wolfram Alpha and obtain the result
$$a(t)=1-\frac{e^{-\zeta\omega_n t}}{\sqrt{\zeta^2-1}}\left[\zeta\sinh\left(t\omega_n\sqrt{\zeta^2-1}\right)+\sqrt{\zeta^2-1}\cosh\left(t\omega_n\sqrt{\zeta^2-1}\right)\right]\tag{2}$$
or, equivalently,
$$a(t)=1-\frac{e^{-\zeta\omega_n t}}{\sqrt{1-\zeta^2}}\left[\zeta\sin\left(t\omega_n\sqrt{1-\zeta^2}\right)+\sqrt{1-\zeta^2}\cos\left(t\omega_n\sqrt{1-\zeta^2}\right)\right]\tag{3}$$
Note that just like your formula for $h(t)$, both $(2)$ and $(3)$ are valid for the underdamped case ($0\le\zeta<1$) as well as for the overdamped case ($\zeta > 1$), but not for the critically damped case ($\zeta=1$).
For the overdamped case, Eq. $(2)$ can be rewritten as
$$a(t)=1-\frac{e^{-\zeta\omega_n t}}{\sqrt{\zeta^2-1}}\sinh\left(t\omega_n\sqrt{\zeta^2-1}+\text{arcosh}(\zeta)\right)\tag{4}$$
(by using the substitution $\zeta=\cosh(y)$).
In a similar manner, Eq. $(3)$ can be rewritten for the underdamped case as
$$a(t)=1-\frac{e^{-\zeta\omega_n t}}{\sqrt{1-\zeta^2}}\sin\left(t\omega_n\sqrt{1-\zeta^2}+\arccos(\zeta)\right)\tag{5}$$
(by the substitution $\zeta=\cos(y)$), just like in your question.
Finally, in the critically damped case $\zeta=1$, the impulse response is
$$h(t)=\omega_n^2te^{-\omega_nt}\tag{6}$$
and the corresponding step response is
$$a(t)=1-e^{-\omega_nt}(\omega_nt+1)\tag{7}$$
Note that $(6)$ and $(7)$ can be obtained from the corresponding general expressions for $\zeta\neq 1$ as limits for $\zeta\rightarrow 1$.