An Electric Potential Glued to a Cube-Shaped Insulator to Replicate a Point Charge: Charge Distribution I have been going back over this problem with a friend for the better part of a day:
A potential is glued to a cube-shaped insulator so that outside of the insulator the field is the same as a point particle. How can we calculate the surface charge distribution and the volume charge distribution?
 A: Introduction
$\def\ph{\varphi}\def\vr{\vec{r}}\def\eps{\varepsilon}\def\rmr{{\rm r}}\def\pd{\partial}\def\l{\left}\def\r{\right}\def\ltag#1{\tag{#1}\label{#1}}\def\div{\operatorname{div}}\def\grad{\operatorname{grad}}\def\nR{{\mathbb{R}}}$
You do not need a combination of a surface charge and a volume charge to re-produce a point-charge field outside a cube.
The surface-charge distribution is uniquely determined if you just use that for your purpose. The construction of this surface charge distribution is described in the first part of this answer.
If you just use a volume charge distribution this is not uniquely determined. An example for that case is given in the second part of the answer.
Surface charge distribution
Let $V$ be the the open interior of the cube containing the origin. For points $\vr$ outside the quader we set the potential $\ph(\vr)$ equal to the potential
$$
\ph_Q(\vr) =\frac Q{4\pi\eps_0|r|}
$$
of a point charge $Q$ at the origin.
This gives you also the boundary condition
$$
\ph(\vr) = \ph_Q(\vr)\text{ for }\vr\in\partial V
$$
for the potential inside the quader. Solve the Dirichlet boundary value problem with the Laplacian equation $\Delta\ph(\vr)=0$ and with this boundary condition for the potential $\ph(r)$ at points $\vr\in V$ inside of the quader.
Let $\partial^+_\nu\ph$ and $\partial^-_\nu\ph$ be the limits of the outer normal derivatives of $\ph$ on $\partial V$ from the outside and the inside of $V$, respectively. The surface charge will be $\sigma(\vr) = -\eps_0(\partial^+_\nu\ph-\eps_{\rm r}\partial^-_\nu\ph)$. Why? Try to explain this with an integral over the surface of a small volume sitting on the surface.
Here we assume that the insulator has a homogeneous and isotropic relative permittivity $\eps_\rmr$.

The following picture shows a section of the solution for the Dirichlet problem in the cube $[-1,1]^3$. Only the part $\bigcup_{z\in[-1,0]}[0,-z]^2\times\{z\}$ of the cube has been modelled that generates the full cube through a sequence of reflections. Natural boundary conditions $\partial_\nu \ph=0$ were used at the symmetry planes. The Dirichtlet boundary condition at $z=-1$ is
$$\ph(x,y,z) = \frac1{\sqrt{x^2+y^2+1}}.$$

The domain has been modelled with gmsh and solution is computed with getdp.
The limit of the outer normal derivative of $\ph$ from the outside at the boundary $z=-1$ is
\begin{align}
\pd^+_\nu \ph(x,y,-1) &= -\pd_z \ph(x,y,-1) = \l.\frac{z}{\l(x^2+y^2+z^2\r)^{3/2}}\r|_{z=-1}\\
\ltag{derPotOutside}
&=\frac{-1}{\l(x^2+y^2+1^2\r)^{3/2}}
\end{align}
The following picture shows the field $\sigma:=-(\pd^+_\nu\ph-\pd^-_\nu\ph)$ in the boundary area $[0,1]^2\times\{-1\}$. Thereby, $\pd^+_\nu\ph$ is given by \eqref{derPotOutside} and $\pd^-_\nu\ph$ results from the numerical solution of the boundary value problem. Note, that this corresponds to a normalized solution of the original problem with $\eps_\rmr=1$.

One can check the numerical result $\sigma$ with the integral $\int_{(x,y)\in[0,1]^2} \sigma(x,y,-1) dA$ which must give $\frac\pi6$. This works within the bounds of numerical precision.
For reproduction of the results I link here the geometry definition and the problem definition for the FEM solver.

The normal component of the field strength is discontinuous at the charged surface but the potential is continuous.
This can be demonstrated with a uniformly charged circular disk in the (x,y)-plane with center at the origin and radius $R$. We denote the charge density with $\sigma$. The complement of the disk is free space with permitivity $\eps_0$. For the demonstration we calculate the potential on the z-axis:
\begin{align}
\ph(z) &= \frac\sigma{4\pi\eps_0}\int_{r=0}^R \frac{ 2\pi r dr}{\sqrt{r^2 + z^2}}\\
&=\frac\sigma{2\eps_0} \left[\sqrt{r^2+z^2}\right]_{r=0}^R\\
&=\frac\sigma{2\eps_0} \left(\sqrt{R^2+z^2}-|z|\right)
\end{align}
The potential can be continuously extended at $z=0$ with the value
\begin{align}
\ph(0) = \frac {\sigma R}{2\eps_0}.
\end{align}
Volume charge distribution
The volume charge distribution is not uniquely determined. For an instance if you have $\eps_\rmr=1$ you can put a ball with uniform volume charge into the cube. That is
$$
\rho(\vr) = \begin{cases}
\frac{Q}{\frac34\pi R^3}&\text{ for }|\vr|<R\\
0&\text{ else}
\end{cases}
$$ 
with $R$ sufficiently small such that the ball fits into the cube.
Another possible space charge distribution is
\begin{align}
\rho(\vr) = \begin{cases}
\frac{\frac87 Q}{\frac34\pi R^3}&\text{ for } \frac12 R<|\vr|<R,\\
0&\text{ else.}
\end{cases}
\end{align}
If the insulator is a quader with constant permitivity $\eps_\rmr\neq1$ then you have given:


*

*the potential $\ph$ on the boundary $\partial V$ of the quader by the outer potential field $\ph^+(\vr)=\frac{Q}{4\pi\eps_0|\vr|}$ of the point charge

*the normal derivative $\partial^-_\nu\ph$ at the boundary of the quader through the normal component of the charge displacement density $$\partial^-_\nu\ph(\vr) = \frac1{\eps_\rmr}\partial^+_\nu\ph(\vr).$$ 


The taks is now to find an at least two times differentiable scalar function $\ph:V\rightarrow\nR$ interpolating the required potential and the required normal derivative at $\partial V$. The Volume charge density can then be calculated by
$$
\rho(\vr) = -\eps_\rmr\eps_0 \Delta \ph(\vr).
$$

The actual problem for varying $\eps_r$ is
\begin{align}
\div\l(\eps(\vr) \grad\ph(\vr)\r) = -\rho.
\end{align}
In your case $\eps(\vr)$ is even discontinous such that you need the weak form of $\eps(\vr)$ and you end up with the problem:
\begin{align}
\int_{\nR^3} (\grad\delta\ph(\vr))\eps(\vr)\grad\ph(\vr) d V = \int_{\nR^3} \delta\ph(\vr)\rho(\vr) dV
\end{align}
for all test functions $\delta\ph$. This can be solved numerically.
I do not know if there are known Green functions for that problem.
