[proofplan]
We verify the Cayley transform by direct algebra. First we check that all displayed rational formulas have non-vanishing denominators on their stated domains, so the maps are holomorphic. Then we compute the defining function $\operatorname{Im}w-|z'|^2$ after applying $C$ and show that it is exactly a positive multiple of $1-|\zeta|^2$. We next prove directly that the proposed inverse formula sends $\mathcal U^n$ into $B^{n+1}$, and only after this domain check do we compose the two maps in both orders.
[/proofplan]
[step:Check that the Cayley transform is holomorphic on the unit ball]
Let
\begin{align*}
C:B^{n+1}\to\mathbb C^n\times\mathbb C,\qquad C(\zeta',\zeta_{n+1})=\left(\frac{\zeta'}{1+\zeta_{n+1}},i\frac{1-\zeta_{n+1}}{1+\zeta_{n+1}}\right)
\end{align*}
be the displayed rational map. If $\zeta=(\zeta',\zeta_{n+1})\in B^{n+1}$, then $|\zeta_{n+1}|\le |\zeta|<1$. Therefore $\zeta_{n+1}\ne -1$, so $1+\zeta_{n+1}\ne0$. Each component of $C$ is a quotient of holomorphic functions with non-vanishing denominator on $B^{n+1}$, hence $C$ is holomorphic on $B^{n+1}$.
[/step]
[step:Compute the Siegel defining function after applying $C$]
Fix $\zeta=(\zeta',\zeta_{n+1})\in B^{n+1}$, and define $(z',w):=C(\zeta',\zeta_{n+1})$. Put $a:=\zeta_{n+1}\in\mathbb C$. Then
\begin{align*}
z'=\frac{\zeta'}{1+a},\qquad w=i\frac{1-a}{1+a}.
\end{align*}
Since $\operatorname{Im}(iq)=\operatorname{Re}(q)$ for every $q\in\mathbb C$, we have
\begin{align*}
\operatorname{Im}w=\operatorname{Re}\left(\frac{1-a}{1+a}\right).
\end{align*}
Multiplying numerator and denominator by $1+\overline a$ gives
\begin{align*}
\frac{1-a}{1+a}=\frac{(1-a)(1+\overline a)}{|1+a|^2}.
\end{align*}
The real part of the numerator is $1-|a|^2$, because $a-\overline a$ is purely imaginary. Hence
\begin{align*}
\operatorname{Im}w=\frac{1-|a|^2}{|1+a|^2}.
\end{align*}
Also,
\begin{align*}
|z'|^2=\frac{|\zeta'|^2}{|1+a|^2}.
\end{align*}
Subtracting these two identities gives
\begin{align*}
\operatorname{Im}w-|z'|^2=\frac{1-|a|^2-|\zeta'|^2}{|1+a|^2}=\frac{1-|\zeta|^2}{|1+\zeta_{n+1}|^2}.
\end{align*}
Since $\zeta\in B^{n+1}$, the numerator $1-|\zeta|^2$ is positive and the denominator is positive. Thus $\operatorname{Im}w>|z'|^2$, so $C(\zeta)\in\mathcal U^n$.
[guided]
The target domain is defined by one real inequality, so the useful quantity to compute is not $w$ alone but the defining expression $\operatorname{Im}w-|z'|^2$. Fix $\zeta=(\zeta',\zeta_{n+1})\in B^{n+1}$ and define $(z',w):=C(\zeta',\zeta_{n+1})$. To simplify notation, set $a:=\zeta_{n+1}$. Then the formula for $C$ says
\begin{align*}
z'=\frac{\zeta'}{1+a},\qquad w=i\frac{1-a}{1+a}.
\end{align*}
We first compute $\operatorname{Im}w$. For any complex number $q$, multiplication by $i$ rotates real and imaginary parts, so $\operatorname{Im}(iq)=\operatorname{Re}(q)$. Applying this to $q=(1-a)/(1+a)$ gives
\begin{align*}
\operatorname{Im}w=\operatorname{Re}\left(\frac{1-a}{1+a}\right).
\end{align*}
To extract the real part, rationalize the denominator:
\begin{align*}
\frac{1-a}{1+a}=\frac{(1-a)(1+\overline a)}{|1+a|^2}.
\end{align*}
Expanding the numerator gives
\begin{align*}
(1-a)(1+\overline a)=1+\overline a-a-|a|^2.
\end{align*}
The term $\overline a-a$ is purely imaginary, so its real part is $0$. Therefore
\begin{align*}
\operatorname{Im}w=\frac{1-|a|^2}{|1+a|^2}.
\end{align*}
Next, from $z'=\zeta'/(1+a)$, the Euclidean norm on $\mathbb C^n$ gives
\begin{align*}
|z'|^2=\frac{|\zeta'|^2}{|1+a|^2}.
\end{align*}
Subtracting the two displayed formulas yields
\begin{align*}
\operatorname{Im}w-|z'|^2=\frac{1-|a|^2-|\zeta'|^2}{|1+a|^2}.
\end{align*}
Since $a=\zeta_{n+1}$ and $|\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*}
The denominator is positive because $1+\zeta_{n+1}\ne0$, and the numerator is positive because $\zeta\in B^{n+1}$. Hence $\operatorname{Im}w>|z'|^2$, which is exactly the condition that $C(\zeta)\in\mathcal U^n$.
[/guided]
[/step]
[step:Check that the proposed inverse is holomorphic on the Siegel domain]
Define
\begin{align*}
D:\mathcal U^n\to\mathbb C^n\times\mathbb C,\qquad D(z',w)=\left(\frac{2iz'}{w+i},\frac{i-w}{w+i}\right).
\end{align*}
If $(z',w)\in\mathcal U^n$, then $\operatorname{Im}w>|z'|^2\ge0$. In particular $\operatorname{Im}w>0$, so $w\ne -i$. Hence $w+i\ne0$ on $\mathcal U^n$. Each component of $D$ is therefore a quotient of holomorphic functions with non-vanishing denominator on $\mathcal U^n$, so $D$ is holomorphic on $\mathcal U^n$.
[/step]
[step:Show that the proposed inverse maps the Siegel domain into the ball]
Fix $(z',w)\in\mathcal U^n$, and define $(\eta',b):=D(z',w)$. Thus
\begin{align*}
\eta'=\frac{2iz'}{w+i},\qquad b=\frac{i-w}{w+i}.
\end{align*}
Since $\operatorname{Im}w>0$, we have $w+i\ne0$. We compute
\begin{align*}
1-|\eta'|^2-|b|^2=1-\frac{4|z'|^2}{|w+i|^2}-\frac{|i-w|^2}{|w+i|^2}.
\end{align*}
Using $|w+i|^2-|i-w|^2=4\operatorname{Im}w$, this becomes
\begin{align*}
1-|\eta'|^2-|b|^2=\frac{4(\operatorname{Im}w-|z'|^2)}{|w+i|^2}.
\end{align*}
Because $(z',w)\in\mathcal U^n$, the numerator is positive and the denominator is positive. Hence $|\eta'|^2+|b|^2<1$, so $D(z',w)=(\eta',b)\in B^{n+1}$.
[guided]
Before composing $C$ with $D$, we must prove that the output of $D$ actually lies in the domain of $C$. Fix $(z',w)\in\mathcal U^n$, and define $(\eta',b):=D(z',w)$, so
\begin{align*}
\eta'=\frac{2iz'}{w+i},\qquad b=\frac{i-w}{w+i}.
\end{align*}
The condition $(z',w)\in\mathcal U^n$ gives $\operatorname{Im}w>|z'|^2\ge0$, hence $\operatorname{Im}w>0$ and $w+i\ne0$. Therefore the displayed formula is defined.
To show $(\eta',b)\in B^{n+1}$, we compute its squared norm. By the Euclidean norm on $\mathbb C^n\times\mathbb C$,
\begin{align*}
|\eta'|^2+|b|^2=\frac{4|z'|^2}{|w+i|^2}+\frac{|i-w|^2}{|w+i|^2}.
\end{align*}
Thus
\begin{align*}
1-|\eta'|^2-|b|^2=1-\frac{4|z'|^2}{|w+i|^2}-\frac{|i-w|^2}{|w+i|^2}.
\end{align*}
Putting the terms over the common positive denominator $|w+i|^2$ gives
\begin{align*}
1-|\eta'|^2-|b|^2=\frac{|w+i|^2-|i-w|^2-4|z'|^2}{|w+i|^2}.
\end{align*}
For a complex number $w=u+iv$, with $u=\operatorname{Re}w$ and $v=\operatorname{Im}w$, we have $|w+i|^2=u^2+(v+1)^2$ and $|i-w|^2=u^2+(1-v)^2$. Hence $|w+i|^2-|i-w|^2=4v=4\operatorname{Im}w$. Therefore
\begin{align*}
1-|\eta'|^2-|b|^2=\frac{4(\operatorname{Im}w-|z'|^2)}{|w+i|^2}.
\end{align*}
The numerator is positive because $\operatorname{Im}w>|z'|^2$, and the denominator is positive because $w+i\ne0$. Hence $1-|\eta'|^2-|b|^2>0$, which is exactly $|\eta'|^2+|b|^2<1$. Therefore $D(z',w)\in B^{n+1}$.
[/guided]
[/step]
[step:Verify that the two formulas are mutual inverses]
First fix $\zeta=(\zeta',a)\in B^{n+1}$, where $a:=\zeta_{n+1}$. Let $(z',w):=C(\zeta',a)$. Then
\begin{align*}
w+i=i\frac{1-a}{1+a}+i=\frac{2i}{1+a}.
\end{align*}
Therefore the first component of $D(C(\zeta))$ is
\begin{align*}
\frac{2iz'}{w+i}=\frac{2i\zeta'/(1+a)}{2i/(1+a)}=\zeta',
\end{align*}
and the second component is
\begin{align*}
\frac{i-w}{w+i}=\frac{i-i(1-a)/(1+a)}{2i/(1+a)}=a.
\end{align*}
Thus $D\circ C=\operatorname{id}_{B^{n+1}}$.
Conversely fix $(z',w)\in\mathcal U^n$ and set
\begin{align*}
(\eta',b):=D(z',w)=\left(\frac{2iz'}{w+i},\frac{i-w}{w+i}\right).
\end{align*}
By the previous step, $(\eta',b)\in B^{n+1}$, so $C(D(z',w))$ is a genuine composition of the stated maps. Then
\begin{align*}
1+b=1+\frac{i-w}{w+i}=\frac{2i}{w+i}.
\end{align*}
The first component of $C(D(z',w))$ is
\begin{align*}
\frac{\eta'}{1+b}=\frac{2iz'/(w+i)}{2i/(w+i)}=z',
\end{align*}
and the second component is
\begin{align*}
i\frac{1-b}{1+b}=i\frac{1-(i-w)/(w+i)}{2i/(w+i)}=w.
\end{align*}
Thus $C\circ D=\operatorname{id}_{\mathcal U^n}$.
[/step]
[step:Conclude that the Cayley transform is a biholomorphism onto the Siegel domain]
The previous steps show that $C:B^{n+1}\to\mathcal U^n$ and $D:\mathcal U^n\to B^{n+1}$ are holomorphic, and that $D\circ C=\operatorname{id}_{B^{n+1}}$ and $C\circ D=\operatorname{id}_{\mathcal U^n}$. Hence $D:\mathcal U^n\to B^{n+1}$ is the inverse map of $C$. Therefore $C$ is a biholomorphism from $B^{n+1}$ onto $\mathcal U^n$, with inverse
\begin{align*}
C^{-1}(z',w)=\left(\frac{2iz'}{w+i},\frac{i-w}{w+i}\right).
\end{align*}
[/step]
