Andy's answer is good, but it can be generalized. For the transfer function:
$$H(s)=\frac{2\zeta\omega_n s}{s^2+2\zeta\omega_n s+\omega_n^2}$$
the inverse Laplace for all of the three mentioned cases have the same exponential decaying term in them:
$$exp(-\zeta\omega_n t)$$
That will be the envelope of the whole response. The rest of the time response is either of the form \$\sin{t}+\cos{t}\$ (oscillating, for underdamped), or \$\sinh{t}+\cosh{t}\$ (for overdamped), or a fixed term as a function of time, involving \$\omega_n\$ and \$\zeta\$.
If you really want them in full, here they are:
$$h(t)=\exp(-\omega_n\zeta t)\left[2\omega_n\zeta\cosh(\sqrt{\zeta^2-1}\omega_nt)-\frac{2\omega_n\zeta^2\sinh(\sqrt{\zeta^2-1}\omega_nt)}{\sqrt{\zeta^2-1}}\right]$$
$$h(t)=\exp(-\omega_n\zeta t)\left[2\omega_n\zeta(1-\omega_n\zeta t)\right]$$
$$h(t)=\exp(-\omega_n\zeta t)\left[2\omega_n\zeta\cos(\sqrt{1-\zeta^2}\omega_nt)-\frac{2\omega_n\zeta^2\sin(\sqrt{1-\zeta^2}\omega_nt)}{\sqrt{1-\zeta^2}}\right]$$
If you need to find the time when the envelope (alone) reaches a certain value x, then it's as simple as:
$$\exp(-\omega_n\zeta t)=x => -\omega_n\zeta t=\log x => t=-\frac{\log x}{\omega_n\zeta}$$
Update: This is only valid for the underdamped case, as there is a clear case of an oscillating term -- cos()-sin() -- multiplied by an inverse exponential. For the other two cases, it's impossible to do it analitically, as the time, t, is in more than one term, and that term is not oscillating, but either linearly variable (critically damped), or exponential (overdamped), so simplification is not possible. Also see the last pargraph.
Don't forget that, in all three cases, there is an additional term that can be factored out: \$2\omega_n\zeta\$, which gives the amplitude of the initial decay.
Note the \$\omega_n\$ in there, which tells you that the higher the frequency the higher the amplitude(!). This is related to an impulse that has the area of 1 (e.g. for a 1ms impulse, the amplitude is 1000). If you want to consider only a unity amplitude impulse, of zero duration (unity amplitude Dirac), then the extra term can be omitted (discarded).
Caveat emptor: the shape of the decay will not correspond to the final result in the case of \$\zeta\$>1, due to the hyperbolic functions! If the sinh() term would have been just like the cosh() term, the whole expression would have been reduced to an exp(-t), but sinh() is multiplied by \$\frac{\zeta}{\sqrt{\zeta^2-1}}\$, which means that the whole cosh()-sinh() term will "explode" if it were all by itself. Being multiplied by the exp() in front limits the response, but the envelope is no longer the same (contrary to the underdamped case). Here is what I mean:
G1, G2, G3 together with C1, C2 form the bandpass filter having V(o) as the output. V(exp) is the exponential decay (together with the \$\omega_n\zeta t\$ term, see Rpar=2*z*wn in B1), V(sinh_cosh) is the cosh()-sinh() term, V(test) is the multiplication of the previous two, compared with V(o) in the lower-most plot (they are the same). The upper two plots show the individual responses of the envelope and of the hyperbolic term -- this is what I meant when I said it "explodes", it goes to \$-\infty\$. And here is why I said that the envelope will not match the response, only by itself, due to the hyperbolic term:
Neither the output (blue trace), or its absolute value (black) will match the envelope alone (red). Which means that the solution above will not work, and the whole expression of the impulse response has to be considered, which means that, sadly, you're out of luck, since the expression contains the time, t, in both the envelope and the hyperbolic term. Even if you were to equivalate this with its equivalents, cosh=[exp(-x)+exp(-x)]/2 (minus for sinh), you cannot apply a logarithm to each exp(), only as a whole, which means the time cannot be extracted analitically, so you're stuck to numerical methods.
