Hi @timwong,
I think my formulation is correct because it is mine. 
Joking aside, I do think that my formulation is correct. The Stackexchange answer is very helpful, but I note that the author has added comment:
I think that I’m missing an ε in the surface integrals. I’ll check it later.
I the missing \varepsilon was lost in the derivation of the weak form. Starting from the strong form (and using Einstein summation notation):
-(\varepsilon \varphi_{,i})_{,i} = \rho
Weak form:
-\int_\Omega {(\varepsilon \varphi_{,i})_{,i} \hat{\varphi} \, dx} = \int_\Omega {\rho \hat{\varphi} \,dx}
Product rule of derivatives:
-\int_\Omega {[(\varepsilon \varphi_{,i} \hat{\varphi})_{,i} - \varepsilon \varphi_{,i} \hat{\varphi}_{,i}] \, dx} = \int_\Omega {\rho \hat{\varphi} \,dx}
We can now imagine splitting the integral over \Omega into the sum of two integrals over \Omega_1 and \Omega_2. Applying the divergence theorem we have that the weak form is the sum of the following:
-\int_{\partial \Omega_1} {\varepsilon_1 \varphi_{,i} \hat{\varphi} n_i\, ds} + \int_{\Omega_1} { \varepsilon_1 \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega_1} {\rho \hat{\varphi} \,dx}
-\int_{\partial \Omega_2} {\varepsilon_2 \varphi_{,i} \hat{\varphi} n_i\, ds} + \int_{\Omega_2} { \varepsilon_2 \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega_2} {\rho \hat{\varphi} \,dx}
Let \Omega_1 be the circle carrying the surface charge, and \Omega_2 be the surroundings. Then \partial \Omega_1 = \Gamma, and the surface integral over \partial \Omega_2 can be split over the interior (\Gamma) and exterior (\partial \Omega) boundaries. The surface integral over \partial \Omega is 0 (because \varepsilon \varphi_{,i} = 0 where Neumann conditions hold, and \hat{\varphi}=0 where Dirichlet conditions hold), leaving:
-\int_{\Gamma} {\varepsilon_1 \varphi_{,i} \hat{\varphi} n_i\, ds} + \int_{\Omega_1} { \varepsilon_1 \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega_1} {\rho \hat{\varphi} \,dx}
\int_{\Gamma} {\varepsilon_2 \varphi_{,i} \hat{\varphi} n_i\, ds} + \int_{\Omega_2} { \varepsilon_2 \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega_2} {\rho \hat{\varphi} \,dx}
where we have introduced a sign change on the integral \int_{\Gamma} {\varepsilon_2 \varphi_{,i} \hat{\varphi} n_i\, ds} so that n_i can unambiguously refer to the outward normal on the interior boundary, as viewed from \Omega_1.
Summing these two equations we find:
-\int_{\Gamma} {(\varepsilon_1 - \varepsilon_2) \varphi_{,i} \hat{\varphi} n_i\, ds} + \int_{\Omega} { \varepsilon \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega} {\rho \hat{\varphi} \,dx}
Recognizing D_i^{(j)} = -\varepsilon_j \varphi_{,i} we write:
-\int_{\Gamma} {(D^{(2)}_i - D^{(1)}_i) \hat{\varphi} n_i\, ds} + \int_{\Omega} { \varepsilon \varphi_{,i} \hat{\varphi}_{,i} \, dx} = \int_{\Omega} {\rho \hat{\varphi} \,dx}
Or
\int_{\Omega} { \varepsilon \nabla\varphi \cdot \nabla\hat{\varphi} \, dx} = \int_{\Omega} {\rho \hat{\varphi} \,dx} + \int_{\Gamma} {(D^{(2)} - D^{(1)}) \cdot n \hat{\varphi}\, ds}
Finally, making use of the relation (D^{(2)} - D^{(1)})\cdot n = \sigma:
\int_{\Omega} { \varepsilon \nabla\varphi \cdot \nabla\hat{\varphi} \, dx} = \int_{\Omega} {\rho \hat{\varphi} \,dx} + \int_{\Gamma} {\sigma \hat{\varphi}\, ds}
(In this derivation, I think I might be using a different convention than you for the numbering in the jump condition. I am taking n to point from D^{(1)} to D^{(2)}.)
I think the discrepancy occurs in nicoguaro’s answer on Stackexchange when he says
Then, the weak form would be
\int_{\Omega} {\varepsilon \nabla \phi \nabla w \, d\Omega} - \int_{\partial \Omega} {w \nabla \phi \cdot \hat{n} \, d\Gamma} = \int_{\Omega} {f w \, d\Omega}
I believe this should read
\int_{\Omega} {\varepsilon \nabla \phi \nabla w \, d\Omega} - \int_{\partial \Omega} {w \varepsilon \nabla \phi \cdot \hat{n} \, d\Gamma} = \int_{\Omega} {f w \, d\Omega}
Note the additional \varepsilon in the surface integral, which would ultimately eliminate the factor of \frac{1}{\varepsilon_1 - \varepsilon_2}.
