Answering my own question.

It is known that the Lindelöf sum of the Mercator series, whose terms are $-(1-z)^n/n$ here, converges for all $z≠-1$ with magnitude 1, as they lie strictly on the interior of the Mittag-Leffler star. Therefore this converges for all orthogonal matrices $Z$ (with no eigenvalue being exactly $-1$) by the spectral theorem. Setting $Z=exp(K)exp(L)$, one can simply Lindelöf sum the BCH-D series. 

This holds if $exp(K)exp(L)$ does not have any eigenvalues equal to $-1$. If it does, one can always use the (trivial) formula $ln(exp(K)exp(L))$ $=$ $ln(exp(K)exp(L+2πiaI))-2πiaI$, each $ln$ being in any suitable branch and choosing real $a$ such that $exp(K)exp(L)exp(2πia)$ does not have any eigenvalues equal to $-1$(such an $a$ trivially exists), and apply the fact that the BCH-D series of $ln(e^Xe^Y)$ is equal to $X+Y$ plus nested commutators. Since commutators are invariant w.r.t. translation re: each component by $2πia*I$, one returns to the Lindelöf sum of the original series, even if $exp(K)exp(L)$ has any eigenvalues equal to $-1$. Q.E.D.

P.S. By the same token, the Lindelöf sum of the BCH-D series re: two finite-dimensional matrix-valued functions, real analytic, is also real analytic as long as its desired exponential is always normal.