Method of Images - Point Charge with Semi-infinite Dielectric

I finally figured out what constraints are missing apart from the boundary conditions: All values of the imaginary charges should be constants, i.e. no dependence on the coordinates e.g. the radial component "r".
For example, solving the boundary conditions only gives the values of the imaginary charges as functions of {the source charge, source charge location, locations of imaginary charges, r }. If one imposes the constraint mentioned before, the solution that all distances are equal naturally drops 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.