In the notion of Jackson 4.4, assume the original charge located at $z = d > 0$. For the right-half plane, the image charge of the left-hand side ($\epsilon_2$) is $q_2$ while for the right-half plane, the image charge of the right-hand side ($\epsilon_1$) is $q_1$.

Let's first pretend that we don't know the exact position of $q_1$ and $q_2$ and therefore denote them as $z=d_1$ and $z=d_2$.

Matching all the boundary condition ($D_z$, $E_\rho$ continuous) yields $$ \frac{qd}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2d_2}{(\sqrt{\rho^2 + d_2^2})^3} = - \frac{q_1d_1}{\sqrt{(\rho^2 + d_1^2})^3} \\ \frac{1}{\epsilon_1} \left( - \frac{q}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2}{(\sqrt{\rho^2 + d_2^2})^3} \right) = \frac{1}{\epsilon_2} \left( - \frac{q_1}{\sqrt{(\rho^2 + d_1^2})^3} \right) $$

Looks like we are running out of criteria to solve these two equations. However, if we want the method of image charge to work, there's a hidden requirement: The set of $q_1, q_2, d_1, d_2$ should work across the entire region of interest, i.e. the final equation we are looking for should not contain $\rho$ and $z$ dependence, otherwise, it means the method of image charge fails because the image charge setup we presumed is simply wrong (or method of image charge is simply not applicable to this problem).

Back to this problem, we need to find a combination of $q_1, q_2, d_1, d_2$ that gets rid of $\rho$ and $z$ dependence. This is the third criterion besides the two boundary conditions. And apparently, setting $d_1=d_2=d$ works so all denominators cancel out.

In the notion of Jackson 4.4, assume the original charge located at $z = d > 0$. For the right-half plane, the image charge of the left-hand side ($\epsilon_2$) is $q_2$ while for the left-half plane, the image charge of the right-hand side ($\epsilon_1$) is $q_1$.

Let's first pretend that we don't know the exact position of $q_1$ and $q_2$ and therefore denote them as $z=d_1$ and $z=d_2$.

Matching all the boundary condition ($D_z$, $E_\rho$ continuous) yields $$ \frac{qd}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2d_2}{(\sqrt{\rho^2 + d_2^2})^3} = - \frac{q_1d_1}{\sqrt{(\rho^2 + d_1^2})^3} \\ \frac{1}{\epsilon_1} \left( - \frac{q}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2}{(\sqrt{\rho^2 + d_2^2})^3} \right) = \frac{1}{\epsilon_2} \left( - \frac{q_1}{\sqrt{(\rho^2 + d_1^2})^3} \right) $$

Looks like we are running out of criteria to solve these two equations. However, if we want the method of image charge to work, there's a hidden requirement: The set of $q_1, q_2, d_1, d_2$ should work across the entire region of interest, i.e. the final equation we are looking for should not contain $\rho$ and $z$ dependence, otherwise, it means the method of image charge fails because the image charge setup we presumed is simply wrong (or method of image charge is simply not applicable to this problem).

Back to this problem, we need to find a combination of $q_1, q_2, d_1, d_2$ that gets rid of $\rho$ and $z$ dependence. This is the third criterion besides the two boundary conditions. And apparently, setting $d_1=d_2=d$ works so all denominators cancel out.