First we note that the right-hand side of (REF ) is well defined. Indeed, since \(|{\langle }\xi ,\eta {\rangle }|\le 1\) 
and \(-1\le z_1, z_2 \le 1\) , we see \(|\sqrt{s}\big ( \tfrac{1-z_1}{2} {\langle }\xi ,\eta {\rangle }+ \tfrac{1+z_1}{2} z_2\big )|\le \sqrt{s}\) , so that both variables in \({\mathbf {P}}_n( U_{{\gamma },{\beta }+\frac{d-1}{2}})\)  are elements of \({\mathbb {U}}\) .
We first need to specify the reproducing kernel \({\mathbf {P}}_n(U_{a,b})\)  on \({\mathbb {U}}\)  when one of its variables is
on the curved boundary of \({\mathbb {U}}\) . Using the well-known identity [1]}
\(\frac{P_m^{(b-\frac{1}{2},b-\frac{1}{2})}(1) P_m^{(b-\frac{1}{2},b-\frac{1}{2})}(\rho )}{h_m^{(b-\frac{1}{2},b-\frac{1}{2})}}= \frac{C_m^b (1) C_m^b(\rho )}{h_m^b} = Z_m^b (\rho ),\) 
it follows from (REF ) and (REF ) that
\( & {\mathbf {P}}_n \left(U_{a,b}; (\sqrt{x_2},x_2), (y_1,y_2)\right) = \sum _{m=0}^n \frac{P_{m,n}^{a,b}(\sqrt{x_2},x_2)P_{m,n}^{a,b}(y_1,y_2)}{h_{m,n}^{a,b}}\\& \qquad = \sum _{m=0}^n \frac{c_{b+m,a}}{c_{b,a}}\frac{P_{n-m}^{(b+m,a)}(1-2x_2)P_{n-m}^{(b+m,a)}(1-2y_2)}{h_{n-m}^{(b+m,a)}} x_2^{\frac{m}{2}}y_2^{\frac{m}{2}}Z_m^b \left(\frac{y_1}{\sqrt{y_2}}\right). \) 

In terms of the orthogonal basis (REF ) and using the addition formula (REF ) for the
spherical harmonics, we obtain, with \({\alpha }= {\beta }+\frac{d-1}{2}\) ,
\({\mathsf {P}}_n(\varpi _{{\beta },{\gamma }}; (x,t),(y,s)) = \sum _{m=0}^n\frac{P_{n-m}^{({\alpha }+m,{\gamma })}(1-2 t)P_{n-m}^{({\alpha }+m,{\gamma })}(1-2 s)}{{\mathsf {h}}_{m,n}^{{\beta },{\gamma }}}t^{\frac{m}{2}} s^{\frac{m}{2}} Z_m^{\frac{d-2}{2}}({\langle }\xi ,\eta {\rangle }),\) 

where \(\xi = \frac{x}{\sqrt{t}} \in {\mathbb {S}^{d-1}}\)  and \(\eta = \frac{y}{\sqrt{s}} \in {\mathbb {S}^{d-1}}\) . For \({\beta }= -\frac{1}{2}\) , the sum in
the right-hand side can be identified, using (REF ), with \({\mathbf {P}}_n(U_{{\gamma },\frac{d-2}{2}})\)  in
(REF ) by setting \(x_2 = t\) , \(y_2 = s\)  and \(y_1 = \sqrt{s} {\langle }\xi ,\eta {\rangle }\) , which proves (REF ).
For \({\beta }> - \frac{1}{2}\) , we need to increase the value of the index in the zonal harmonic by using the following identity,
proved recently in [2]},
\(Z_m^{\lambda }(t) = c_{{\lambda },\sigma -1} c_{\sigma } \int _{-1}^1 \int _{-1}^1 & Z_m^{{\lambda }+\sigma } \left( \tfrac{1-z_1}{2} t + \tfrac{1+z_1}{2} z_2\right)\\& \times (1-z_1)^{\lambda }(1+z_1)^{\sigma -1} (1-z_2^2)^{\sigma -\frac{1}{2}} {\mathrm {d}}z_1 {\mathrm {d}}z_2 \) 

with \({\lambda }= \frac{d-2}{2}\)  and \(\sigma = {\beta }+\frac{1}{2}\) . This shows that
\({\mathsf {P}}_n \big (\varpi _{{\beta },{\gamma }};& (x,t), (y,s)\big ) =c_{\frac{d-2}{2},{\beta }-\frac{1}{2}} c_{{\beta }+\frac{1}{2}} \int _{-1}^1 \int _{-1}^1 \sum _{m=0}^n\frac{P_{n-m}^{({\alpha }+m,{\gamma })}(1-2 t)P_{n-m}^{({\alpha }+m,{\gamma })}(1-2 s)}{{\mathsf {h}}_{m,n}^{{\beta },{\gamma }}}\\& \times t^{\frac{m}{2}} s^{\frac{m}{2}} Z_m^{{\alpha }}\big ( \tfrac{1-z_1}{2} {\langle }\xi ,\eta {\rangle }+ \tfrac{1+z_1}{2} z_2\big )(1-z_1)^\frac{d-2}{2} (1+z_1)^{{\beta }-\frac{1}{2}} (1-z_2^2)^{{\beta }} {\mathrm {d}}z_1 {\mathrm {d}}z_2.\) 

The sum in the right-hand side can be identified, using (REF ), with the reproducing kernel
\({\mathbf {P}}_n(U_{{\gamma }, {\beta }+\frac{d-1}{2}})\)  in (REF ), which gives (REF ).
