[proofplan] Fix two distinct poles $x,y \in \Omega$ and remove small closed balls around them. On the punctured domain the two functions $G(\cdot,x)$ and $G(\cdot,y)$ are harmonic, so Green's second identity converts the equality of Laplacians into a boundary integral. The outer boundary contribution vanishes because both Green functions have zero Dirichlet trace. The two small spherical boundary pieces converge to $-G(x,y)$ and $G(y,x)$, respectively, because the [standard fundamental solution](/theorems/26) has total inward flux $1$ through a small deleted sphere. Letting the radius tend to zero gives $G(y,x)-G(x,y)=0$. [/proofplan]

[step:Remove small disjoint balls around the two poles]Fix $x,y \in \Omega$ with $x \neq y$. Let $\operatorname{dist}(a,E) := \inf\{|a-w| : w \in E\}$ denote the Euclidean distance from a point $a \in \mathbb{R}^n$ to a set $E \subset \mathbb{R}^n$, let $\mathcal{L}^n$ denote $n$-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}^n$, and let $\mathcal{H}^{n-1}$ denote $(n-1)$-dimensional [Hausdorff measure](/page/Hausdorff%20Measure) on hypersurfaces in $\mathbb{R}^n$. Let $\mathbb{R}^n_0 := \mathbb{R}^n \setminus \{0\}$ denote the punctured Euclidean space. Let $\Gamma$ denote the [standard fundamental solution](/theorems/26) for $-\Delta$, viewed as the map \begin{align*} \Gamma: \mathbb{R}^n_0 \to \mathbb{R}. \end{align*} For $\xi \in \mathbb{R}^n_0$, the real number $\Gamma(\xi)$ denotes its value at $\xi$. It is normalized so that the inward normal derivative through a deleted sphere satisfies the flux identity used below. Define \begin{align*} d := \min\left\{\operatorname{dist}(x,\partial\Omega), \operatorname{dist}(y,\partial\Omega), \frac{|x-y|}{3}\right\}. \end{align*} Then $d>0$. For $0<r<d$, define the punctured domain \begin{align*} D_r := \Omega \setminus \left(\overline{B}(x,r) \cup \overline{B}(y,r)\right). \end{align*} The balls $\overline{B}(x,r)$ and $\overline{B}(y,r)$ are disjoint and contained in $\Omega$, so \begin{align*} \partial D_r = \partial\Omega \cup \partial B(x,r) \cup \partial B(y,r). \end{align*} Let $\nu_r: \partial D_r \to \mathbb{R}^n$ denote the outward unit normal to $D_r$. On the deleted sphere $\partial B(x,r)$ this normal is $\nu_r(z)=(x-z)/r$, and on $\partial B(y,r)$ it is $\nu_r(z)=(y-z)/r$. Define \begin{align*} u_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,x). \end{align*} Define \begin{align*} v_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,y). \end{align*} By the hypotheses on $G$, both $u_r$ and $v_r$ are $C^2$ on $D_r$ and harmonic there: \begin{align*} \Delta u_r = 0 \quad \text{and} \quad \Delta v_r = 0 \quad \text{on } D_r. \end{align*}[/step]

[guided]Fix distinct points $x,y \in \Omega$. We first choose the puncture radius small enough that the two deleted balls stay inside $\Omega$ and do not touch each other. Define \begin{align*} d := \min\left\{\operatorname{dist}(x,\partial\Omega), \operatorname{dist}(y,\partial\Omega), \frac{|x-y|}{3}\right\}. \end{align*} The distance terms are positive because $x,y \in \Omega$ and $\Omega$ is open, and $|x-y|/3>0$ because $x\neq y$; hence $d>0$. For $0<r<d$, define \begin{align*} D_r := \Omega \setminus \left(\overline{B}(x,r) \cup \overline{B}(y,r)\right). \end{align*} The choice $r<\operatorname{dist}(x,\partial\Omega)$ and $r<\operatorname{dist}(y,\partial\Omega)$ puts both closed balls inside $\Omega$, and the choice $r<|x-y|/3$ makes them disjoint. Therefore the boundary of the punctured domain is the union of the original boundary and the two deleted spherical boundaries: \begin{align*} \partial D_r = \partial\Omega \cup \partial B(x,r) \cup \partial B(y,r). \end{align*} Let $\nu_r: \partial D_r \to \mathbb{R}^n$ denote the outward unit normal to $D_r$. On the deleted sphere $\partial B(x,r)$, outward from $D_r$ points inward toward the centre $x$, so \begin{align*} \nu_r(z)=\frac{x-z}{r}. \end{align*} On $\partial B(y,r)$ the same orientation convention gives \begin{align*} \nu_r(z)=\frac{y-z}{r}. \end{align*} Now define the two functions to which Green's second identity will be applied: \begin{align*} u_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,x). \end{align*} Define also \begin{align*} v_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,y). \end{align*} The two poles have been removed from $D_r$, so the hypotheses on $G$ imply that $u_r$ and $v_r$ are $C^2$ on $D_r$. Since $G(\cdot,x)$ is harmonic away from $x$ and $G(\cdot,y)$ is harmonic away from $y$, we have \begin{align*} \Delta u_r = 0 \quad \text{and} \quad \Delta v_r = 0 \quad \text{on } D_r. \end{align*}[/guided]

[step:Apply Green's second identity on the punctured domain]The assumed validity of [Green's identities](/theorems/30) on punctured domains, including the existence of the required normal traces on each boundary component, applies to $u_r$ and $v_r$ on $D_r$. Since both functions are harmonic on $D_r$, it gives \begin{align*} 0 = \int_{D_r} \left(u_r(z)\Delta v_r(z)-v_r(z)\Delta u_r(z)\right)\,d\mathcal{L}^n(z) \end{align*} and hence \begin{align*} 0 = \int_{\partial D_r} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z). \end{align*} On $\partial\Omega$, the zero Dirichlet boundary trace gives $u_r=0$ and $v_r=0$, so the contribution from $\partial\Omega$ is zero. Therefore \begin{align*} 0 = I_x(r)+I_y(r), \end{align*} where \begin{align*} I_x(r) := \int_{\partial B(x,r)} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z) \end{align*} and \begin{align*} I_y(r) := \int_{\partial B(y,r)} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z). \end{align*}[/step]

[guided]We use the theorem statement's hypothesis that [Green's identities](/theorems/30) may be applied on punctured domains, including the normal traces on $\partial\Omega$, $\partial B(x,r)$, and $\partial B(y,r)$, and we apply it to the two maps \begin{align*} u_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,x) \end{align*} and \begin{align*} v_r: D_r \to \mathbb{R}, \quad z \mapsto G(z,y). \end{align*} The reason for deleting the two balls is that $u_r$ is singular at $x$ and $v_r$ is singular at $y$, while both functions are $C^2$ away from those points. On $D_r$, both singularities have been removed, so the hypotheses of Green's second identity are satisfied. Because $\Delta u_r=0$ and $\Delta v_r=0$ on $D_r$, the volume integral in Green's second identity is \begin{align*} \int_{D_r} \left(u_r(z)\Delta v_r(z)-v_r(z)\Delta u_r(z)\right)\,d\mathcal{L}^n(z)=0. \end{align*} Thus the boundary integral must vanish: \begin{align*} 0 = \int_{\partial D_r} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z). \end{align*} The boundary $\partial D_r$ has three pieces: the original boundary $\partial\Omega$, the small sphere $\partial B(x,r)$, and the small sphere $\partial B(y,r)$. On $\partial\Omega$ both Green functions have zero continuous boundary trace, so the integrand is zero there because each term contains one of the factors $u_r$ or $v_r$. Therefore only the two small spheres remain: \begin{align*} 0 = I_x(r)+I_y(r), \end{align*} with \begin{align*} I_x(r) := \int_{\partial B(x,r)} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z) \end{align*} and \begin{align*} I_y(r) := \int_{\partial B(y,r)} \left(u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)-v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\right)\,d\mathcal{H}^{n-1}(z). \end{align*} This identity is the whole mechanism of the proof: the outer boundary contributes nothing, and the two deleted balls record the values of one Green function at the other pole.[/guided]

[step:Compute the small sphere limit around $x$]Let \begin{align*} \sigma_{n-1} := \mathcal{H}^{n-1}(\partial B(0,1)). \end{align*} For non-negative functions $A,B:(0,d)\to[0,\infty)$, the notation $A(r)=O(B(r))$ as $r\downarrow0$ means that there are constants $C>0$ and $r_0\in(0,d)$ such that $A(r)\leq C B(r)$ for all $0<r<r_0$. The notation $A(r)=o(1)$ as $r\downarrow0$ means $\lim_{r\downarrow0} A(r)=0$. Near $x$, the singular expansion gives \begin{align*} u_r(z)=\Gamma(z-x)+R_x(z) \end{align*} with $R_x$ of class $C^1$ near $x$. Since the outward unit normal to $D_r$ on $\partial B(x,r)$ is $\nu_r(z)=(x-z)/r$, the fundamental solution satisfies \begin{align*} \frac{\partial}{\partial \nu_r}\Gamma(z-x)=\frac{1}{\sigma_{n-1}r^{n-1}} \end{align*} for $z \in \partial B(x,r)$, both for $n\geq 3$ and for $n=2$. Since $R_x$ is $C^1$ near $x$, there is a constant $M_x>0$ and a radius $r_x>0$ such that \begin{align*} \left|\frac{\partial R_x}{\partial \nu_r}(z)\right| \leq M_x \end{align*} for all $z \in \partial B(x,r)$ and $0<r<r_x$. The function $v_r(z)=G(z,y)$ is $C^1$ in a neighbourhood of $x$, because $x \neq y$. Hence \begin{align*} \lim_{r\downarrow 0} \sup_{z\in \partial B(x,r)} |v_r(z)-G(x,y)|=0 \end{align*} and there are constants $M_v>0$ and $r_v>0$ such that \begin{align*} \left|\frac{\partial v_r}{\partial \nu_r}(z)\right| \leq M_v \end{align*} for all $z \in \partial B(x,r)$ and $0<r<r_v$. On $\partial B(x,r)$, the singular part satisfies $|\Gamma(z-x)|=O(r^{2-n})$ when $n\geq 3$ and $|\Gamma(z-x)|=O(|\log r|)$ when $n=2$, while $R_x$ is bounded. Therefore the estimate for the product with the bounded normal derivative is \begin{align*} \left|\int_{\partial B(x,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)\right| = O(r^{2-n})O(1)O(r^{n-1})+O(r^{n-1}) \end{align*} for $n\geq 3$, which is $O(r)$, and \begin{align*} \left|\int_{\partial B(x,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)\right| = O(|\log r|)O(1)O(r)+O(r) \end{align*} for $n=2$, which tends to $0$. Hence \begin{align*} \lim_{r\downarrow 0} \int_{\partial B(x,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=0. \end{align*} Also, \begin{align*} \int_{\partial B(x,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=\int_{\partial B(x,r)} v_r(z)\frac{1}{\sigma_{n-1}r^{n-1}}\,d\mathcal{H}^{n-1}(z)+o(1). \end{align*} Since $\mathcal{H}^{n-1}(\partial B(x,r))=\sigma_{n-1}r^{n-1}$ and $v_r(z)\to G(x,y)$ uniformly on $\partial B(x,r)$, we estimate \begin{align*} \left|\frac{1}{\sigma_{n-1}r^{n-1}}\int_{\partial B(x,r)} v_r(z)\,d\mathcal{H}^{n-1}(z)-G(x,y)\right| \leq \sup_{z\in \partial B(x,r)} |v_r(z)-G(x,y)|. \end{align*} The right-hand side tends to $0$, so the last expression converges to $G(x,y)$. Hence \begin{align*} \lim_{r\downarrow 0} I_x(r) = -G(x,y). \end{align*}[/step]

[guided]We now compute the contribution from the small sphere around $x$. The map $u_r: D_r \to \mathbb{R}$, $z \mapsto G(z,x)$ has its singularity at $x$, while the map $v_r: D_r \to \mathbb{R}$, $z \mapsto G(z,y)$ is regular near $x$ because $x \neq y$. Near $x$, the singular expansion is \begin{align*} u_r(z)=\Gamma(z-x)+R_x(z), \end{align*} where $R_x$ is a $C^1$ function near $x$. On $\partial B(x,r)$ the outward unit normal to $D_r$ is $\nu_r(z)=(x-z)/r$, and the fundamental solution is normalized so that \begin{align*} \frac{\partial}{\partial \nu_r}\Gamma(z-x)=\frac{1}{\sigma_{n-1}r^{n-1}} \end{align*} for $z \in \partial B(x,r)$, where $\sigma_{n-1}=\mathcal{H}^{n-1}(\partial B(0,1))$. For non-negative functions $A,B:(0,d)\to[0,\infty)$, the notation $A(r)=O(B(r))$ as $r\downarrow0$ means that there are constants $C>0$ and $r_0\in(0,d)$ such that $A(r)\leq C B(r)$ for all $0<r<r_0$, and $A(r)=o(1)$ means $\lim_{r\downarrow0}A(r)=0$. The sign is important: $\nu_r$ points outward from the punctured domain, hence inward toward the deleted ball. Because $R_x$ is $C^1$ near $x$, its normal derivative is bounded on all sufficiently small spheres. Because $v_r$ is $C^1$ near $x$, the normal derivative $\partial v_r/\partial \nu_r$ is also bounded there. The term \begin{align*} \int_{\partial B(x,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z) \end{align*} therefore tends to zero: $|\Gamma(z-x)|=O(r^{2-n})$ for $n\geq 3$ and $|\Gamma(z-x)|=O(|\log r|)$ for $n=2$, while $\mathcal{H}^{n-1}(\partial B(x,r))=\sigma_{n-1}r^{n-1}$. Thus the integral is $O(r)$ for $n\geq 3$ and $O(r|\log r|)$ for $n=2$. For the other term, the boundedness of $\partial R_x/\partial \nu_r$ gives an $o(1)$ contribution after integration over $\partial B(x,r)$. Hence \begin{align*} \int_{\partial B(x,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=\int_{\partial B(x,r)} v_r(z)\frac{1}{\sigma_{n-1}r^{n-1}}\,d\mathcal{H}^{n-1}(z)+o(1). \end{align*} Because $v_r(z)=G(z,y)$ is continuous at $x$, it converges uniformly to $G(x,y)$ on $\partial B(x,r)$. The elementary estimate \begin{align*} \left|\frac{1}{\sigma_{n-1}r^{n-1}}\int_{\partial B(x,r)} v_r(z)\,d\mathcal{H}^{n-1}(z)-G(x,y)\right| \leq \sup_{z\in \partial B(x,r)} |v_r(z)-G(x,y)| \end{align*} shows that the averaged integral converges to $G(x,y)$. Since $I_x(r)$ contains this term with a minus sign, we obtain \begin{align*} \lim_{r\downarrow 0} I_x(r)=-G(x,y). \end{align*}[/guided]

[step:Compute the small sphere limit around $y$]Near $y$, the singular expansion gives \begin{align*} v_r(z)=\Gamma(z-y)+R_y(z) \end{align*} with $R_y$ of class $C^1$ near $y$. On $\partial B(y,r)$, the outward unit normal to $D_r$ is $\nu_r(z)=(y-z)/r$, so \begin{align*} \frac{\partial}{\partial \nu_r}\Gamma(z-y)=\frac{1}{\sigma_{n-1}r^{n-1}}. \end{align*} Since $R_y$ is $C^1$ near $y$, there are constants $M_y>0$ and $r_y>0$ such that \begin{align*} \left|\frac{\partial R_y}{\partial \nu_r}(z)\right| \leq M_y \end{align*} for all $z \in \partial B(y,r)$ and $0<r<r_y$. The function $u_r(z)=G(z,x)$ is $C^1$ in a neighbourhood of $y$, because $y \neq x$. Hence \begin{align*} \lim_{r\downarrow 0} \sup_{z\in \partial B(y,r)} |u_r(z)-G(y,x)|=0 \end{align*} and there are constants $M_u>0$ and $r_u>0$ such that \begin{align*} \left|\frac{\partial u_r}{\partial \nu_r}(z)\right| \leq M_u \end{align*} for all $z \in \partial B(y,r)$ and $0<r<r_u$. On $\partial B(y,r)$, the singular part satisfies $|\Gamma(z-y)|=O(r^{2-n})$ when $n\geq 3$ and $|\Gamma(z-y)|=O(|\log r|)$ when $n=2$, while $R_y$ is bounded. Therefore \begin{align*} \left|\int_{\partial B(y,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)\right| = O(r^{2-n})O(1)O(r^{n-1})+O(r^{n-1}) \end{align*} for $n\geq 3$, which is $O(r)$, and \begin{align*} \left|\int_{\partial B(y,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)\right| = O(|\log r|)O(1)O(r)+O(r) \end{align*} for $n=2$, which tends to $0$. Hence \begin{align*} \lim_{r\downarrow 0} \int_{\partial B(y,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=0. \end{align*} Also, \begin{align*} \int_{\partial B(y,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=\int_{\partial B(y,r)} u_r(z)\frac{1}{\sigma_{n-1}r^{n-1}}\,d\mathcal{H}^{n-1}(z)+o(1). \end{align*} Since $\mathcal{H}^{n-1}(\partial B(y,r))=\sigma_{n-1}r^{n-1}$ and $u_r(z)\to G(y,x)$ uniformly on $\partial B(y,r)$, we estimate \begin{align*} \left|\frac{1}{\sigma_{n-1}r^{n-1}}\int_{\partial B(y,r)} u_r(z)\,d\mathcal{H}^{n-1}(z)-G(y,x)\right| \leq \sup_{z\in \partial B(y,r)} |u_r(z)-G(y,x)|. \end{align*} The right-hand side tends to $0$, so the last expression converges to $G(y,x)$. Therefore \begin{align*} \lim_{r\downarrow 0} I_y(r) = G(y,x). \end{align*}[/step]

[guided]The computation around $y$ is the same boundary-flux computation with the roles of $u_r$ and $v_r$ exchanged, but we spell it out because the guided proof must stand alone. Near $y$, the function $v_r: D_r \to \mathbb{R}$, $z \mapsto G(z,y)$ has the singular expansion \begin{align*} v_r(z)=\Gamma(z-y)+R_y(z), \end{align*} where $R_y$ is $C^1$ near $y$. On $\partial B(y,r)$ the outward unit normal to $D_r$ is $\nu_r(z)=(y-z)/r$, so the normalization of $\Gamma$ gives \begin{align*} \frac{\partial}{\partial \nu_r}\Gamma(z-y)=\frac{1}{\sigma_{n-1}r^{n-1}}. \end{align*} Because $R_y$ is $C^1$ near $y$, $\partial R_y/\partial \nu_r$ is bounded on sufficiently small spheres. Because $u_r(z)=G(z,x)$ is $C^1$ near $y$, since $x\neq y$, the normal derivative $\partial u_r/\partial \nu_r$ is also bounded there. Therefore \begin{align*} \int_{\partial B(y,r)} v_r(z)\frac{\partial u_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z) \to 0. \end{align*} Indeed, the singular part contributes $O(r)$ for $n\geq 3$ and $O(r|\log r|)$ for $n=2$, while the regular part contributes $O(r^{n-1})$. For the remaining term, the singular normal derivative gives \begin{align*} \int_{\partial B(y,r)} u_r(z)\frac{\partial v_r}{\partial \nu_r}(z)\,d\mathcal{H}^{n-1}(z)=\int_{\partial B(y,r)} u_r(z)\frac{1}{\sigma_{n-1}r^{n-1}}\,d\mathcal{H}^{n-1}(z)+o(1). \end{align*} Since $u_r(z)=G(z,x)$ is continuous at $y$, it converges uniformly to $G(y,x)$ on $\partial B(y,r)$. The estimate \begin{align*} \left|\frac{1}{\sigma_{n-1}r^{n-1}}\int_{\partial B(y,r)} u_r(z)\,d\mathcal{H}^{n-1}(z)-G(y,x)\right| \leq \sup_{z\in \partial B(y,r)} |u_r(z)-G(y,x)| \end{align*} shows that the averaged integral converges to $G(y,x)$. Since $I_y(r)$ contains this term with a plus sign and the other term tends to $0$, we obtain \begin{align*} \lim_{r\downarrow 0} I_y(r)=G(y,x). \end{align*}[/guided]

[step:Let the puncture radius tend to zero and identify the two values]For every sufficiently small $r>0$, Green's second identity gave \begin{align*} 0 = I_x(r)+I_y(r). \end{align*} Taking the limit as $r\downarrow 0$ and using the two small-sphere limits yields \begin{align*} 0 = -G(x,y)+G(y,x). \end{align*} Thus \begin{align*} G(x,y)=G(y,x). \end{align*} Since $x,y \in \Omega$ with $x\neq y$ were arbitrary, the Dirichlet Green function is symmetric off the diagonal.[/step]

[guided]For every sufficiently small $r>0$, the boundary identity from Green's second identity is \begin{align*} 0=I_x(r)+I_y(r). \end{align*} The previous small-sphere computations give the two limits \begin{align*} \lim_{r\downarrow 0}I_x(r)=-G(x,y) \end{align*} and \begin{align*} \lim_{r\downarrow 0}I_y(r)=G(y,x). \end{align*} Because limits preserve finite sums, passing to the limit in $0=I_x(r)+I_y(r)$ yields \begin{align*} 0=-G(x,y)+G(y,x). \end{align*} Rearranging gives $G(x,y)=G(y,x)$. The points $x,y\in\Omega$ with $x\neq y$ were arbitrary, so the symmetry holds for every off-diagonal pair.[/guided]