Is there a "more rigorous" derivation of the electrostatic boundary conditions? When I first saw a derivation of the electrostatic boundary conditions it wasn't quite rigorous. It was essentially the argument used by Griffiths in his book:

Suppose we draw a wafer-thin Gaussian pillbox, extending just barely over the edge in each direction. Gauss' law states that:
$$\oint_S \mathbf{E}\cdot d\mathbf{ a} = \dfrac{1}{\epsilon_0}Q_{\text{enc}}=\dfrac{1}{\epsilon_0}\sigma A,$$
where $A$ is the area of the pillbox lid. (If $\sigma$ varies from point to point or the surface is curved, we must pick $A$ to be extremly small.) Now, the sides of the pillbox contribute nothing to the flux, in the limit as the thickness $\epsilon$ goes to zero, so we are left with:
$$E^\perp_{\text{above}}-E^\perp_{\text{below}}=\dfrac{1}{\epsilon_0}\sigma,$$
where $E^\perp_{\text{above}}$ denotes the component of $\mathbf{E}$ that is perpendicular to the surface immediately above, and $E^\perp_{\text{below}}$ is the same, only just below the surface. For consistency, we let "upward" be the positive direction for both. Conclusion: The normal component of $\mathbf{E}$ is discontinuous by an amount $\sigma/\epsilon_0$ at any boundary.

One analogous argument is also used for the continuity of the tangential component.
Now, on the one hand this argument is intuitive and easy to follow. It allows one to have one intuition on what is going on. On the other hand, I find it a quite "hand waving" argument.
Is there a more rigorous, less hand waving way to derive the boundary conditions: both for the normal and tangential components?
I though on something along the lines of expanding $\mathbf{E}$ around some point on the boundary, but it didn't work very much.
How can one derive in a little more precise manner these boundary conditions?
 A: Today I think I can answer the question myself, and considering the renewed interest on it, I shall do so. If anything is found to be mistaken, comments are welcome.
The reason for the question is that Griffiths was imprecise, he gave the idea of the proof, but not the mathematical rigor details. First, let us state better what is the result we want to prove:

Theorem: Let $\mathbf{E}$ be an electrostatic field, i.e. in particular obeying Gauss's law in integral form. Let $S\subset \mathbb{R}^3$ be one charged surface with surface charge $\sigma : S\to \mathbb{R}$ and with normal $\mathbf{n}_0 : S\to T\mathbb{R}^3$. Then if $p\in S$ it holds
$$\lim_{\epsilon\to 0}\mathbf{E}(p+\epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p)-\mathbf{E}(p-\epsilon \mathbf{n}_0(p))\cdot\mathbf{n}_0(p)=\dfrac{\sigma(p)}{\epsilon_0}.$$

We are going to show this using Griffith's idea, but filling in the details. First let's construct the "wafer-thin gaussian pillbox" and use Gauss' law.
Take $p\in S$ a point on the surface. Let $\mathbf{n}_0 : S\to T\mathbb{R}^3$ be the surface's normal vector field. Then let $U\subset S$ be open containing $p$, and consider the following set of points of $\mathbb{R}^3$
$$\mathcal{D}_U(\epsilon)=\{q\in \mathbb{R}^3 : q = q_0 + \lambda \mathbf{n}(q_0),\quad q_0\in U,\lambda \in [-\epsilon,\epsilon]\}$$
In other words: pick a neighborhood of $p$ on the surface. Take every point of the neighborhood and "go out of $S$ from it" in both directions. This is the "gaussian pillbox".
Define also $\mathcal{D}_U^\lambda(\epsilon)$ to be the subset of $\mathcal{D}_U(\epsilon)$ corresponding to one specific $\lambda\in [-\epsilon,\epsilon]$. It should be clear that the area satisfies $A(\mathcal{D}_U^\lambda(\epsilon))=A(U)$.
Now, $\mathcal{D}_U(\epsilon)$ is a three-dimensional volume with a boundary $\partial \mathcal{D}_U(\epsilon)$. Let's analyze its boundary. By construction its boundary can be split up in three parts
$$\partial \mathcal{D}_U(\epsilon)= \mathcal{D}_U^\epsilon(\epsilon)\cup \mathcal{D}_U^{-\epsilon}(\epsilon)\cup \mathcal{W}_U(\epsilon)$$
the first and second terms are respectivelly the top and bottom lids and the third, $\mathcal{W}_U(\epsilon)$ is the wall made up by picking the line bounding $U$ and moving it straight up and down. The wall has area $A(\mathcal{W}_U(\epsilon))=2\ell\epsilon$, where $\ell$ is the lenght of the curve bounding $U$.
Bringing these together, and remembering the only charge inside of $\mathcal{D}_U(\epsilon)$ is that one from the surface, Gauss's law
$$\int_{\partial \mathcal{D}_U(\epsilon)} \mathbf{E}\cdot d\mathbf{a}=\dfrac{1}{\epsilon_0}\int_{U} \sigma dA$$
can be recast as
$$\int_{\partial \mathcal{D}_U^\epsilon(\epsilon)}\mathbf{E}\cdot d\mathbf{a}+\int_{\partial \mathcal{D}_U^{-\epsilon}(\epsilon)} \mathbf{E}\cdot d\mathbf{a} + \int_{\mathcal{W}_U(\epsilon)}\mathbf{E}\cdot d\mathbf{a}=\dfrac{1}{\epsilon_0}\int_U \sigma dA.$$
Now we apply the the mean value theorem for integrals on both sides. This yields $p^\pm \in \mathcal{D}_U^{\pm \epsilon}(\epsilon)$ on the top/bottom lids, $q\in \mathcal{W}_U(\epsilon)$ and $p'\in U$ such that the equation becomes
$$\mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) A(\mathcal{D}_U^{\epsilon}(\epsilon))-\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) A(\mathcal{D}_U^{-\epsilon}(\epsilon))+\mathbf{E}(q)\cdot \mathbf{n}_W(q) A(\mathcal{W}_U^{\epsilon}(\epsilon))=\dfrac{\sigma(p')}{\epsilon_0}A(U).$$
where we have used that the normals of $\mathcal{D}_U^{\pm \epsilon}(\epsilon)$ are just $\mathbf{n}_0$ copied over, and with direction flipped on the bottom lid, and the normal of the wall is $\mathbf{n}_W$. Inserting what we know about the areas gives
$$\mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) A(U)-\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) A(U)+\mathbf{E}(q)\cdot \mathbf{n}_W(q) 2\ell \epsilon =\dfrac{\sigma(p')}{\epsilon_0}A(U).$$
Finally taking $\epsilon\to 0$ yields
$$\lim_{\epsilon \to 0} \mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) -\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) =\dfrac{\sigma(p')}{\epsilon_0}.$$
Now notice that since $p^{\pm} \in \mathcal{D}^{\pm \epsilon}_U(\epsilon)$ it must be $p^+ = p_0 + \epsilon \mathbf{n}_0(p_0)$ and $p^- = p_0' - \epsilon \mathbf{n}_0(p_0')$. This in turn gives a result
$$\lim_{\epsilon \to 0} \mathbf{E}(p_0 + \epsilon \mathbf{n}_0(p_0))\cdot \mathbf{n}_0(p_0) -\mathbf{E}(p_0' - \epsilon \mathbf{n}_0(p_0'))\cdot \mathbf{n}_0(p_0') =\dfrac{\sigma(p')}{\epsilon_0}.$$
Up to now we have left $U$ arbitrary and thus it holds for any such $U$. Consider now $U = B(p,1/n) \cap S$ the open ball centered in $p$ of radius $1/n$. The results will depend on $n$ of course. Thus we shall have $p_0 = x_n$, $p_0' = y_n$ and $p' = z_n$. The equation then becomes
$$\lim_{\epsilon \to 0} \mathbf{E}(x_n + \epsilon \mathbf{n}_0(p_0))\cdot \mathbf{n}_0(x_n) -\mathbf{E}(y_n - \epsilon \mathbf{n}_0(y_n))\cdot \mathbf{n}_0(y_n) =\dfrac{\sigma(z_n)}{\epsilon_0}.$$
Finally take $n\to \infty$ shrinking down to $p$. Since $x_n,y_n,z_n\in B(p,1/n)$ these sequences all converge to $p$. Assuming smoothnes of everything to be able to exchange the limits, we obtain the result
$$\lim_{\epsilon \to 0} \mathbf{E}(p + \epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p) -\mathbf{E}(p- \epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p) =\dfrac{\sigma(p)}{\epsilon_0},$$
as we wished to prove.
