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, $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 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, 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,E) = L^{1,2}(M,S)$$H^1(M,S) = L^{1,2}(M,S)$ of square-integrable sections of $E$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).