[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$.[/step]

[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: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]