Clarification of a procedure to compute the product of the exponential of two matrices

In trying to understand a method outlined here (page 3, subroutine 1). Consider $$R_3 = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} .$$ Let $$A$$ be a square matrix over $$\mathbb{C}$$. Define $$X_3(A) = R_3 \otimes A + R_3^{\dagger} \otimes A^{\dagger}$$, a.k.a $$X_3(A) = \begin{bmatrix} 0_n & 0_n & A \\ 0_n & 0_n & 0_n \\ A^{\dagger} & 0_n & 0_n \end{bmatrix}$$ where $$0_n$$ is the $$n$$-dimensional zero matrix. $$X_3(A)$$ is Hermitian, and the author refers to this as embedding $$A$$ in a Hermitian matrix.

Next the author makes the assumption that for two square matrices $$A_1, A_2$$ of the same dimension, we have access to unitary operators $$e^{iX_3(A_1) \tau},e^{iX_3(A_2)\tau}$$ (this is possible since $$X_3(A_i)$$ is Hermitian), for $$t$$ some "simulation time", and $$n$$ a positive integer designated as the number of applications. This is described in "input assumption #1" (with $$\tau$$ a time-parameter). In order to obtain an estimate of $$e^{iX_3(A_1 + A_2)t}$$ , the procedure is described as:

Procedure: Sequentially apply $$e^{iA_1 t/n}, e^{iA_2 t/n}$$ for a total of $$n$$ consecutive times , defining

$$u_{add}(t) = (e^{iX_3(A_1)t/n} e^{iX_3(A_2) t/n})^n$$ where the number of applications of the unitaries with $$\tau = t/n$$ is proportional to $$n = O(t^2 / \epsilon)$$, where $$\epsilon$$ is the error term.

From what I understand for two matrices $$A,B$$, $$e^A e^B = e^{A+B}$$ is generally true only if $$A,B$$ commute, so the above expression say for $$n=2$$ would be $$e^{iX_3(A_1)t/2} e^{iX_3(A_2)t/2} e^{iX_3(A_1)t/2} e^{iX_3(A_2)t/2}$$, what I'm imagining the procedure shows is :

$$e^{iX_3(A_1)t/2} e^{iX_3(A_2)t/2} e^{iX_3(A_1)t/2} e^{iX_3(A_2)t/2} = e^{iX_3(A_1)t/2} e^{iX_3(A_1)t/2}e^{iX_3(A_2)t/2} e^{iX_3(A_2)t/2} = e^{iX_3(A_1)t} e^{iX_3(A_2)t}=e^{iX_3(A_1+A_2)t}$$,

but this assumes that $$X_3(A_1)$$ and $$X_3(A_2)$$ commute? I'm also not sure what the author means by "apply" these operators.

So, what you should really be thinking about is taking terms in the opposite order: $$e^{iH_1t/2}e^{iH_2t/2}e^{iH_1t/2}e^{iH_2t/2}\approx(e^{i(H_1+H_2)t/2})(e^{i(H_1+H_2)t/2})\approx e^{i(H_1+H_2)t}.$$ To understand this, be more explicit. Take the first term and compare it to the claimed approximation $$e^{iH_1\delta t}e^{iH_2t\delta t}-e^{i(H_1+H_2)\delta t}$$ in the limit that $$\delta t$$ is small. You can perform a Taylor expansion. Thus, $$(1+iH_1\delta t+\frac12H_1^2\delta t^2+\ldots)(1+iH_2\delta t+\frac12H_1^2\delta t^2+\ldots)-(1+i(H_1+H_2)\delta t+\frac12(H_1+H_2)^2\delta t^2+\ldots).$$ The terms of order $$\delta t^0$$ and $$\delta t^1$$ all cancel, so this approximation is accurate up to $$O(\delta t^2)$$. Now, to get a full evolution of time $$t$$, I need $$t/\delta t$$ time steps. All those $$\delta t^2$$ sized errors could add up, of which there are $$t/\delta t$$, so the overall error is $$O(\delta t)$$ for the full sequence. Thus, make $$\delta t$$ small enough and you've got a fairly accurate simulation.