As a first step and in response to OP's comments, in this post I aim to explain equation (1) in the arXiv version of Tao's paper and the resulting iterative scheme
To kind of circumvent the non-separability of a Hamiltonian H, Tao defines the still non-separable augmented Hamiltonian
\begin{equation*}
\tilde{H}
:= H_A + H_B + \omega H_C,
\end{equation*}
for some coupling constant $\omega \ge 0$, where
\begin{equation*}
H_A(q, p, x, y)
:= H(q, y), \quad
H_B(q, p, x, y)
:= H(x, p), \quad
H_C(q, p, x, y)
:= \frac{1}{2} \| q - x \|_2^2 + \frac{1}{2} \| p - y \|_2^2.
\end{equation*}
The key idea is now to approximate the flow of $\tilde{H}$ (which is essentially the flow of $H$) by composing the flows of $H_A$, $H_B$ and $H_C$.
With respect to the new variables $\tilde{q} := (q, x)$ and $\tilde{p} := (p, y)$, the Hamiltonian flow of $H_A$ is
\begin{equation*}
\begin{cases}
\dot{q} = \partial_p H(q, p, x, y) = 0, \\
\dot{p} = - \partial_q H(q, p, x, y) = - \partial_q H(q, y), \\
\dot{x} = \partial_y H(q, p, x, y) = \partial_y H(q, y), \\
\dot{y} = - \partial_x H(q, p, x, y) = 0,
\end{cases}
\end{equation*}
Note that this system is easily solvable in closed form:
1. Since $\dot{q} = \dot{y} = 0$, the variables $q$ and $y$ are constant in time, so $q(t) = q(0)$ and $y(t) = y(0)$.
2. For the other ODEs, the right hand side does not depend on the argument of the left hand side, they are essentially of the for $f'(t) = x$, so integrating yields $f(t) = t x + \text{const}$.
For some $\delta > 0$, the time $\delta$ flow map associated to $H_A$, i.e. the symplectomorphism with $\phi_{H_A}^{\delta}( q(0), p(0), x(0), y(0)) = (q(\delta), p(\delta), x(\delta), y(\delta))$, where $(( q(t), p(t), x(t), y(t)))_{t \ge 0}$ is the Hamiltonian flow of $H_A$ is thus given by
\begin{equation*}
\phi_{H_A}^{\delta}(q, p, x, y)
= \big(q, p - \delta \partial_q H(q, y), x + \delta \partial_y H(q, y), y\big).
\end{equation*}
By the same argument, the Hamiltonian flow of $H_B$ is
\begin{equation*}
\begin{cases}
\dot{q} = \nabla_p H_B = \nabla_p H(x, p), \\
\dot{p} = \nabla_y H_B = 0, \\
\dot{p} = - \nabla_q H_B = 0, \\
\dot{y} = - \nabla_x H_B = - \nabla_x H(x, p)
\end{cases}
\end{equation*}
and thus the flow map is
\begin{equation*}
\phi_{H_B}^{\delta}(q, p, x, y)
= \big(q + \delta \nabla_p H(x, p), p, x, y - \delta \nabla_x H(x, p)\big).
\end{equation*}
The Hamiltonian flow of $\omega H_C$ is
$$
\begin{cases}
\dot{q} = \omega (p - y), \\
\dot{x} = \omega (y - p), \\
\dot{p}= \omega (x - q), \\
\dot{y}= \omega (q - x)
\end{cases}
$$
This is a coupled system of linear first order ODEs of the form $\dot{w} = A w$, which has the solution $w(t) = \exp(t A) w_0$, where $w = (q, x, p, y)$. Calculating the matrix exponential gives
$$
\exp(t A)
= \frac{1}{2}\begin{pmatrix}
\cos(2 t) + 1 & \sin(2 t) & 1 - \cos(2 t) & -\sin(2 t) \\
-\sin(2 t) & \cos(2 t) + 1 & \sin(2 t) & 1 - \cos(2 t) \\
1 - \cos(2 t) & -\sin(2 t) & \cos(2 t) + 1 & \sin(2 t) \\
\sin(2 t) &1 - \cos(2 t) & -\sin(2 t) & \cos(2 t) + 1 \\
\end{pmatrix}
$$
and in this way you can derive
\begin{equation*}
\phi_{\omega H_C}^{\delta}(q, p, x, y)
= \frac{1}{2} \left(\begin{pmatrix} q + x \\ p + y \end{pmatrix} + R(\delta) \begin{pmatrix} q - x \\ p - y \end{pmatrix}, \begin{pmatrix} q + x \\ p + y \end{pmatrix} - R(\delta) \begin{pmatrix} q - x \\ p - y \end{pmatrix}\right),
\end{equation*}
where
\begin{equation*}
R(\delta) := \begin{pmatrix} \cos(2 \omega \delta) I & \sin(2 \omega \delta) I \\ - \sin(2 \omega \delta) I & \cos(2 \omega \delta) I \end{pmatrix}
\end{equation*}
is a rotation matrix.
One way to compose the flows of $H_A$, $H_B$ and $H_C$ is given by Strang splitting $$\phi^{\delta} := \phi_{H_A}^{\frac{\delta}{2}} \circ \phi_{H_B}^{\frac{\delta}{2}} \circ \phi_{\omega H_C}^{\delta} \circ \phi_{H_B}^{\frac{\delta}{2}} \circ \phi_{H_A}^{\frac{\delta}{2}}.$$
This means that our approximation of $q$ and $p$ at time $N \delta$ is given by iteratively applying $\phi^{\delta}$ to the initial conditions, that is, by (the first two components of) $$(\phi^{\delta})^N(q(0), p(0), x(0), y(0)).$$
One update step is hence of the form
$$
(q_{k + 1}, p_{k + 1}, x_{k + 1}, y_{k + 1})
= \phi^{\tau}(q_k, p_k, x_k, y_k),
$$
where $k$ are the discrete time steps and $\tau > 0$ is the step size.
Composing all the maps and writing these steps out, one could summarize this as follows:
\begin{align*}
\text{(First } H_A\text{ and } H_B \text{ steps) } & \begin{cases}
p_{k + \frac{1}{3}}
\gets p_k - \frac{\tau}{2} \delta_1 H(q_k, y_k), \\
x_{k + \frac{1}{3}}
\gets x_k + \frac{\tau}{2} \delta_2 H(q_k, y_k), \\
q_{k + \frac{1}{3}}
\gets q_k + \frac{\tau}{2} \delta_2 H(x_{k + \frac{1}{2}}, p_{k + \frac{1}{2}}), \\
y_{k + \frac{1}{3}}
\gets y_k - \frac{\tau}{2} \delta_2 H(x_{k + \frac{1}{2}}, p_{k + \frac{1}{2}}), \\
\end{cases} \\
\text{(} H_C \text{ step) }
& \begin{cases}
q_{k + \frac{2}{3}}
\gets \frac{1}{2} \left( q_{k + \frac{1}{3}} + x_{k + \frac{1}{3}} + \cos(2 \tau \omega) \left(q_{k + \frac{1}{3}} - x_{k + \frac{1}{3}} \right) + \sin(2 t) \left(y_{k + \frac{1}{3}} - p_{k + \frac{1}{3}}\right) \right) , \\
p_{k + \frac{2}{3}}
\gets \frac{1}{2} \left( p_{k + \frac{1}{3}} + y_{k + \frac{1}{3}} + \cos(2 \tau \omega) (p_{k + \frac{1}{3}} - y_{k
+ \frac{1}{3}}) - \sin(2 t) (q_{k + \frac{1}{3}} - x_{k + \frac{1}{3}}) \right), \\
x_{k + \frac{2}{3}}
\gets \frac{1}{2} \left( q_{k + \frac{1}{3}} + x_{k + \frac{1}{3}} - \cos(2 \tau \omega) \left(q_{k + \frac{1}{3}} - x_{k + \frac{1}{3}} \right) - \sin(2 t) \left(y_{k + \frac{1}{3}} - p_{k + \frac{1}{3}}\right) \right), \\
y_{k + \frac{2}{3}}
\gets \frac{1}{2} \left( p_{k + \frac{1}{3}} + y_{k + \frac{1}{3}} - \cos(2 \tau \omega) (p_{k + \frac{1}{3}} + y_{k
+ \frac{1}{3}}) - \sin(2 t) (q_{k + \frac{1}{3}} - x_{k + \frac{1}{3}}) \right), \\
\end{cases} \\
\text{(Second } H_A\text{ and } H_B \text{ steps) } & \begin{cases}
q_{k + 1}
\gets q_{k + \frac{2}{3}} + \frac{\tau}{2} \delta_2 H(x_{k + \frac{2}{3}}, p_{k + \frac{2}{3}}), \\
y_{k + 1}
\gets y_{k + \frac{2}{3}} - \frac{\tau}{2} \delta_1 H(x_{k + \frac{2}{3}}, p_{k + \frac{2}{3}}), \\
p_{k + 1}
\gets p_{k + \frac{2}{3}} - \frac{\tau}{2} \delta_1 H(q_{k + 1}, y_{k + 1}), \\
x_{k + 1}
\gets x_{k + \frac{2}{3}} + \frac{\tau}{2} \delta_2 H(q_{k + 1}, y_{k + 1}).
\end{cases}
\end{align*}
But I think one could also replace the $H_C$ step by implicit steps, which are still available in closed-form.
