[proofplan]
We Taylor-expand both $e^{t(B+C)}$ and $e^{tB}e^{tC}$ to second order in $t$ using the matrix exponential series $e^X = I + X + \frac{1}{2}X^2 + O(\|X\|^3)$. Expanding $(B+C)^2 = B^2 + BC + CB + C^2$ in the first expression and collecting the product of two truncated series in the second, we find that the two expressions agree to first order and differ at second order by $\frac{t^2}{2}(BC - CB)$. Rearranging gives the stated formula. When $BC = CB$, the commutator vanishes, the $O(t^3)$ remainder can be shown to vanish to all orders by the Cauchy product of the two exponential series, and the exact factorisation $e^{t(B+C)} = e^{tB}e^{tC}$ follows.
[/proofplan]
[step:Expand $e^{t(B+C)}$ to second order]
The matrix exponential is defined by the convergent power series
\begin{align*}
e^{tX} = \sum_{n=0}^{\infty} \frac{t^n}{n!} X^n = I + tX + \frac{t^2}{2}X^2 + O(t^3).
\end{align*}
Applying this with $X = B + C$:
\begin{align*}
e^{t(B+C)} = I + t(B + C) + \frac{t^2}{2}(B + C)^2 + O(t^3).
\end{align*}
Expanding $(B + C)^2 = B^2 + BC + CB + C^2$:
\begin{align*}
e^{t(B+C)} = I + t(B + C) + \frac{t^2}{2}(B^2 + BC + CB + C^2) + O(t^3).
\end{align*}
[/step]
[step:Expand $e^{tB}e^{tC}$ to second order]
Expanding each factor to second order:
\begin{align*}
e^{tB} &= I + tB + \frac{t^2}{2}B^2 + O(t^3), \\
e^{tC} &= I + tC + \frac{t^2}{2}C^2 + O(t^3).
\end{align*}
Multiplying and collecting terms by powers of $t$. The zeroth-order term is $I \cdot I = I$. The first-order term is $I \cdot tC + tB \cdot I = t(B + C)$. The second-order term is $I \cdot \frac{t^2}{2}C^2 + tB \cdot tC + \frac{t^2}{2}B^2 \cdot I = \frac{t^2}{2}(B^2 + 2BC + C^2)$. Therefore:
\begin{align*}
e^{tB}e^{tC} = I + t(B + C) + \frac{t^2}{2}(B^2 + 2BC + C^2) + O(t^3).
\end{align*}
[guided]
When multiplying two power series in $t$, terms of a given order in $t$ arise from all products whose exponents sum to that order. For the $t^2$ term:
\begin{align*}
[e^{tB}e^{tC}]_{t^2} = \underbrace{I \cdot \frac{t^2}{2}C^2}_{\text{order } 0 \times 2} + \underbrace{tB \cdot tC}_{\text{order } 1 \times 1} + \underbrace{\frac{t^2}{2}B^2 \cdot I}_{\text{order } 2 \times 0} = \frac{t^2}{2}C^2 + t^2 BC + \frac{t^2}{2}B^2.
\end{align*}
Collecting: $\frac{t^2}{2}(B^2 + 2BC + C^2)$. Note the coefficient of $BC$ is $2$, not $1$ — this is because the cross term $tB \cdot tC = t^2 BC$ contributes a full $t^2 BC$, while there is no corresponding $t^2 CB$ term (the order of multiplication matters since $B$ and $C$ need not commute).
[/guided]
[/step]
[step:Subtract to obtain the splitting error]
Subtracting the expansion of $e^{t(B+C)}$ from that of $e^{tB}e^{tC}$:
\begin{align*}
e^{tB}e^{tC} - e^{t(B+C)} &= \frac{t^2}{2}(B^2 + 2BC + C^2) - \frac{t^2}{2}(B^2 + BC + CB + C^2) + O(t^3) \\
&= \frac{t^2}{2}(2BC - BC - CB) + O(t^3) \\
&= \frac{t^2}{2}(BC - CB) + O(t^3).
\end{align*}
Rearranging:
\begin{align*}
e^{t(B+C)} = e^{tB}e^{tC} - \frac{t^2}{2}(BC - CB) + O(t^3) = e^{tB}e^{tC} + \frac{t^2}{2}(CB - BC) + O(t^3),
\end{align*}
which is the stated formula.
[guided]
The second-order terms in $e^{t(B+C)}$ contain $BC + CB$ (both orderings appear in the square $(B+C)^2$), while $e^{tB}e^{tC}$ contains $2BC$ (only the ordering $BC$ appears, because $B$ is always applied after $C$ in the product $e^{tB}e^{tC}$). The difference is
\begin{align*}
2BC - (BC + CB) = BC - CB = [B, C],
\end{align*}
the commutator of $B$ and $C$. This is the fundamental reason why matrix exponentials do not factor in general: the non-commutativity of matrix multiplication introduces a splitting error proportional to the commutator $[B, C]$.
The $O(t^3)$ remainder comes from all higher-order terms in the Taylor expansions. A precise bound is $\|O(t^3)\| \leq C t^3 e^{t(\|B\| + \|C\|)}$ for a constant depending on the operator norms, but for the stated result the asymptotic form suffices.
[/guided]
[/step]
[step:Deduce exact factorisation when $BC = CB$]
When $BC = CB$, the commutator $CB - BC = 0$, so the splitting formula reduces to
\begin{align*}
e^{t(B+C)} = e^{tB}e^{tC} + O(t^3).
\end{align*}
To show that the equality is exact (not merely to second order), we use the Cauchy product of the exponential series. Since $BC = CB$, all powers of $B$ and $C$ commute, and the binomial theorem holds for matrices:
\begin{align*}
(B + C)^n = \sum_{j=0}^{n} \binom{n}{j} B^j C^{n-j}.
\end{align*}
Therefore:
\begin{align*}
e^{t(B+C)} = \sum_{n=0}^{\infty} \frac{t^n}{n!}(B+C)^n = \sum_{n=0}^{\infty} \frac{t^n}{n!} \sum_{j=0}^{n} \binom{n}{j} B^j C^{n-j} = \sum_{n=0}^{\infty} \sum_{j=0}^{n} \frac{t^j}{j!} B^j \cdot \frac{t^{n-j}}{(n-j)!} C^{n-j}.
\end{align*}
Re-indexing with $k = n - j$:
\begin{align*}
e^{t(B+C)} = \left(\sum_{j=0}^{\infty} \frac{t^j}{j!} B^j\right)\left(\sum_{k=0}^{\infty} \frac{t^k}{k!} C^k\right) = e^{tB} e^{tC}.
\end{align*}
Setting $t = 1$ gives $e^{B+C} = e^B e^C$.
[guided]
When $B$ and $C$ commute, the matrix algebra they generate is commutative, and all the usual identities from scalar exponentials carry over. In particular, the binomial theorem $(B+C)^n = \sum_{j=0}^n \binom{n}{j} B^j C^{n-j}$ holds — for non-commuting matrices, the expansion of $(B+C)^n$ involves all $2^n$ ordered products $X_1 X_2 \cdots X_n$ with each $X_i \in \{B, C\}$, grouped by the number of $B$'s and $C$'s but not by order. Commutativity collapses all products with the same numbers of $B$'s and $C$'s into a single term $B^j C^{n-j}$, weighted by $\binom{n}{j}$.
With the binomial expansion in hand, the Cauchy product formula for the double series gives
\begin{align*}
e^{t(B+C)} = \sum_{n=0}^{\infty} \sum_{j=0}^{n} \frac{(tB)^j}{j!} \cdot \frac{(tC)^{n-j}}{(n-j)!} = \left(\sum_{j=0}^{\infty} \frac{(tB)^j}{j!}\right)\left(\sum_{k=0}^{\infty} \frac{(tC)^k}{k!}\right) = e^{tB}e^{tC}.
\end{align*}
The convergence of both exponential series in operator norm justifies the rearrangement. This factorisation is exact — no remainder term — because commutativity eliminates the splitting error at every order, not just at order $t^2$.
[/guided]
[/step]
