[proofplan]
We use the explicit Cayley formula and its inverse. The only possible pole of the boundary formula is where $1+\zeta_{n+1}=0$, which on the unit sphere forces the omitted point. A direct boundary identity gives
\begin{align*}
\operatorname{Im}(w)-|z|^2=\frac{1-|\zeta|^2}{|1+\zeta_{n+1}|^2}
\end{align*}
and shows that the punctured sphere maps into $\partial\mathcal U^n$, while the inverse formula gives a smooth inverse on the whole Siegel boundary. Finally, because the Cayley transform and its inverse are holomorphic near the relevant boundary points, their differentials preserve the $(1,0)$ complex tangent bundles, so the boundary diffeomorphism is CR.
[/proofplan]
[step:Extend the Cayley formula smoothly away from the south pole]
Let $S:=\partial B^{n+1}$ and let
\begin{align*}
p_-:=(0,\dots,0,-1)\in\mathbb C^n\times\mathbb C.
\end{align*}
Define
\begin{align*}
C_{\partial}:S\setminus\{p_-\}&\longrightarrow\mathbb C^n\times\mathbb C
\end{align*}
\begin{align*}
(\zeta',\zeta_{n+1})&\longmapsto \left(\frac{\zeta'}{1+\zeta_{n+1}}, i\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}\right).
\end{align*}
If $\zeta\in S$ and $1+\zeta_{n+1}=0$, then $\zeta_{n+1}=-1$. Since $|\zeta|^2=|\zeta'|^2+|\zeta_{n+1}|^2=1$, this forces $|\zeta'|^2=0$, hence $\zeta'=0$ and $\zeta=p_-$. Thus the denominator is nonzero on $S\setminus\{p_-\}$, and each component of $C_{\partial}$ is a smooth quotient of smooth functions there. Therefore $C_{\partial}$ is smooth on the punctured sphere.
[/step]
[step:Use the boundary identity to locate the image in the Siegel boundary]
Write
\begin{align*}
(z,w):=C_{\partial}(\zeta',\zeta_{n+1})
\end{align*}
for $\zeta\in S\setminus\{p_-\}$. Then
\begin{align*}
z=\frac{\zeta'}{1+\zeta_{n+1}}
\end{align*}
and
\begin{align*}
w=i\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}.
\end{align*}
Since $\operatorname{Im}(ia)=\operatorname{Re}(a)$ for every $a\in\mathbb C$, we have
\begin{align*}
\operatorname{Im}(w)=\operatorname{Re}\left(\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}\right)=\frac{1-|\zeta_{n+1}|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Also
\begin{align*}
|z|^2=\frac{|\zeta'|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Subtracting gives
\begin{align*}
\operatorname{Im}(w)-|z|^2=\frac{1-|\zeta_{n+1}|^2-|\zeta'|^2}{|1+\zeta_{n+1}|^2}=\frac{1-|\zeta|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Because $\zeta\in S$, $|\zeta|=1$, so
\begin{align*}
\operatorname{Im}(w)=|z|^2.
\end{align*}
Hence $C_{\partial}(S\setminus\{p_-\})\subset\partial\mathcal U^n$.
[guided]
The point of this computation is to check the target boundary equation with the exact normalization of the Cayley transform. Let $\zeta=(\zeta',\zeta_{n+1})\in S\setminus\{p_-\}$, and define $(z,w):=C_{\partial}(\zeta',\zeta_{n+1})$. Thus
\begin{align*}
z=\frac{\zeta'}{1+\zeta_{n+1}}
\end{align*}
and
\begin{align*}
w=i\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}.
\end{align*}
The Siegel boundary is the set of points satisfying $\operatorname{Im}(w)=|z|^2$, so we compute these two quantities from the formula.
For the second coordinate, use the elementary identity $\operatorname{Im}(ia)=\operatorname{Re}(a)$ for $a\in\mathbb C$. This gives
\begin{align*}
\operatorname{Im}(w)=\operatorname{Re}\left(\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}\right).
\end{align*}
Multiplying numerator and denominator by $1+\overline{\zeta}_{n+1}$ yields
\begin{align*}
\operatorname{Re}\left(\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}\right)=\frac{1-|\zeta_{n+1}|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
For the first coordinate,
\begin{align*}
|z|^2=\left|\frac{\zeta'}{1+\zeta_{n+1}}\right|^2=\frac{|\zeta'|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Therefore
\begin{align*}
\operatorname{Im}(w)-|z|^2=\frac{1-|\zeta_{n+1}|^2-|\zeta'|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Since $|\zeta|^2=|\zeta'|^2+|\zeta_{n+1}|^2$, this becomes
\begin{align*}
\operatorname{Im}(w)-|z|^2=\frac{1-|\zeta|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Now $\zeta\in S=\partial B^{n+1}$ means $|\zeta|=1$, and so the numerator is zero. Hence $\operatorname{Im}(w)=|z|^2$, which is precisely the condition that $(z,w)\in\partial\mathcal U^n$.
[/guided]
[/step]
[step:Restrict the inverse Cayley formula to obtain a smooth inverse on the Siegel boundary]
By [citetheorem:9200], the inverse Cayley transform on $\mathcal U^n$ is
\begin{align*}
C^{-1}(z,w)=\left(\frac{2iz}{w+i},\frac{i-w}{w+i}\right).
\end{align*}
Let
\begin{align*}
\Theta:\mathbb C^n\times\mathbb R&\longrightarrow\partial\mathcal U^n
\end{align*}
\begin{align*}
(z,t)&\longmapsto (z,t+i|z|^2)
\end{align*}
be the boundary parametrisation. For $(z,t)\in\mathbb C^n\times\mathbb R$, the denominator in the inverse formula is
\begin{align*}
w+i=t+i(|z|^2+1).
\end{align*}
This complex number is never zero, because its imaginary part is $|z|^2+1>0$. First define
\begin{align*}
\Psi:\mathbb C^n\times\mathbb R&\longrightarrow \mathbb C^n\times\mathbb C
\end{align*}
\begin{align*}
(z,t)&\longmapsto \left(\frac{2iz}{t+i(|z|^2+1)},\frac{i-t-i|z|^2}{t+i(|z|^2+1)}\right).
\end{align*}
The non-vanishing denominator shows that $\Psi$ is smooth. To verify that $\Psi$ takes values in $S$, set
\begin{align*}
D:=t+i(|z|^2+1).
\end{align*}
Then
\begin{align*}
|D|^2=t^2+(|z|^2+1)^2.
\end{align*}
For $\Psi(z,t)=(\eta',\eta_{n+1})$, direct computation gives
\begin{align*}
|\eta'|^2=\frac{4|z|^2}{|D|^2}
\end{align*}
and
\begin{align*}
|\eta_{n+1}|^2=\frac{|-t+i(1-|z|^2)|^2}{|D|^2}=\frac{t^2+(1-|z|^2)^2}{|D|^2}.
\end{align*}
Adding these identities yields
\begin{align*}
|\eta'|^2+|\eta_{n+1}|^2=\frac{4|z|^2+t^2+(1-|z|^2)^2}{t^2+(|z|^2+1)^2}=1,
\end{align*}
so $\Psi(z,t)\in S$. Also $\eta_{n+1}\ne -1$, because $\eta_{n+1}=-1$ would imply $i-t-i|z|^2=-t-i(|z|^2+1)$, hence $i=-i$, impossible. Thus $\Psi$ takes values in $S\setminus\{p_-\}$.
The identity $C\circ C^{-1}=\operatorname{id}_{\mathcal U^n}$ and $C^{-1}\circ C=\operatorname{id}_{B^{n+1}}$ from [citetheorem:9200] holds by equality of rational holomorphic maps wherever the displayed formulas are defined. The same rational formulas are defined at the boundary points under consideration because $1+\zeta_{n+1}\ne 0$ on $S\setminus\{p_-\}$ and $w+i\ne 0$ on $\partial\mathcal U^n$. Restricting these identities to those boundary points gives
\begin{align*}
C_{\partial}\circ \Psi=\Theta
\end{align*}
as maps $\mathbb C^n\times\mathbb R\to\partial\mathcal U^n$, and
\begin{align*}
\Psi\circ \Theta^{-1}\circ C_{\partial}=\operatorname{id}_{S\setminus\{p_-\}}.
\end{align*}
Thus $C_{\partial}:S\setminus\{p_-\}\to\partial\mathcal U^n$ is bijective, and its inverse is the smooth map induced by $\Psi$ through the parametrisation $\Theta$.
[/step]
[step:Check that the boundary map preserves complex tangent spaces]
Let $M:=S\setminus\{p_-\}$ and $N:=\partial\mathcal U^n$. For $q\in M$, define
\begin{align*}
T^{1,0}_qM:=T^{1,0}_q\mathbb C^{n+1}\cap (T_qM\otimes_{\mathbb R}\mathbb C),
\end{align*}
and for $r\in N$, define
\begin{align*}
T^{1,0}_rN:=T^{1,0}_r\mathbb C^{n+1}\cap (T_rN\otimes_{\mathbb R}\mathbb C).
\end{align*}
These are the induced CR structures on the two embedded hypersurfaces. Fix $\zeta\in M$ and set $(z,w):=C_{\partial}(\zeta)$. Since $1+\zeta_{n+1}\ne 0$, the same rational formula defines a holomorphic map on an open neighbourhood of $\zeta$ in $\mathbb C^{n+1}$. Let
\begin{align*}
\rho_B:\mathbb C^{n+1}&\longrightarrow\mathbb R
\end{align*}
\begin{align*}
\eta&\longmapsto |\eta|^2-1
\end{align*}
be a defining function for $S$, and let
\begin{align*}
\rho_U:\mathbb C^n\times\mathbb C&\longrightarrow\mathbb R
\end{align*}
\begin{align*}
(\xi,\omega)&\longmapsto \operatorname{Im}(\omega)-|\xi|^2
\end{align*}
be a defining function for $N$.
The identity from the previous step gives $\rho_U\circ C=0$ on $M$. If $Z\in T^{1,0}_{\zeta}M$, then $Z$ annihilates every smooth function that vanishes on $M$ near $\zeta$. Therefore
\begin{align*}
0=Z(\rho_U\circ C)=d\rho_U{}_{(z,w)}(dC_{\zeta}Z).
\end{align*}
Because $C$ is holomorphic, $dC_{\zeta}$ is complex-linear and sends $(1,0)$ vectors to $(1,0)$ vectors. Hence
\begin{align*}
dC_{\zeta}(Z)\in T^{1,0}_{(z,w)}N.
\end{align*}
Let $r\in N$ and let $W\in T^{1,0}_rN$. The inverse formula is holomorphic on a neighbourhood of $r$ because $w+i\ne 0$ on $N$, and $C_{\partial}^{-1}(N)\subset M$. Therefore $\rho_B\circ C^{-1}=0$ on $N$ near $r$. Since $W$ is tangent to $N$ and of type $(1,0)$, we have
\begin{align*}
0=W(\rho_B\circ C^{-1})=d\rho_B{}_{C_{\partial}^{-1}(r)}(dC^{-1}_rW).
\end{align*}
Holomorphicity of $C^{-1}$ implies that $dC^{-1}_rW$ is again a $(1,0)$ vector. Hence
\begin{align*}
dC^{-1}_rW\in T^{1,0}_{C_{\partial}^{-1}(r)}M.
\end{align*}
Thus $C_{\partial}$ and its smooth inverse both preserve the CR bundles, so $C_{\partial}$ is a CR diffeomorphism.
[/step]
[step:Identify the punctured sphere with the Heisenberg boundary model]
The map $\Theta:\mathbb C^n\times\mathbb R\to\partial\mathcal U^n$ is a smooth bijection with smooth inverse
\begin{align*}
\Theta^{-1}:\partial\mathcal U^n&\longrightarrow\mathbb C^n\times\mathbb R
\end{align*}
\begin{align*}
(z,w)&\longmapsto (z,\operatorname{Re}(w)).
\end{align*}
Indeed, on $\partial\mathcal U^n$ one has $\operatorname{Im}(w)=|z|^2$, so $w=\operatorname{Re}(w)+i|z|^2$. The standard identification of $\mathbb H^n$ with $\mathbb C^n\times\mathbb R$ equips $\mathbb H^n$ with the boundary CR structure pulled back by $\Theta$. Since $C_{\partial}:S\setminus\{p_-\}\to\partial\mathcal U^n$ is a CR diffeomorphism, the composition
\begin{align*}
\Theta^{-1}\circ C_{\partial}:S\setminus\{p_-\}\longrightarrow\mathbb H^n
\end{align*}
is a CR diffeomorphism. This proves that the Cayley boundary map identifies the punctured sphere with $\mathbb H^n$ as a CR manifold.
[/step]
