[proofplan]
We prove that the displayed kernel is a positive normalized approximate identity on the boundary hyperplane and that each kernel slice is harmonic in the interior variable. Harmonicity follows by differentiating under the integral sign after checking uniform domination of the first two derivatives on compact subsets of the half-space. The boundary convergence is the approximate-identity argument: the kernel mass is one, continuity controls the near part, and an explicit tail estimate forces the far part to vanish as $x_n \downarrow 0$ and $x' \to y_0$.
[/proofplan]
[step:Verify that the Poisson kernel has total mass one]
Fix $x = (x',x_n) \in \mathbb{R}^n_+$. Define the translation-scaling map $T_x: \mathbb{R}^{n-1} \to \mathbb{R}^{n-1}$ by $T_x(y)=(y-x')/x_n$.
The [Change of Variables Theorem](/theorems/22), applied to the affine $C^1$ diffeomorphism $T_x^{-1}: \mathbb{R}^{n-1} \to \mathbb{R}^{n-1}$, gives $y = x' + x_n z$ and
\begin{align*}
d\mathcal{L}^{n-1}(y) = x_n^{n-1}\,d\mathcal{L}^{n-1}(z).
\end{align*}
Therefore
\begin{align*}
\int_{\mathbb{R}^{n-1}} P(x,y)\,d\mathcal{L}^{n-1}(y) = \frac{2}{n\omega_n}\int_{\mathbb{R}^{n-1}}\frac{1}{(1 + |z|^2)^{n/2}}\,d\mathcal{L}^{n-1}(z).
\end{align*}
It remains to compute the last integral. For each $m \in \mathbb{N}$, let $\omega_m := \mathcal{L}^m(B(0,1))$ denote the volume of the unit ball in $\mathbb{R}^m$. The polar-coordinate formula for [Lebesgue measure](/page/Lebesgue%20Measure) in $\mathbb{R}^{n-1}$ applies because the integrand is non-negative and radial, and gives
\begin{align*}
\int_{\mathbb{R}^{n-1}}\frac{1}{(1 + |z|^2)^{n/2}}\,d\mathcal{L}^{n-1}(z) = (n-1)\omega_{n-1}\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r).
\end{align*}
We compute this one-dimensional integral by the substitution $s=r^2/(1+r^2)$, which maps $(0,\infty)$ bijectively onto $(0,1)$. Its inverse is $r=(s/(1-s))^{1/2}$, and the one-dimensional change-of-variables formula gives
\begin{align*}
\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r)=\frac{1}{2}\int_0^1 s^{(n-3)/2}(1-s)^{-1/2}\,d\mathcal{L}^1(s)
\end{align*}
Since $n \geq 2$, the beta parameters $(n-1)/2$ and $1/2$ are positive. Let $\Gamma:(0,\infty)\to(0,\infty)$ denote the Gamma function. The beta-gamma identity therefore gives
\begin{align*}
\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r)=\frac{\sqrt{\pi}\,\Gamma((n-1)/2)}{2\Gamma(n/2)}.
\end{align*}
Together with the unit-ball volume formula $\omega_m=\pi^{m/2}/\Gamma(m/2+1)$, obtained from the [Gaussian integral](/theorems/1140) in $\mathbb{R}^m$, the right-hand side equals
\begin{align*}
\frac{n\omega_n}{2}.
\end{align*}
Hence
\begin{align*}
\int_{\mathbb{R}^{n-1}} P(x,y)\,d\mathcal{L}^{n-1}(y) = 1.
\end{align*}
[guided]
Fix an interior point $x = (x',x_n) \in \mathbb{R}^n_+$. The purpose of this step is to prove that $P(x,\cdot)$ is a probability density on the boundary hyperplane. This normalization is the reason the Poisson integral can reproduce constant boundary data and later allows us to write $u(x) - g(y_0)$ as an average of $g(y) - g(y_0)$.
Define $T_x: \mathbb{R}^{n-1} \to \mathbb{R}^{n-1}$ by $T_x(y)=(y-x')/x_n$. Since $x_n > 0$, this is a bijective affine $C^1$ map with $C^1$ inverse. By the [Change of Variables Theorem](/theorems/22), under the substitution $z = T_x(y)$, we have $y = x' + x_n z$ and the Jacobian determinant in dimension $n-1$ is $x_n^{n-1}$. Thus
\begin{align*}
d\mathcal{L}^{n-1}(y) = x_n^{n-1}\,d\mathcal{L}^{n-1}(z).
\end{align*}
Also,
\begin{align*}
|x' - y|^2 + x_n^2 = x_n^2|z|^2 + x_n^2 = x_n^2(1+|z|^2).
\end{align*}
Substituting these two identities into the integral of $P(x,\cdot)$ yields
\begin{align*}
\int_{\mathbb{R}^{n-1}} P(x,y)\,d\mathcal{L}^{n-1}(y) = \frac{2}{n\omega_n}\int_{\mathbb{R}^{n-1}}\frac{1}{(1 + |z|^2)^{n/2}}\,d\mathcal{L}^{n-1}(z).
\end{align*}
We now compute the remaining radial integral. For each $m \in \mathbb{N}$, let $\omega_m := \mathcal{L}^m(B(0,1))$ denote the volume of the unit ball in $\mathbb{R}^m$. The polar-coordinate formula for Lebesgue measure in $\mathbb{R}^{n-1}$ applies here because the integrand is non-negative and depends only on $|z|$. It gives
\begin{align*}
\int_{\mathbb{R}^{n-1}}\frac{1}{(1 + |z|^2)^{n/2}}\,d\mathcal{L}^{n-1}(z) = (n-1)\omega_{n-1}\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r).
\end{align*}
To evaluate the one-dimensional integral, set $s=r^2/(1+r^2)$. This substitution sends $(0,\infty)$ onto $(0,1)$, has inverse $r=(s/(1-s))^{1/2}$, and transforms the measure by the one-dimensional change-of-variables formula. Hence
\begin{align*}
\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r) = \frac{1}{2}\int_0^1 s^{(n-3)/2}(1-s)^{-1/2}\,d\mathcal{L}^1(s).
\end{align*}
Since $n \geq 2$, the beta parameters $(n-1)/2$ and $1/2$ are positive, so the beta-gamma identity applies. Let $\Gamma:(0,\infty)\to(0,\infty)$ denote the Gamma function. We obtain
\begin{align*}
\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r) = \frac{\sqrt{\pi}\,\Gamma((n-1)/2)}{2\Gamma(n/2)}.
\end{align*}
Using the unit-ball volume formula, $(n-1)\omega_{n-1} = 2\pi^{(n-1)/2}/\Gamma((n-1)/2)$ and $n\omega_n/2 = \pi^{n/2}/\Gamma(n/2)$, so the product is
\begin{align*}
(n-1)\omega_{n-1}\int_0^\infty \frac{r^{n-2}}{(1+r^2)^{n/2}}\,d\mathcal{L}^1(r) = \frac{n\omega_n}{2}.
\end{align*}
Therefore
\begin{align*}
\int_{\mathbb{R}^{n-1}} P(x,y)\,d\mathcal{L}^{n-1}(y) = 1.
\end{align*}
[/guided]
[/step]
[step:Show that each kernel slice is harmonic in the interior variable]
Fix $y \in \mathbb{R}^{n-1}$ and define
\begin{align*}
K_y: \mathbb{R}^n_+ \to \mathbb{R}, \quad x \mapsto P(x,y).
\end{align*}
Define the displacement map $X_y: \mathbb{R}^n_+ \to \mathbb{R}^n$ by $X_y(x)=(x' - y,x_n)$. Define the radius function $\rho_y: \mathbb{R}^n_+ \to (0,\infty)$ by $\rho_y(x)=|X_y(x)|$.
Then
\begin{align*}
K_y(x) = \frac{2}{n\omega_n}x_n\rho_y(x)^{-n}.
\end{align*}
Since $x_n > 0$, the vector $X_y(x)$ is never zero, so $K_y$ is smooth on $\mathbb{R}^n_+$. Write $X_{y,i}(x)$ for the $i$-th coordinate of $X_y(x)$. For each $i\in\{1,\dots,n\}$,
\begin{align*}
\partial_{x_i}(\rho_y^{-n}) = -nX_{y,i}\rho_y^{-n-2}
\end{align*}
and
\begin{align*}
\partial_{x_i}^2(\rho_y^{-n}) = -n\rho_y^{-n-2}+n(n+2)X_{y,i}^2\rho_y^{-n-4}.
\end{align*}
Summing over $i$ and using $\sum_{i=1}^n X_{y,i}^2=\rho_y^2$ gives
\begin{align*}
\Delta(\rho_y^{-n}) = 2n\rho_y^{-n-2}.
\end{align*}
Also, since $X_{y,n}(x)=x_n$,
\begin{align*}
\partial_{x_n}(\rho_y^{-n}) = -n x_n\rho_y^{-n-2}.
\end{align*}
Using the product rule for $\Delta(x_n\rho_y^{-n})$,
\begin{align*}
\Delta(x_n\rho_y^{-n}) = x_n\Delta(\rho_y^{-n}) + 2\partial_{x_n}(\rho_y^{-n}) = 2n x_n\rho_y^{-n-2} - 2n x_n\rho_y^{-n-2} = 0.
\end{align*}
Multiplying by the constant $2/(n\omega_n)$ gives
\begin{align*}
\Delta_x P(x,y) = 0.
\end{align*}
Thus $K_y$ is harmonic in $\mathbb{R}^n_+$.
[guided]
Fix $y\in\mathbb{R}^{n-1}$. We regard $P(x,y)$ as a function of the interior variable $x\in\mathbb{R}^n_+$. Define the map $K_y: \mathbb{R}^n_+\to\mathbb{R}$ by $K_y(x)=P(x,y)$. Define the displacement map $X_y: \mathbb{R}^n_+\to\mathbb{R}^n$ by $X_y(x)=(x'-y,x_n)$, and define the radius function $\rho_y: \mathbb{R}^n_+\to(0,\infty)$ by $\rho_y(x)=|X_y(x)|$. Since $x_n>0$, $X_y(x)\neq 0$ for every $x\in\mathbb{R}^n_+$, so the powers of $\rho_y$ used below are smooth functions on the half-space.
With this notation,
\begin{align*}
K_y(x)=\frac{2}{n\omega_n}x_n\rho_y(x)^{-n}.
\end{align*}
We verify the Laplacian identity rather than quote it. Write $X_{y,i}(x)$ for the $i$-th coordinate of $X_y(x)$. For each $i\in\{1,\dots,n\}$, the chain rule gives
\begin{align*}
\partial_{x_i}(\rho_y^{-n})=-nX_{y,i}\rho_y^{-n-2}.
\end{align*}
Differentiating once more gives
\begin{align*}
\partial_{x_i}^2(\rho_y^{-n})=-n\rho_y^{-n-2}+n(n+2)X_{y,i}^2\rho_y^{-n-4}.
\end{align*}
When we sum over all $i$, the identity $\sum_{i=1}^n X_{y,i}^2=\rho_y^2$ yields
\begin{align*}
\Delta(\rho_y^{-n})=2n\rho_y^{-n-2}.
\end{align*}
Also, differentiating $\rho_y^{-n}$ in the vertical coordinate and using $X_{y,n}(x)=x_n$ gives
\begin{align*}
\partial_{x_n}(\rho_y^{-n})=-nx_n\rho_y^{-n-2}.
\end{align*}
The point of writing the kernel as $x_n\rho_y^{-n}$ is that the product rule produces exactly two terms, and their signs cancel. Namely,
\begin{align*}
\Delta(x_n\rho_y^{-n})=x_n\Delta(\rho_y^{-n})+2\partial_{x_n}(\rho_y^{-n})=0.
\end{align*}
Multiplying by the constant $2/(n\omega_n)$ preserves the equality, hence
\begin{align*}
\Delta_xP(x,y)=0.
\end{align*}
Therefore $K_y$ is harmonic on $\mathbb{R}^n_+$.
[/guided]
[/step]
[step:Differentiate under the integral to prove harmonicity of the Poisson integral]
Let $K \subset \mathbb{R}^n_+$ be compact. Define
\begin{align*}
a_K := \operatorname{dist}(K,\partial\mathbb{R}^n_+) > 0.
\end{align*}
For every multi-index $\alpha \in \mathbb{N}_0^n$ with $|\alpha| \leq 2$, differentiating the rational function $x_n(|x'-y|^2+x_n^2)^{-n/2}$ finitely many times expresses $\partial_x^\alpha P(x,y)$ as a finite sum of terms whose numerators are polynomials in $x_n$ and $x'-y$ and whose denominators are powers of $|x'-y|^2+x_n^2$.
We now record the domination constant. Since $K$ is compact, there is $R_K>0$ such that $|x'|\leq R_K$ and $a_K\leq x_n\leq R_K$ for every $x\in K$. On the bounded region $|y|\leq 2R_K+1$, continuity of each derivative on the compact set $K\times \overline{B}(0,2R_K+1)$ gives a finite bound, and this bound is absorbed into $C_{K,\alpha}(1+|y|)^{-n}$ after enlarging $C_{K,\alpha}$. On the region $|y|>2R_K+1$, the inequalities $|x'-y|\geq |y|/2$ and $x_n\leq R_K$ show that each finite-sum term is bounded by a constant depending only on $K$ and $\alpha$ times $|y|^{-n}$. Therefore there is a constant $C_{K,\alpha}>0$ such that
\begin{align*}
|\partial_x^\alpha P(x,y)| \leq C_{K,\alpha}(1+|y|)^{-n}
\end{align*}
for all $x \in K$ and all $y \in \mathbb{R}^{n-1}$. Define
\begin{align*}
Y: \mathbb{R}^{n-1} \to [0,\infty), \quad y \mapsto (1+|y|)^{-n}.
\end{align*}
The function $Y$ belongs to $L^1(\mathbb{R}^{n-1})$ because $n > n-1$. Since $g$ is bounded, the functions $\partial_x^\alpha P(x,y)g(y)$ are dominated on $K \times \mathbb{R}^{n-1}$ by $C_{K,\alpha}\|g\|_\infty Y(y)$, which is integrable with respect to $\mathcal{L}^{n-1}$.
By the [Dominated Convergence Theorem](/theorems/4) applied to the difference quotients in each coordinate direction, with the integrable domination just established, differentiation under the integral sign is valid for all $|\alpha| \leq 2$, and $u \in C^2(\mathbb{R}^n_+)$. Using the harmonicity of each kernel slice from the previous step,
\begin{align*}
\Delta u(x) = \int_{\mathbb{R}^{n-1}} \Delta_x P(x,y)g(y)\,d\mathcal{L}^{n-1}(y) = 0
\end{align*}
for every $x \in \mathbb{R}^n_+$. Hence $u$ is harmonic in $\mathbb{R}^n_+$.
[guided]
The issue in this step is not the pointwise identity $\Delta_xP(x,y)=0$; that was proved for each fixed $y$. The issue is whether we are allowed to move the derivatives through the boundary integral defining $u$. Fix a compact set $K\subset\mathbb{R}^n_+$ and define
\begin{align*}
a_K := \operatorname{dist}(K,\partial\mathbb{R}^n_+) > 0.
\end{align*}
For every multi-index $\alpha\in\mathbb{N}_0^n$ with $|\alpha|\leq 2$, differentiating the rational function $x_n(|x'-y|^2+x_n^2)^{-n/2}$ produces finitely many terms with polynomial numerators in $x_n$ and $x'-y$ and denominator a power of $|x'-y|^2+x_n^2$.
Because $K$ is compact, choose $R_K>0$ such that $|x'|\leq R_K$ and $a_K\leq x_n\leq R_K$ for all $x\in K$. On the bounded boundary region $|y|\leq 2R_K+1$, each derivative $\partial_x^\alpha P$ is continuous on the compact set $K\times\overline{B}(0,2R_K+1)$, hence bounded there. On the exterior region $|y|>2R_K+1$, the triangle inequality gives $|x'-y|\geq |y|/2$, and the denominator estimate for each finite-sum term gives a bound by a constant depending only on $K$ and $\alpha$ times $|y|^{-n}$. Enlarging the constant to include the bounded region, there is $C_{K,\alpha}>0$ such that
\begin{align*}
|\partial_x^\alpha P(x,y)| \leq C_{K,\alpha}(1+|y|)^{-n}
\end{align*}
for all $x\in K$ and all $y\in\mathbb{R}^{n-1}$.
Define
\begin{align*}
Y: \mathbb{R}^{n-1}\to[0,\infty), \quad y\mapsto (1+|y|)^{-n}.
\end{align*}
The function $Y$ lies in $L^1(\mathbb{R}^{n-1})$ because the decay exponent $n$ is strictly larger than the boundary dimension $n-1$. Since $g$ is bounded, $|\partial_x^\alpha P(x,y)g(y)|\leq C_{K,\alpha}\|g\|_\infty Y(y)$, and the right-hand side is integrable with respect to $\mathcal{L}^{n-1}$. These are precisely the domination hypotheses needed to apply the [Dominated Convergence Theorem](/theorems/4) to the difference quotients in each coordinate direction. Therefore differentiation under the integral sign is valid through second order, so $u\in C^2(\mathbb{R}^n_+)$ and
\begin{align*}
\Delta u(x)=\int_{\mathbb{R}^{n-1}}\Delta_xP(x,y)g(y)\,d\mathcal{L}^{n-1}(y).
\end{align*}
The previous step proved $\Delta_xP(x,y)=0$ for every $x\in\mathbb{R}^n_+$ and every $y\in\mathbb{R}^{n-1}$, hence
\begin{align*}
\Delta u(x)=0
\end{align*}
for every $x\in\mathbb{R}^n_+$. Thus $u$ is harmonic in the upper half-space.
[/guided]
[/step]
[step:Prove that the kernel mass away from the boundary point vanishes]
Fix $y_0 \in \mathbb{R}^{n-1}$ and $\delta > 0$. Define the exterior boundary region
\begin{align*}
E_\delta := \{y \in \mathbb{R}^{n-1} : |y-y_0| \geq \delta\}.
\end{align*}
Assume $x = (x',x_n) \in \mathbb{R}^n_+$ satisfies $|x' - y_0| < \delta/2$. Then every $y \in E_\delta$ satisfies
\begin{align*}
|x' - y| \geq |y-y_0| - |x'-y_0| \geq \delta/2.
\end{align*}
Hence
\begin{align*}
\int_{E_\delta} P(x,y)\,d\mathcal{L}^{n-1}(y) \leq \frac{2x_n}{n\omega_n}\int_{\{y: |x'-y| \geq \delta/2\}} |x'-y|^{-n}\,d\mathcal{L}^{n-1}(y).
\end{align*}
Using the polar-coordinate formula for Lebesgue measure in $\mathbb{R}^{n-1}$ centered at $x'$, which applies because the integrand is non-negative and radial about $x'$, the last integral is
\begin{align*}
(n-1)\omega_{n-1}\int_{\delta/2}^{\infty} r^{-2}\,d\mathcal{L}^1(r) = \frac{2(n-1)\omega_{n-1}}{\delta}.
\end{align*}
Therefore
\begin{align*}
\int_{E_\delta} P(x,y)\,d\mathcal{L}^{n-1}(y) \leq \frac{4(n-1)\omega_{n-1}}{n\omega_n\delta}x_n.
\end{align*}
Thus
\begin{align*}
\lim_{x=(x',x_n)\to(y_0,0),\, x_n>0}\int_{E_\delta} P(x,y)\,d\mathcal{L}^{n-1}(y) = 0.
\end{align*}
[/step]
[step:Combine continuity near $y_0$ with the vanishing tail estimate]
Fix $\varepsilon > 0$. Since $g$ is continuous at $y_0$, choose $\delta > 0$ such that
\begin{align*}
|g(y)-g(y_0)| < \varepsilon
\end{align*}
whenever $|y-y_0| < \delta$. Define
\begin{align*}
M := \sup_{y \in \mathbb{R}^{n-1}} |g(y)|.
\end{align*}
The number $M$ is finite because $g$ is bounded. For $x = (x',x_n) \in \mathbb{R}^n_+$, the normalization of the kernel gives
\begin{align*}
u(x)-g(y_0) = \int_{\mathbb{R}^{n-1}} P(x,y)(g(y)-g(y_0))\,d\mathcal{L}^{n-1}(y).
\end{align*}
Splitting the integral into the set $\{y: |y-y_0| < \delta\}$ and its complement $E_\delta$, we obtain
\begin{align*}
|u(x)-g(y_0)| \leq \varepsilon\int_{\mathbb{R}^{n-1}}P(x,y)\,d\mathcal{L}^{n-1}(y) + 2M\int_{E_\delta}P(x,y)\,d\mathcal{L}^{n-1}(y).
\end{align*}
By normalization, the first integral equals $1$. By the previous step, the second integral tends to $0$ as $x=(x',x_n)\to(y_0,0)$ with $x_n>0$. Therefore
\begin{align*}
\limsup_{x=(x',x_n)\to(y_0,0),\, x_n>0} |u(x)-g(y_0)| \leq \varepsilon.
\end{align*}
Since $\varepsilon > 0$ was arbitrary,
\begin{align*}
\lim_{x=(x',x_n)\to(y_0,0),\, x_n>0} u(x) = g(y_0).
\end{align*}
This completes the proof.
[guided]
Fix $\varepsilon>0$. Since $g$ is continuous at $y_0$, choose $\delta>0$ such that
\begin{align*}
|g(y)-g(y_0)|<\varepsilon
\end{align*}
whenever $|y-y_0|<\delta$. Define
\begin{align*}
M:=\sup_{y\in\mathbb{R}^{n-1}}|g(y)|.
\end{align*}
This number is finite because $g$ is bounded. The normalization of the kernel lets us subtract $g(y_0)$ inside the integral: since $\int_{\mathbb{R}^{n-1}}P(x,y)\,d\mathcal{L}^{n-1}(y)=1$,
\begin{align*}
u(x)-g(y_0)=\int_{\mathbb{R}^{n-1}}P(x,y)(g(y)-g(y_0))\,d\mathcal{L}^{n-1}(y).
\end{align*}
Now split the boundary into the near set $\{y:|y-y_0|<\delta\}$ and the exterior set $E_\delta=\{y\in\mathbb{R}^{n-1}:|y-y_0|\geq\delta\}$. On the near set, continuity gives $|g(y)-g(y_0)|<\varepsilon$. On the exterior set, the boundedness estimate gives $|g(y)-g(y_0)|\leq 2M$. Therefore
\begin{align*}
|u(x)-g(y_0)|\leq\varepsilon\int_{\mathbb{R}^{n-1}}P(x,y)\,d\mathcal{L}^{n-1}(y)+2M\int_{E_\delta}P(x,y)\,d\mathcal{L}^{n-1}(y).
\end{align*}
The first integral equals $1$ by normalization. The previous step proved that the second integral tends to $0$ as $x=(x',x_n)\to(y_0,0)$ with $x_n>0$. Hence
\begin{align*}
\limsup_{x=(x',x_n)\to(y_0,0),\,x_n>0}|u(x)-g(y_0)|\leq\varepsilon.
\end{align*}
Since $\varepsilon>0$ was arbitrary, this proves
\begin{align*}
\lim_{x=(x',x_n)\to(y_0,0),\,x_n>0}u(x)=g(y_0).
\end{align*}
Together with harmonicity in the interior, this completes the proof of the theorem.
[/guided]
[/step]
