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