Greens function in EM with boundary conditions confusion

Good questions; I'm sure a lot of people are confused on this stuff (as I was the first time I used Jackson).

Essentially your confusion boils down to being careful to consider the following fact:

The Green's function for a particular boundary value problem depends on the boundary conditions.

In particular, let's say you have a Dirichlet boundary value problem. Then, as Jackson shows on page 39, the appropriate Green's function for such a boundary value problem must (a) satisfy Poisson's equation with a delta function source in that region and (b) vanish on the boundary (see eq. 1.43) of that region. If you can find a function that has these two properties in the region you are considering, then you have found the Green's function for the Dirichlet problem.

So if we consider the half space ($z>0$) with dirichlet boundary conditions at $z=0$, then we are looking for a function that satisfies Poisson's equation in the upper half space with unit source and which vanishes on the boundary, which in this case is $z=0$ plus the "boundary at infinity."

You can check yourself that the function $$ G(\mathbf x, \mathbf x') = \frac{1}{|\mathbf x - \mathbf x'|} - \frac{1}{|\mathbf x + \mathbf x'|} $$ has these properties. The intuition for this, and the reason why most people say you can just immediately write this down, is that the first term clearly satisfies Poisson's equation with unit source in the upper half space where $\mathbf x'$ is being taken to have $z>0$, and the second term corresponds to having a unit source in the lower half space with opposite sign. Our intuition about potentials of points charges indicates that this will cause that combination to vanish at $z=0$. Also, the Laplacian of the second term vanishes in the upper half space, so it doesn't affect the fact that in the upper half space, this function satisfies Poisson's equation with unit source at $\mathbf x'$.

I hope that was clear? I can definitely try to clean it up or expand on this. I know from personal experience that it's a confusing topic!

Addendum in response to comments.

Green's functions are associated with a set of two data (1) A region (2) boundary conditions on that region. The function $1/|\mathbf x-\mathbf x'|$ is the Green's function for (1) All of space with (2) Dirichlet boundary conditions. This is because it (a) satisfies Poisson's equation with unit source in that region and (b) vanishes at the boundary of that region (which in this case is at infinity). In general, for any region $R$, for Dirichlet boundary conditions, as long as we simply find a function $G(\mathbf x,\mathbf x')$ that, (a) satisfies Poisson's equation with a unit source placed in that region, and that (b) vanishes on the boundary of the region, then we have found the Green's function for that Dirichlet problem (by the definition of a Green's function).

When we are trying to find the Green's function for the Dirichlet problem on the upper half space, we first imagine putting a point charge in the upper half space so that condition (a) is satisfied, this leaves us with the function $1/|\mathbf x-\mathbf x'|$ where $\mathbf x'$ is a point in the upper half space. Then, we notice that although this function is an appropriate solution to Poisson's equation, it does not vanish for $x$ on the boundary, so this can't be the Green's function for this Dirichlet problem. We need to do something to this function which does not spoil the fact that it satisfies the unit source Poisson equation in the upper half space but such that the resulting function additionally satisfies the appropriate boundary condition.

So we ask ourselves "what can we do to this function so that (a) remains satisfied, but so that (b) is also satisfied in the upper half space. Well we notice that if we add any function that satisfies Laplace's equation (Laplacian equaling zero) in the upper half space to $1/|\mathbf x-\mathbf x'|$, then the resulting function will still satisfy (a).

Now what sorts of functions satisfy Laplace's equation in the upper half space?

The answer is any charge distribution whose charge density is only nonzero outside of the upper half space will create a potential that satisfies Laplace's equation in the upper half space.

So if we can find a charge distribution that when placed in the lower half space produces a potential that when added to $1/|\mathbf x-\mathbf x'|$ causes their sum to vanish on the boundary, then their sum will satisfy the properties required of a Green's function. This is where we notice that an "image" point charge will do exactly this!

All we are doing with this point charges is an intuitive way of finding a function that satisfies the appropriate mathematical properties (a) and (b) that a Green's function for a Dirichlet problem must satisfy.

Just wanted to add a final point to clear this up for you. The potential of a point charge and the Green's function for your problem are the same, up to the normalization constant. In a comment, you say this is by coincidence; It is not, it's physical!

Take your equation $\nabla^2 \Phi = -\frac{\rho}{\epsilon}$, and suppose you wanted the charge density in that equation to be that of a point charge. What would you put in for the charge density of a point charge? Well, it only exists at one point in space, and the way we express densities that exist at only one point is to use delta functions. So for a point charge, $\rho\sim q\,\delta(|{\bf x-x'}|)$. The exact proportionality constants don't matter so much as the fact that there is a delta function in the charge density. Now you can see that a point charge's charge density inserted into the equation get's you the differential equation for the Green's function, up to the constants like $\pi$ and $\epsilon$.

To reiterate, point charges have densities which are proportional to their charge and a delta function, so talking about Green's functions and responses to point charges is equivalent.