Let $f\in L^2\mathbb{R}^d$ such that $\nabla f$ (in the distributional sense) coincides with an $L^2$-function outside of $0$. In which dimensions $d$ do we automatically have that $\nabla f\in (L^2\mathbb{R}^d)^d$?
For $d=1$ it is wrong, e.g. the Heaviside functions provides a counterexample.
For $d\ge3$, I can prove the statement: We have equality of $\nabla f$ with some $L^2$-function when testet against a test function which vanishes at $0$. In order to extend this result to arbitrary test functions, I want to use a collection of cutoff functions $\phi_\epsilon\colon \mathbb{R}^d\rightarrow B(0,1)$ with support in $B(0,\epsilon)$. The crucial point that would allow to pass to the limit $\epsilon \rightarrow 0$ is that $$ \nabla\phi_\epsilon \rightarrow 0 \quad \text{weakly in } L^2. $$ For that I considered $\psi(t)=\exp(-t^2/(1-t^2))$ and $\phi_\epsilon(x)=\psi(\vert x \vert/\epsilon)$, which is the standard way of producing bumps. Then $$ \vert \nabla\phi_\epsilon(x)\vert \lesssim 1/\epsilon, \quad \text{supp}(\nabla \phi_\epsilon)\subset B(0,\epsilon), $$ hence for any $L^2$-function $g$ we have $$ \vert\int g \nabla \phi_\epsilon\vert \lesssim \frac{1}{\epsilon}\int_{B(0,\epsilon)} \vert g\vert \le \frac{1}{\epsilon} \Vert g\Vert_{L^2}\cdot \vert B(0,\epsilon)\vert^{1/2} = O(\epsilon^{d/2 -1}), $$ thus the result follows if $d \ge 3$.
Question: Is the result true for $d=2$ or does there exist a counterexample?
Your final estimate is a bit rough. You can refine it as follows:
$$ \frac{1}{\epsilon}\int_{B_\epsilon} |g|dx \leq \frac{1}{\epsilon} \|g\|_{L^2(B_\epsilon)} \cdot |B_\epsilon|^{1/2}. $$ Now, for $d=2$, you have something that is bounded multiplied with $\|g\|_{L^2(B_\epsilon)}$. But the dominated convergence theorem shows $\|g\|_{L^2(B_\epsilon)} \to 0$ as $\epsilon \to 0$ (to see this, rewrite everything using indicator functions, so that the integrand converges to 0 pointwise a.e.).
This should allow you to get a positive answer also for $d=2$.