If I understand correctly, the Euclidean Dirac operator is given by
$$D=\sum_{i=1}^n e_i \frac{\partial}{\partial x_i},$$
where $e_i$ are bases for $Cl_{0,n}(\mathbb{R})$, i.e., the $n$-dimensional Clifford algebra with negative-definite signature over the reals (so $e_i^2=-1$), and $x_i$ are the corresponding coordinates. Several sources state that $D^2 = -\Delta_n$ where $\Delta_n$ is the standard Euclidean Laplace operator
$$\Delta_n = \sum_{i=1}^n \frac{\partial^2}{\partial x_i^2}.$$
When I write out $D^2 f$ explicitly for some function $f:\mathbb{R}^n \rightarrow \mathbb{R}$, scalar terms from the Laplacian certainly appear, e.g.,
$$e_1 \frac{\partial}{\partial x_1}\left( e_1 \frac{\partial}{\partial x_1} f \right) = e_1 \left( e_1 \frac{\partial^2}{\partial x_1^2}f + \left(\frac{\partial}{\partial x_1}e_1\right)\frac{\partial}{\partial x_1}f \right)=e_1^2 \frac{\partial^2}{\partial x_1^2}f = -\frac{\partial^2}{\partial x_1^2}f.$$
But I also end up with bivector cross terms that shouldn't be there:
$$e_1 \frac{\partial}{\partial x_1}\left( e_2 \frac{\partial}{\partial x_2} f \right) = e_1 \left( e_2 \frac{\partial^2}{\partial x_1 \partial x_2}f + \left(\frac{\partial}{\partial x_1}e_2\right)\frac{\partial}{\partial x_2}f \right)=e_1 e_2 \frac{\partial^2}{\partial x_1 \partial x_2}f = e_{12}\frac{\partial^2}{\partial x_1 \partial x_2}f.$$
Should I only be considering the scalar part of $D^2$, or am I simply doing something wrong here?