Recall that a densely defined operator $T$ on a Hilbert space $H$ is *essentially self-adjoint* if and only if its minimal closure $\overline{T}$ is self-adjoint, if and only if $T$ is symmetric and has a unique extension to a self-adjoint operator (i.e., the minimal closure $\overline{T}$).

In general, you can consider the spin Dirac operator $D$ on a complete Riemannian spin manifold $(M,g)$ with spinor bundle $S$ (e.g., the Dirac operator you consider on Euclidean $\mathbb{R}^4$). By [a theorem of Wolf](https://math.berkeley.edu/~jawolf/publications.pdf/paper_050.pdf), $D$ (as well as its square $D^2$) is essentially self-adjoint on the subspace $C_c^\infty(M,S)$ of compactly supported smooth spinor fields. More generally, it follows from [a theorem of Chernoff](https://www.sciencedirect.com/science/article/pii/0022123673900037) that $D^k$ is essentially self-adjoint on $C_c^\infty(M,S)$ for every $k \in \mathbb{N}$.

What about the domain of the minimal closure $\overline{D}$? [By a theorem of Gromov and Lawson](http://www.numdam.org/item/PMIHES_1983__58__83_0.pdf), if the curvature of the spinor Levi-Civita connection is uniformly bounded in the appropriate sense, then the domain of $\overline{D}$ is the first $L^2$ Sobolev space $H^1(M,S) = L^{1,2}(M,S)$ of square-integrable spinor fields with square-integrable weak covariant derivatives (w.r.t. the spinor Levi-Civita connection). In your case, this means that the domain of the unique self-adjoint extension of the Dirac operator on Euclidean $\mathbb{R}^4$ is precisely the space $H^1(\mathbb{R}^4) \otimes \mathbb{C}^4$, where $H^1(\mathbb{R}^4)$ is the Sobolev space of square integrable functions on $\mathbb{R}^4$ with square-integrable weak first-order partial derivatives (or equivalently square-integrable weak gradient).