$\newcommand{\bm}[1]{\boldsymbol{#1}}$ $\newcommand{\img}{\operatorname{img}}$
The Scalar Setting
When looking for a solution $u:\mathbb{R}\to\mathbb{R}$ of the following linear homogeneous second order ODE with constant coefficients:
\begin{equation} au''(t) + bu'(t) + cu(t) = 0, \end{equation}
one typically sets $u(t)=\exp(tx)\gamma$, $\gamma\in\mathbb{R}$, and substitutes in the equation:
\begin{equation} \begin{aligned} au''(t) + bu'(t) + cu(t) &= 0 \\ a\biggl(\exp(tx)\gamma\biggr)'' + b\biggl(\exp(tx)\gamma\biggr)' + c\biggl(\exp(tx)\gamma\biggr) &= 0 \\ \biggl(ax^2 + bx + c\biggr)\exp(tx)\gamma &= 0. \end{aligned} \end{equation}
Setting $\gamma=1$ and multiplying by $\exp(-tx)$ both sides yields a polynomial equation:
\begin{equation} \begin{aligned} ax^2 + bx + c = 0. \end{aligned} \end{equation}
If $a=b=c=0$ then any $x$ is a solution. If $a=b=0\ne c$ then this equation has no solution, but the one it was derived from has a solution for $x = -\infty$ (the discrepancy is because I assumed that I could divide both sides by $\exp(tx)$ which doesn't hold for $\exp(-\infty)=0$). If $a=0$ then the solution is $x = -(b^{-1}c)$ where we should interpret $0/0=0$ (we already took care of the case $c/0$, when we discussed $a=b=0\ne c$ by setting $x=-\infty$). In this case we have a first order ODE and the solution is:
\begin{equation} u(t) = \exp\left(-t(b^{-1}c)\right) u(0). \end{equation}
Finally if $a\ne 0$ then the roots are
\begin{equation} x_{1,2} = \frac{-(a^{-1}b) \pm \sqrt{(a^{-1}b)^2-4c}}{2} \end{equation}
If $(a^{-1}b)^2-4c<0$ then denote the conjugate roots $x_{1,2} = \alpha \pm i\beta$, $\alpha,\beta : \mathbb{R}\to\mathbb{R}$. Now there are three cases:
\begin{equation} u(t) = \begin{cases} \exp(t\lambda_1) \gamma_1 + \exp(t\lambda_2)\gamma_2, &\text{for }\,\, (a^{-1}b)^2-4c>0, \\ \exp(t\lambda_1) (\gamma_1 + t\gamma_2), &\text{for } \,\, (a^{-1}b)^2-4c=0, \\ \exp(t\alpha)(\cos(t\beta)\gamma_1 + \sin(t\beta)\gamma_2), &\text{for } \,\, (a^{-1}b)^2-4c<0, \\ \end{cases} \end{equation}
where the constants are picked as follows for the different cases: \begin{alignat}{3} (a^{-1}b)^2-4c>0 &\implies&& \begin{cases} \gamma_1 = (x_2-x_1)^{-1}(x_2u(0)-u'), \\ \gamma_2 = (x_2-x_1)^{-1}(u'(0)-x_1u(0)). \end{cases} \\ (a^{-1}b)^2-4c=0 &\implies&& \begin{cases} \gamma_{1} = u(0), \\ \gamma_{2} = u'(0) - x_1u(0). \end{cases} \\ (a^{-1}b)^2-4c<0 &\implies&& \begin{cases} \gamma_{1} = u(0), \\ \gamma_{2} = \beta^{-1}(u'(0)-\alpha u(0)). \end{cases} \end{alignat}
In the special case we have $x_1 = -x_2$, (which happens for $b=0$) we may also write the solution as:
\begin{equation} \begin{aligned} \exp(tx_1)\kappa_1 + \exp(-tx_1)\kappa_2 &= \exp(-i(itx_1))\kappa_1 + \exp(i(itx_1))\kappa_2 \\ &= \cos(itx_1)\kappa_1 - \sin(itx_1)\kappa_1 + \cos(itx_1)\kappa_2 + \sin(itx_1)\kappa_2 \\ &= \cos(itx_1)\gamma_1 + \sin(itx_1)\gamma_2 \\ &\implies \gamma_1 = u(0), \, \gamma_2 = -ix_1^{-1}u'(0). \end{aligned} \end{equation}
The Matrix Setting
Now I am looking for a vector function in time $\bm{u}:\mathbb{R}\to\mathbb{R}^n$ that solves the following second order matrix ODE:
$$\bm{A}\bm{u}''(t) + \bm{B}\bm{u}'(t) + \bm{C}\bm{u}(t) = \bm{0}.$$
I mimic the ansatz from the scalar-valued case and set $\bm{u} = \bm{\exp}(t\bm{X})\bm{\gamma}$ where $\bm{\gamma}\in\mathbb{R}^n$. Substituting in the ODE and differentiating yields:
\begin{equation} \begin{aligned} \bm{A}\bm{u}''(t) + \bm{B}\bm{u}'(t) + \bm{C}\bm{u}(t) &= \bm{0} \\ \bm{A}(\bm{\exp}(t\bm{X})\bm{\gamma})'' + \bm{B}(\bm{\exp}(t\bm{X})\bm{\gamma})' + \bm{C}(\bm{\exp}(t\bm{X})\bm{\gamma}) &= \bm{0} \\ (\bm{A}\bm{X}^2+\bm{B}\bm{X}+\bm{C})\bm{\exp}(t\bm{X})\bm{\gamma} &= \bm{0} \\ \end{aligned} \end{equation}
Here I have used that $(\bm{D}(t)\bm{\gamma})' = \bm{D}'(t)\bm{\gamma}$ and also the fact that if $\bm{D}(t)$ commutes with $\bm{D}'(t)$ (and $t\bm{X}$ indeed commutes with $\bm{X}$) then I may compute the derivative of the matrix exponential as:
\begin{equation*} \begin{aligned} \bigl(\bm{\exp}(\bm{D}(t))\bigr)' &= \sum_{k=0}^{\infty} \frac{(\bm{D}^k(t))'}{k!} \\ &= \sum_{k=1}^{\infty} \frac{\bm{D}'(t)\bm{D}^{k-1}(t) + \bm{D}(t)\bm{D}'(t)\bm{D}^{k-2}(t) + \ldots + \bm{D}^{k-1}(t)\bm{D}'(t)}{k!} \\ &= \sum_{k=1}^{\infty} \frac{k\bm{D}'(t)\bm{D}^{k-1}(t)}{k!} = \bm{D}'(t)\sum_{k=1}^{\infty} \frac{\bm{D}^{k-1}}{(k-1)!} = \bm{D}'(t) \bm{\exp}(\bm{D}(t)). \end{aligned} \end{equation*}
Now since $\bm{\gamma}$ was arbitrary I can consecutively set $\bm{\gamma} = \bm{e}_1$, $\ldots$, $\bm{\gamma}=\bm{e}_n$, and all of those together are equivalent to the matrix equation:
\begin{equation} \begin{aligned} (\bm{A}\bm{X}^2+\bm{B}\bm{X} + \bm{C})\bm{\exp}(t\bm{X}) &= \bm{0} \end{aligned} \end{equation}
I can now multiply both sides by $\bm{\exp}(-t\bm{X})$ to get the following quadratic matrix equation:
\begin{equation} \bm{A}\bm{X}^2+\bm{B}\bm{X}+\bm{C} = \bm{0}. \end{equation}
The Matrix Setting: Shared Eigenvectors
Let $\bm{A}, \, \bm{B}, \, \bm{C}$ be non-defective matrices having the same eigenvectors:
\begin{equation} \bm{A} = \bm{P}\bm{\Lambda}_{\bm{A}}\bm{P}^{-1}, \,\, \bm{B} = \bm{P}\bm{\Lambda}_{\bm{B}}\bm{P}^{-1}, \,\, \bm{C} = \bm{P}\bm{\Lambda}_{\bm{C}}\bm{P}^{-1}. \end{equation}
Then substituting $\bm{u}(t) = \bm{P}\bm{w}(t)$ and multiplying the system by $\bm{P}^{-1}$ from the left results in:
\begin{equation} \begin{aligned} \bm{A}\bm{u}''(t) + \bm{B}\bm{u}'(t) + \bm{C}\bm{u}(t) &= 0, \\ \bm{\Lambda}_{\bm{A}}\bm{w}''(t) + \bm{\Lambda}_{\bm{B}}\bm{w}'(t) + \bm{\Lambda}_{\bm{C}}\bm{w}(t) &= \bm{0}. \end{aligned} \end{equation}
Since these equations are decoupled we can just solve each row-equation separately as in the scalar case:
\begin{equation} \lambda_{A,i} w_i''(t) + \lambda_{B,i} w_{i}'(t) + \lambda_{C,i} w_i(t) = 0. \end{equation}
The Matrix Setting: Special Case $\bm{B}=\bm{0}$
One special case I came up with was if $\bm{A}$ is invertible and $\bm{B}=\bm{0}$ then I get the following solution:
\begin{align} \bm{A}\bm{u}'' + \bm{C}\bm{u} &= \bm{0} \\ \bm{u}'' &= -\bm{A}^{-1}\bm{C}\bm{u} \\ \bm{u}(t) &= \bm{\exp}(it\sqrt{\bm{A}^{-1}\bm{C}})\bm{\kappa}_1 + \bm{\exp}(-it\sqrt{\bm{A}^{-1}\bm{C}})\bm{\kappa}_2 \\ \bm{u}(t) &= \bm{\cos}(t\sqrt{\bm{A}^{-1}\bm{C}})\bm{\gamma}_1 + \bm{\sin}(t\sqrt{\bm{A}^{-1}\bm{C}})\bm{\gamma}_2 \\ \bm{u}(t) &= \bm{\cos}(t\sqrt{\bm{A}^{-1}\bm{C}})\bm{u}(0) + \bm{\sin}(t\sqrt{\bm{A}^{-1}\bm{C}})\bigl(\sqrt{\bm{A}^{-1}\bm{C}}\bigr)^{D}\bm{u}'(0). \end{align}
Here $\bm{L}^D = \bm{P}\bm{J}^{D}\bm{P}^{-1}$ denotes the Drazin inverse of $\bm{L} = \bm{P}\bm{J}\bm{P}^{-1}$ (I think I need this inverse so that its product with the $\bm{\sin}$ would not have a larger kernel than just $\bm{\sin}$). First of all, is the above derivation correct? What would I do if $\bm{A}$ is not invertible? My guess was to take the singular value decomposition of $\bm{A}$:
\begin{equation} \bm{A} = \bm{U}\begin{bmatrix} \bm{\Sigma} & \bm{0} \\ \bm{0} & \bm{0} \end{bmatrix} \bm{V}^T = \begin{bmatrix}\bm{U}_{\img} & \bm{U}_{\ker} \end{bmatrix} \begin{bmatrix} \bm{\Sigma} & \bm{0} \\ \bm{0} & \bm{0} \end{bmatrix} \begin{bmatrix}\bm{V}_{\img}^T \\ \bm{V}_{\ker}^T \end{bmatrix}, \end{equation}
and then perform the substitution $\bm{u} = \bm{V}_{\img} \bm{w}_{\img} + \bm{V}_{\ker}\bm{w}_{\ker}$, and multiply the equation with $\bm{U}^T$:
\begin{align} \bm{A}\bm{u}''&=-\bm{C}\bm{u} \\ \bm{U}^T\bm{A}(\bm{V}_{\img}\bm{w}_{\img}'' + \bm{V}_{\ker}\bm{w}_{\ker}'') &= -\bm{U}^T\bm{C}\bm{V}\bm{w} \\ \begin{bmatrix} \bm{\Sigma}\bm{w}''_{\img} \\ \bm{0} \end{bmatrix} &= -\begin{bmatrix} \bm{U}_{\img}^T \\ \bm{U}_{\ker}^T\end{bmatrix}\bm{C}\bm{V}\bm{w} \end{align} which can be decomposed into two subproblems: \begin{align} \bm{w}''_{\img} &= -\bm{\Sigma}^{-1}\bm{U}^T_{\img}\bm{C}(\bm{V}_{\img}\bm{v}_{\img} + \bm{V}_{\ker}\bm{w}_{\ker}) \\ \bm{0} &= \bm{U}^T_{\ker}\bm{C}\bm{V}\bm{w} \end{align}
The first is a second order equation with a non-homogenity, and the second looks like some kind of constraint. I have no idea how to write the solution of this analytically though.
The Matrix Setting: General Case
Now consider the case when $\bm{B}\ne\bm{0}$ in
\begin{equation} \bm{A}\bm{u}''(t) + \bm{B}\bm{u}'(t) + \bm{C}\bm{u}(t) = \bm{0}. \end{equation}
Once again we take $\bm{u}(t) = \bm{\exp}(t\bm{X})\bm{\gamma}$ for $\bm{\gamma}\in\mathbb{R}^n$. Then substituting in the above ODE we get:
\begin{equation} (\bm{A}\bm{X}^2 + \bm{B}\bm{X} + \bm{C})\bm{\exp}(t\bm{X})\bm{\gamma} = \bm{0}. \end{equation}
But since $\bm{\gamma}$ was arbitrary, we can set $\bm{\gamma} = \bm{e}_i$ consecutively for $1\leq i \leq n$ and we get:
\begin{equation} \begin{aligned} (\bm{A}\bm{X}^2 + \bm{B}\bm{X}+\bm{C})\bm{\exp}(-t\bm{X}) &= \bm{0}, \\ \bm{A}\bm{X}^2 + \bm{B}\bm{X}+\bm{C} &= \bm{0}. \end{aligned} \end{equation}
If I want to mimic the scalar case, then I need to solve a general quadratic matrix equation. I tried to complete the square:
\begin{equation} \bm{0} =\bm{A}\bm{X}^2 + \bm{B}\bm{X} + \bm{C} = \bm{A}(\bm{X}-\bm{H})^2 + \bm{K} = \bm{A}\bm{X}^2 - \bm{A}(\bm{X}\bm{H}+\bm{H}\bm{X}) + \bm{A}\bm{H}^2 + \bm{K}. \end{equation}
For $\bm{A}^{-1}\bm{B}$ I have some kind of Sylvester equation $-\bm{A}^{-1}\bm{B}\bm{X} = \bm{X}\bm{H}+\bm{H}\bm{X}$, which I don't know how to deal with. On the other hand if I assume that $\bm{H}$ and $\bm{X}$ commute, then I get $\bm{B} = -2\bm{A}\bm{H}$ and thus $\bm{H} = -\frac{1}{2}\bm{A}^{-1}\bm{B}$ similar to the scalar case. I can use this to find $\bm{K}$ also: $$\bm{K} = \bm{C}-\bm{A}\bm{H}^2 = \frac{4\bm{C} - \bm{B}\bm{A}^{-1}\bm{B}}{4}.$$ Now the solution for $\bm{X}$ is given as: $$\bm{X}_{1,2} = \bm{H}\pm\sqrt{-\bm{A}^{-1}\bm{K}} = \frac{1}{2}\left(-\bm{A}^{-1}\bm{B}\pm\sqrt{(\bm{A}^{-1}\bm{B})^{2} - 4\bm{A}^{-1}\bm{C}}\right)$$
If I assume that the solution proceeds similar to the scalar case, then for $\bm{K} =\bm{0}$ I have: $$\bm{u}(t) = \bm{\exp}(t\bm{X})(\bm{\gamma}_1 + t\bm{\gamma}_2),$$ and otherwise I have $$\bm{u}(t) = \bm{\exp}(t\bm{X}_1)\bm{\gamma}_1 + \bm{\exp}(t\bm{X}_2)\bm{\gamma}_2.$$
I feel like this doesn't quite match the scalar case, as the condition $\bm{K}=\bm{0}$ is missing cases such as just some rows being equal to zero (imagine diagonal matrices, where some equations have a single root and others double roots).
Does the above look correct? A major constraint here is that $\bm{X}$ and $\bm{H}$ must commute (e.g. this happens for $\bm{B}=\bm{0}$ or $\bm{K}=\bm{0}$). Is there an obvious way to generalise this to the case where those do not commute? What can I do about the case where $\bm{A}$ is not invertible?
Order Reduction Approach
I am aware that in practice I could also reduce the order of the ODE. As an example if I have the following ODE:
$$\sum_{j=0}^{m}\bm{A}_j\bm{u}^{(j)}(t) = \bm{0}$$
I can substitute $\bm{u}_j = \bm{u}^{(j)}$ and set up the coupled system $\bm{u}'_j =\bm{u}_{j+1}$. Then I may rewrite this as: \begin{align} \bm{u}_0' &= \bm{u}_1 \\ &\ldots \\ \bm{u}'_{m-2} &= \bm{u}_{m-1} \\ \bm{A}_m\bm{u}_{m-1}' &= -\sum_{j=0}^{m-1}\bm{A}_j\bm{u}_{j}. \end{align} If I assume that $\bm{A}_m$ is invertible then I can rewrite this as: \begin{equation} \bm{v}' = \bm{M}\bm{v} \iff \begin{bmatrix} \bm{u}_0' \\ \vdots \\ \bm{u}_{m-2}' \\ \bm{u}_{m-1}' \end{bmatrix} = \begin{bmatrix} \bm{0}& \bm{I} & \bm{0} & \bm{0} \\ \vdots& \ddots& \ddots & \bm{0}\\ \bm{0}&\ldots &\bm{0} & \bm{I} \\ -\bm{A}_m^{-1}\bm{A}_0 & \ldots & -\bm{A}_{m}^{-1}\bm{A}_{m-2} & -\bm{A}_{m}^{-1}\bm{A}_{m-1} \end{bmatrix} \begin{bmatrix} \bm{u}_0 \\ \bm{u}_1 \\ \vdots \\ \bm{u}_{m-1} \end{bmatrix} \end{equation}
Then the solution is given as $\bm{v}(t) = \bm{\exp}(t\bm{M})\bm{v}(0)$. This feels a bit unsatisfying though compared to the special structure I got e.g. for the solution of $\bm{A}\bm{u}'' + \bm{C}\bm{u} = \bm{0}$.