The strategy is to verify the properties from the definition. Part (i) is immediate. Part (ii) uses the maximum principle. Part (iii) uses Green's second identity on a doubly-excised domain. Part (iv) follows from the [Green's Representation Formula](/theorems/567) by substituting the Green's function for the fundamental solution.

**Step 1: Proof of (i).** Since $\Phi(y - x)$ is harmonic in $y$ for $y \neq x$ and $\phi^x$ is harmonic in $y$ for all $y \in \Omega$, the difference $G(x,y) = \Phi(y-x) - \phi^x(y)$ is harmonic in $y$ on $\Omega \setminus \{x\}$. The boundary condition $G(x,y) = 0$ for $y \in \partial\Omega$ holds by definition of the corrector: $\phi^x(y) = \Phi(y-x)$ on $\partial\Omega$.

**Step 2: Proof of (ii).**

[claim:Positivity Of The Green Function] $G(x,y) > 0$ for all $x, y \in \Omega$ with $x \neq y$. [/claim]

[proof] Fix $x \in \Omega$. For $y$ near $x$, $\Phi(y-x) \to +\infty$ while $\phi^x(y)$ remains bounded (being harmonic, hence smooth, on $\Omega$), so $G(x,y) > 0$ in a neighbourhood of $x$. On $\partial\Omega$, $G(x,y) = 0$. The function $y \mapsto G(x,y)$ is harmonic on $\Omega \setminus \{x\}$ and continuous on $\overline{\Omega} \setminus \{x\}$. By the [Strong Maximum Principle](/theorems/32) applied to $-G$ (which is also harmonic), if $G$ achieved a non-positive minimum at an interior point $y_0 \in \Omega \setminus \{x\}$, then $G$ would be constant on the connected component containing $y_0$, contradicting $G(x,y) \to +\infty$ as $y \to x$. Hence $G(x,y) > 0$ for all $y \in \Omega \setminus \{x\}$. [/proof]

**Step 3: Proof of (iii).**

[claim:Symmetry Of The Green Function] $G(x,y) = G(y,x)$ for all $x, y \in \Omega$ with $x \neq y$. [/claim]

[proof] Fix distinct $x, z \in \Omega$. Choose $\varepsilon > 0$ small enough that $\overline{B(x,\varepsilon)}$, $\overline{B(z,\varepsilon)}$ are disjoint and contained in $\Omega$. Define $\Omega_\varepsilon := \Omega \setminus (\overline{B(x,\varepsilon)} \cup \overline{B(z,\varepsilon)})$. The functions $u(y) := G(x,y)$ and $v(y) := G(z,y)$ are both harmonic in $y$ on $\Omega_\varepsilon$. Green's second identity gives \begin{align*} \int_{\Omega_\varepsilon} (u\,\Delta v - v\,\Delta u)\,d\mathcal{L}^n &= \int_{\partial\Omega_\varepsilon} \Bigl(u\,\frac{\partial v}{\partial\nu} - v\,\frac{\partial u}{\partial\nu}\Bigr)\,d\mathcal{H}^{n-1}. \end{align*} The left-hand side is zero since $\Delta u = \Delta v = 0$ on $\Omega_\varepsilon$. The boundary $\partial\Omega_\varepsilon$ consists of three pieces: $\partial\Omega$, $\partial B(x,\varepsilon)$, and $\partial B(z,\varepsilon)$. On $\partial\Omega$: $u = G(x,\cdot) = 0$ and $v = G(z,\cdot) = 0$, so the integrand vanishes.

On $\partial B(x,\varepsilon)$ (with $\nu$ pointing into $B(x,\varepsilon)$): $u(y) = G(x,y) = \Phi(y-x) - \phi^x(y)$ has a singularity like $\Phi$ as $y \to x$, while $v(y) = G(z,y)$ is smooth near $x$. The same computation as in the proof of the [Green's Representation Formula](/theorems/567) gives, as $\varepsilon \to 0$: \begin{align*} \int_{\partial B(x,\varepsilon)} \Bigl(u\,\frac{\partial v}{\partial\nu} - v\,\frac{\partial u}{\partial\nu}\Bigr)\,d\mathcal{H}^{n-1} &\to -v(x) = -G(z,x). \end{align*} By the identical argument on $\partial B(z,\varepsilon)$ with the roles of $u$ and $v$ exchanged: \begin{align*} \int_{\partial B(z,\varepsilon)} \Bigl(u\,\frac{\partial v}{\partial\nu} - v\,\frac{\partial u}{\partial\nu}\Bigr)\,d\mathcal{H}^{n-1} &\to u(z) = G(x,z). \end{align*} Hence $0 = -G(z,x) + G(x,z)$, giving $G(x,z) = G(z,x)$. [/proof]

**Step 4: Proof of (iv).** Starting from the [Green's Representation Formula](/theorems/567): \begin{align*} u(x) &= -\int_{\partial\Omega} \Phi(x-y)\,\frac{\partial u}{\partial\nu}(y)\,d\mathcal{H}^{n-1} + \int_{\partial\Omega} u(y)\,\frac{\partial\Phi}{\partial\nu_y}(x-y)\,d\mathcal{H}^{n-1} - \int_\Omega \Phi(x-y)\,\Delta u(y)\,d\mathcal{L}^n. \end{align*} Since $\phi^x$ is harmonic in $\Omega$ and smooth on $\overline{\Omega}$, applying Green's second identity to $\phi^x$ and $u$ on $\Omega$ gives \begin{align*} 0 &= \int_{\partial\Omega} \Bigl(\phi^x\,\frac{\partial u}{\partial\nu} - u\,\frac{\partial\phi^x}{\partial\nu}\Bigr)\,d\mathcal{H}^{n-1} - \int_\Omega \phi^x(y)\,\Delta u(y)\,d\mathcal{L}^n. \end{align*} Subtracting the second equation from the first and using $G(x,y) = \Phi(y-x) - \phi^x(y)$ and $G(x,y) = 0$ on $\partial\Omega$ (which eliminates the $\frac{\partial u}{\partial\nu}$ boundary integral): \begin{align*} u(x) &= -\int_\Omega G(x,y)\,\Delta u(y)\,d\mathcal{L}^n(y) + \int_{\partial\Omega} u(y)\,\frac{\partial G}{\partial\nu_y}(x,y)\,d\mathcal{H}^{n-1}(y). \end{align*} Substituting $-\Delta u = f$ and $u = g$ on $\partial\Omega$ gives the claimed formula.