Where does the $t$ come from in the solution to a critically damped differential equation? When solving homogeneous second order constant coefficient linear ODEs ($ay''+by'+cy=0$), there are three 'cases' that solutions fall into, based on the roots $r_1$, $r_2$ of the characteristic equation $ar^2+br+c=0$:


*

*$r_1, r_2 \in \Bbb R,\; r_1 \ne r_2$, aka "overdamped", for which the general solution is of the form $y(t) = c_1 e^{r_1 t} + c_2 e^{r_2 t}$.

*$r_1, r_2 \in \Bbb C,\; r_1=\overline{r_2}$, aka "underdamped", for which the general solution is $y(t) = c_1 e^{\Re(r)\,t}
   \cos\Im(r)\,t + c_2 e^{\Re(r)\,t} \sin\Im(r)\,t$.

*$r_1, r_2 \in \Bbb R,\; r_1 = r_2$, aka "critically damped", for which the general solution is $y(t) = c_1 e^{rt} + c_2 t e^{rt}$.


I understand the underdamped solution as essentially the same thing as the overdamped solution, equivalent to picking complex $c_1, c_2$ such that $y(t)$ is always real. But the critically damped case doesn't seem to fit neatly as a generalization of the other cases - that extra 't' makes it at least appear like it's something else entirely.
Where does the "$t$" actually come from in the critically damped case?
Why don't the other cases require or allow this?
 A: For a critically damped system, you have a double root so that the usual exponentials $e^{r_1t}$ and $e^{r_2t}$ merge in $e^{rt}$ and no longer form two linearly independent functions. Hence they cannot generate all solutions.
Now if you plug $te^{rt}$ in the equation, you get
$$a(2r+r^2t)e^{rt}+b(1+rt)e^{rt}+cte^{rt}=2ar+b+(ar^2+br+c)t.$$ The $t$ term obviously cancels out (whatever the damping), and so does the constant one because the double root verifies $2ar+b=0$.

A heuristic reasoning can justify the choice of such a term.
Let the damping converge to the critical one, with roots $r$ and $r+h$ getting closer and closer. We consider two independent solutions $e^{rt}$ and $e^{(r+h)t}-e^{rt}$, and amplify the second to avoid cancellation:
$$\lim_{h\to0}\frac{e^{(r+h)t}-e^{rt}}h=e^{rt}\lim_{h\to0}\frac{e^{ht}-1}h=e^{rt}t.$$
A: The general theory guarantees a two-dimensional solution space, whatever the given coefficients $a\ne0$, $b$, $c$. If the characteristic equation has a double root we obtain just one solution of the form $t\mapsto e^{\lambda t}$. There have to be other solutions, and it so happens that the $t$ thrown in works.
But you want to know why $t$, and not $\log t$, or somethig even more complicated. In this regard consider the following example:
$$y'' -\lambda y=0\ .$$
If $\lambda\ne0$ we have two different characteristic values, namely $\lambda$ and $0$, and the general solution is
$$y(t)=c_1 e^{\lambda t}+c_2\ .$$
The solution space ${\cal L}$ contains in particular the function
$$y_{\rm special}(t):={e^{\lambda t}-1\over\lambda}\ .$$
Now, if the underlying dynamical system is continuously modified (by turning a knob) such that $\lambda\to0$ then this special solution will converge to $y_*(t)=t$. It is then plausible that this $y_*(\cdot)$ will in fact be a solution in the limiting case.
A: We lose nothing by starting with the equation
$$ y'' + 2by' + cy = 0, $$
since the coefficient of $y''$ must be nonzero. We can proceed in a couple of different ways:
Reduction of Order
The simplest way is to suppose that $\alpha$ is a root of the characteristic equation $k^2+2bk+c=0$, and substitute $ y=ue^{\alpha x} $. This gives
$$ 0 = y'' + 2by' + cy = e^{\alpha x} \left( (\alpha^2+2b\alpha+c)u + u''+ 2(\alpha+b) u' \right). $$
The first term in the bracket vanishes since $\alpha$ is a root, so we end up having to solve
$$ u'' + 2(\alpha+b)u' = 0, $$
a first-order equation. Now, suppose that $\beta$ is the other root of the characteristic equation. We have
$$ k^2+2bk+c = (k-\alpha)(k-\beta) = k^2 - (\alpha+\beta)k+\alpha\beta, $$
so $2b = -\alpha-\beta$. Hence the new differential equation can be written as
$$ u'' + (\alpha-\beta)u' = 0. $$
There are now two possibilities:


*

*if $\alpha=\beta$, this is $u''=0$, which has solution $At+B$

*on the other hand, if $\alpha \neq \beta$, this equation can be solved by using the integrating factor $e^{(\alpha-\beta)x}$, so
$$ (e^{(\alpha-\beta)x}u')' = 0 \implies u' = Ae^{(\beta-\alpha)x}, $$
and hence
$$ u = A\frac{e^{(\beta-\alpha)x}}{\beta-\alpha} +B, $$
which obviously only works if $\alpha-\beta \neq 0$.


Limit of the Initial Value Problem
A differential equation normally comes with boundary conditions of some kind; the type that are easiest to use here for understanding are initial conditions, say $y(0) = q$, $y'(0) = p $. We can write down a general solution with unequal roots to this,
$$ y = Ae^{-bx}\cos{\omega t} + B e^{-bx} \frac{\sin{\omega t}}{\omega}, $$
where $\omega = \sqrt{c-b^2}$ is the square root of the discriminant. (This solution works for both real and imaginary $\omega$, and is always real, because $(\sin{\omega t})/\omega$ is real in both cases.) Note that this works because the quadratic formula gives the roots as $-b \pm i\omega$.
Now, the solution with the right initial conditions can be found to be
$$ y = q e^{-bx}\cos{\omega t} + (p+bq) e^{-bx} \frac{\sin{\omega t}}{\omega}  $$
by solving the equations simultaneously, and this is a good reason for writing the solution in the form we did in the first place. Now, we can rewrite this in terms of the roots:
$$ b = -\frac{\alpha+\beta}{2}, i\omega = \frac{\alpha-\beta}{2}, $$
where we chose which is $\alpha$ and which $\beta$ so that the sign makes sense. (Since $\cos{\omega t}$ and $(\sin{\omega t})/\omega$ are both even, it actually makes no difference.)
The upshot of all this is that taking $\beta \to \alpha$ is the same as taking $\omega$ to $0$. We have the limits
$$ \cos{\omega t} \to 1, \quad \frac{\sin{\omega t}}{\omega} \to t $$
as $\omega \to 0$, so the solution assumes the form we wanted.
The moral of the latter story is that $(\sin{\omega t})/\omega$ is the important function to consider in initial value problems, and has a limit that gives the answer we expect. (This function is also important in Fourier theory, and geometry of surfaces with constant curvature, as it happens.)
A: By dividing through by the coefficient $a$, we may as well assume our o.d.e. has the form
$$y'' + \beta y' + \gamma y = 0 .$$
As you say, critical damping means that the characteristic polynomial,
$$r^2 + \beta r + \gamma,$$
has a double root, say, $r_0$, so that
$r^2 + \beta r + \gamma = (r - r_0)^2$; in particular, $\beta = - 2 r_0$ and $\gamma = r_0^2$.
Now, expanding the l.h.s. of the differential equation $$(y e^{-r_0 t})'' = 0$$ by applying the Leibniz Rule twice yields
$$e^{-r_0 t} (y'' - 2 r_0 y' + r_0^2 y) = 0 ,$$
and multiplying through by the exponential factor and the above expressions for $\beta, \gamma$ shows that this differential equation is equivalent to our original one. Integrating twice both sides of $$(y e^{-r_0 t})'' = 0$$ gives that $y e^{-r_0 t}$ is some affine function $a t + b$, and multiplying through gives the general solution
$$y = e^{r_0 t}(a t + b)$$
as desired.
A: We have $a \ddot y + b \dot y + c y = 0$. Let $x_1 := y$ and $x_2 := \dot y$. Hence,
$$\begin{bmatrix} \dot x_1\\ \dot x_2\end{bmatrix} = \begin{bmatrix} 0 & 1\\ -\frac ca & -\frac ba\end{bmatrix} \begin{bmatrix} x_1\\ x_2\end{bmatrix}$$
More succinctly, 
$$\dot{\mathrm x} = \mathrm A \mathrm x$$
Let a Jordan decomposition of $\mathrm A$ be $\mathrm A = \mathrm P \mathrm J \mathrm P^{-1}$. Let $\eta := \mathrm P^{-1} x$. Hence, we have
$$\dot{\eta} = \mathrm J \eta$$
Integrating, we obtain
$$\eta (t) = \exp(t \mathrm J) \, \eta_0$$
where $\eta_0$ is a vector of initial conditions. In the critically damped case, $\mathrm A$ is non-diagonalizable and, thus, $\mathrm J$ has a single $2 \times 2$ Jordan block
$$\mathrm J = \begin{bmatrix} r & 1\\ 0 & r\end{bmatrix}$$
and, thus,
$$\exp(t \mathrm J) = \begin{bmatrix} e^{r t} & t \, e^{r t}\\ 0 & e^{r t}\end{bmatrix}$$
If $r < 0$ and assuming nonzero initial conditions, then cascading two leaky integrators with the same time constant leads to (non-oscillatory) resonance. Hence, the $t$ multiplying $e^{r t}$ in $t \, e^{r t}$.
