[proofplan]
We compute the Lie brackets directly in the real coordinates $(x_1,\dots,x_n,y_1,\dots,y_n,t)$. The only non-constant coefficients in $X_j$ and $Y_k$ are the $t$-coefficients $2y_j$ and $-2x_k$, so the bracket $[X_j,Y_k]$ is entirely determined by differentiating those coefficients. Once the real commutators are known, the complex commutators follow from complex bilinearity of the Lie bracket after substituting $Z_j=(X_j-iY_j)/2$ and $\bar Z_j=(X_j+iY_j)/2$.
[/proofplan]
[step:Compute the bracket of $X_j$ and $Y_k$ from the coordinate coefficients]
Fix indices $j,k\in\{1,\dots,n\}$. Let
\begin{align*}
f:\mathbb H^n\to\mathbb C
\end{align*}
be a smooth function. Using the definition $[A,B]f=A(Bf)-B(Af)$ for smooth vector fields $A$ and $B$, we compute
\begin{align*}
X_j(Y_k f)=\left(\partial_{x_j}+2y_j\partial_t\right)\left(\partial_{y_k}f-2x_k\partial_t f\right).
\end{align*}
Expanding by the product rule gives
\begin{align*}
X_j(Y_k f)=\partial_{x_j}\partial_{y_k}f-2\delta_{jk}\partial_t f-2x_k\partial_{x_j}\partial_t f+2y_j\partial_t\partial_{y_k}f-4x_ky_j\partial_t^2 f.
\end{align*}
Similarly,
\begin{align*}
Y_k(X_j f)=\left(\partial_{y_k}-2x_k\partial_t\right)\left(\partial_{x_j}f+2y_j\partial_t f\right),
\end{align*}
and the product rule gives
\begin{align*}
Y_k(X_j f)=\partial_{y_k}\partial_{x_j}f+2\delta_{jk}\partial_t f+2y_j\partial_{y_k}\partial_t f-2x_k\partial_t\partial_{x_j}f-4x_ky_j\partial_t^2 f.
\end{align*}
Since $f$ is smooth, the mixed coordinate derivatives commute. Subtracting the second displayed identity from the first cancels all second-derivative terms and leaves
\begin{align*}
[X_j,Y_k]f=-4\delta_{jk}\partial_t f.
\end{align*}
Because this identity holds for every smooth function $f:\mathbb H^n\to\mathbb C$, we have
\begin{align*}
[X_j,Y_k]=-4\delta_{jk}T.
\end{align*}
[guided]
Fix $j,k\in\{1,\dots,n\}$. To prove an identity between vector fields, we test both sides on an arbitrary smooth function. Let
\begin{align*}
f:\mathbb H^n\to\mathbb C
\end{align*}
be smooth. The Lie bracket is defined by
\begin{align*}
[A,B]f=A(Bf)-B(Af)
\end{align*}
for smooth vector fields $A$ and $B$. Here
\begin{align*}
X_j=\partial_{x_j}+2y_j\partial_t
\end{align*}
and
\begin{align*}
Y_k=\partial_{y_k}-2x_k\partial_t.
\end{align*}
First compute $X_j(Y_k f)$:
\begin{align*}
X_j(Y_k f)=\left(\partial_{x_j}+2y_j\partial_t\right)\left(\partial_{y_k}f-2x_k\partial_t f\right).
\end{align*}
The only coefficient differentiated by $\partial_{x_j}$ is $-2x_k$, and
\begin{align*}
\partial_{x_j}(-2x_k)=-2\delta_{jk}.
\end{align*}
Therefore the product rule gives
\begin{align*}
X_j(Y_k f)=\partial_{x_j}\partial_{y_k}f-2\delta_{jk}\partial_t f-2x_k\partial_{x_j}\partial_t f+2y_j\partial_t\partial_{y_k}f-4x_ky_j\partial_t^2 f.
\end{align*}
Now compute $Y_k(X_j f)$:
\begin{align*}
Y_k(X_j f)=\left(\partial_{y_k}-2x_k\partial_t\right)\left(\partial_{x_j}f+2y_j\partial_t f\right).
\end{align*}
This time the only coefficient differentiated by $\partial_{y_k}$ is $2y_j$, and
\begin{align*}
\partial_{y_k}(2y_j)=2\delta_{jk}.
\end{align*}
Again using the product rule,
\begin{align*}
Y_k(X_j f)=\partial_{y_k}\partial_{x_j}f+2\delta_{jk}\partial_t f+2y_j\partial_{y_k}\partial_t f-2x_k\partial_t\partial_{x_j}f-4x_ky_j\partial_t^2 f.
\end{align*}
Now subtract. Since $f$ is smooth, the mixed coordinate derivatives commute, so the terms
$\partial_{x_j}\partial_{y_k}f$ and $\partial_{y_k}\partial_{x_j}f$ cancel, as do the corresponding mixed $t$-derivative terms and the $\partial_t^2 f$ terms. The remaining contribution is exactly the derivative of the two non-constant coefficients:
\begin{align*}
[X_j,Y_k]f=X_j(Y_k f)-Y_k(X_j f)=-4\delta_{jk}\partial_t f.
\end{align*}
Since $T=\partial_t$, this says
\begin{align*}
[X_j,Y_k]f=(-4\delta_{jk}T)f.
\end{align*}
The function $f$ was arbitrary, so the vector fields are equal:
\begin{align*}
[X_j,Y_k]=-4\delta_{jk}T.
\end{align*}
[/guided]
[/step]
[step:Show the remaining real coordinate brackets vanish]
Let $j,k\in\{1,\dots,n\}$ and let
\begin{align*}
f:\mathbb H^n\to\mathbb C
\end{align*}
be smooth. We compute each bracket by applying the vector fields to $f$.
First, for $X_j=\partial_{x_j}+2y_j\partial_t$ and $X_k=\partial_{x_k}+2y_k\partial_t$,
\begin{align*}
X_j(X_k f)=\partial_{x_j}\partial_{x_k}f+2y_k\partial_{x_j}\partial_t f+2y_j\partial_t\partial_{x_k}f+4y_jy_k\partial_t^2 f.
\end{align*}
Interchanging $j$ and $k$ gives
\begin{align*}
X_k(X_j f)=\partial_{x_k}\partial_{x_j}f+2y_j\partial_{x_k}\partial_t f+2y_k\partial_t\partial_{x_j}f+4y_jy_k\partial_t^2 f.
\end{align*}
Since $f$ is smooth, mixed coordinate derivatives commute, so $[X_j,X_k]f=0$.
Next, for $Y_j=\partial_{y_j}-2x_j\partial_t$ and $Y_k=\partial_{y_k}-2x_k\partial_t$,
\begin{align*}
Y_j(Y_k f)=\partial_{y_j}\partial_{y_k}f-2x_k\partial_{y_j}\partial_t f-2x_j\partial_t\partial_{y_k}f+4x_jx_k\partial_t^2 f.
\end{align*}
Interchanging $j$ and $k$ gives
\begin{align*}
Y_k(Y_j f)=\partial_{y_k}\partial_{y_j}f-2x_j\partial_{y_k}\partial_t f-2x_k\partial_t\partial_{y_j}f+4x_jx_k\partial_t^2 f.
\end{align*}
Again the mixed derivatives commute, hence $[Y_j,Y_k]f=0$.
Finally, since $T=\partial_t$,
\begin{align*}
X_j(Tf)=\partial_{x_j}\partial_t f+2y_j\partial_t^2 f.
\end{align*}
Also,
\begin{align*}
T(X_j f)=\partial_t\partial_{x_j}f+2(\partial_t y_j)\partial_t f+2y_j\partial_t^2 f=\partial_{x_j}\partial_t f+2y_j\partial_t^2 f,
\end{align*}
because $\partial_t y_j=0$. Thus $[X_j,T]f=0$. Similarly,
\begin{align*}
Y_j(Tf)=\partial_{y_j}\partial_t f-2x_j\partial_t^2 f
\end{align*}
and
\begin{align*}
T(Y_j f)=\partial_t\partial_{y_j}f-2(\partial_t x_j)\partial_t f-2x_j\partial_t^2 f=\partial_{y_j}\partial_t f-2x_j\partial_t^2 f,
\end{align*}
because $\partial_t x_j=0$. Hence $[Y_j,T]f=0$.
Since these identities hold for every smooth $f:\mathbb H^n\to\mathbb C$, we obtain
\begin{align*}
[X_j,X_k]=[Y_j,Y_k]=[X_j,T]=[Y_j,T]=0.
\end{align*}
[/step]
[step:Convert the real commutators into complex commutators]
Fix $j,k\in\{1,\dots,n\}$. By definition,
\begin{align*}
Z_j=\frac{1}{2}(X_j-iY_j)
\end{align*}
and
\begin{align*}
\bar Z_k=\frac{1}{2}(X_k+iY_k).
\end{align*}
Using complex bilinearity of the Lie bracket after extending the real vector fields complex-linearly, we get
\begin{align*}
[Z_j,\bar Z_k]=\frac{1}{4}[X_j-iY_j,X_k+iY_k].
\end{align*}
Expanding the bracket gives
\begin{align*}
[Z_j,\bar Z_k]=\frac{1}{4}\left([X_j,X_k]+i[X_j,Y_k]-i[Y_j,X_k]+[Y_j,Y_k]\right).
\end{align*}
The real commutator identities already proved give
\begin{align*}
[X_j,X_k]=[Y_j,Y_k]=0.
\end{align*}
Also,
\begin{align*}
[Y_j,X_k]=-[X_k,Y_j]=4\delta_{jk}T.
\end{align*}
Therefore
\begin{align*}
[Z_j,\bar Z_k]=\frac{1}{4}\left(i(-4\delta_{jk}T)-i(4\delta_{jk}T)\right)=-2i\delta_{jk}T.
\end{align*}
[/step]
[step:Verify that brackets among fields of the same complex type vanish]
For $j,k\in\{1,\dots,n\}$, complex bilinearity gives
\begin{align*}
[Z_j,Z_k]=\frac{1}{4}[X_j-iY_j,X_k-iY_k].
\end{align*}
Expanding and substituting the real commutator identities,
\begin{align*}
[Z_j,Z_k]=\frac{1}{4}\left([X_j,X_k]-i[X_j,Y_k]-i[Y_j,X_k]-[Y_j,Y_k]\right).
\end{align*}
Since $[X_j,X_k]=[Y_j,Y_k]=0$, $[X_j,Y_k]=-4\delta_{jk}T$, and $[Y_j,X_k]=4\delta_{jk}T$, the two middle terms cancel:
\begin{align*}
[Z_j,Z_k]=\frac{1}{4}\left(4i\delta_{jk}T-4i\delta_{jk}T\right)=0.
\end{align*}
Likewise,
\begin{align*}
[\bar Z_j,\bar Z_k]=\frac{1}{4}[X_j+iY_j,X_k+iY_k].
\end{align*}
Expanding and using the same real commutators gives
\begin{align*}
[\bar Z_j,\bar Z_k]=\frac{1}{4}\left([X_j,X_k]+i[X_j,Y_k]+i[Y_j,X_k]-[Y_j,Y_k]\right)=0.
\end{align*}
Thus all asserted real and complex Heisenberg commutation relations hold.
[/step]
