Factoring the Laplace operator $\Delta$ in dimensions $D \geq 3$ Consider the Laplace operator in 2 dimensions
\begin{equation}
\Delta = \frac{\partial^2}{\partial  x^2} + \frac{\partial^2}{\partial  y^2} = \partial^2_x + \partial^2_y
\end{equation}
By defining the complex variables $z$ and $\bar{z}$
\begin{equation}
\begin{split}
z &= x+iy\\
\bar{z} &= x-iy
\end{split}
\end{equation}
we can "factor" the Laplace operator into the product of 2 first order differential operators using freshman calculus
\begin{equation}
\Delta = 4\ \partial_z \ \partial_{\bar{z}}
\end{equation}
The other day I was watching the lecture Part 1 | Introduction to conformal field theory: Liouville model | Leon Takhtajan | Лекториум, when the lecturer remarked that this "factoring" property of the Laplace operator in 2 dimensions lies in the heart of Conformal Field Theory and Complex Analysis. He then said that this factoring bussiness can also be realized in dimensions $D \geq 3$, albeit in a much more involved manner, using spinors and Dirac matrices.
My question is, how can we factor the Laplace operator in dimensions $D \geq 3$, using spinors and Dirac matrices? What does this tell us about Conformal Field Theories in $D \geq 3$?
 A: @AccidentalFourierTransform drew your attention to the Dirac operator, which, of course, can go Euclidean and extend to all dimensions. It has little to do with conformal FT, if that's where you want to go, because the conformal group is finite in all dimensions different than d=2. So disaggregate the conformal gig.
First appreciate that Leon's construction can be superfluously replicated to 2×2 matrices (Pauli, duh!), through defining the Euclidean 2d Dirac operator,
$$
D=\sigma _x\partial_x + \sigma_y \partial_y =    \begin{pmatrix}
      0&\partial_x -i\partial_y\\
      \partial_x+i\partial_y&0
    \end{pmatrix} ~~~\leadsto \\
D^2 = \Delta_2 ~~{\mathbb 1}_2. 
$$
So far superfluous, but you may extend this to 3d, again with 2×2 matrices since their spinors are the same for an even dimension and its odd-one-higher one,
$$
D=\sigma _x\partial_x + \sigma_y \partial_y + \sigma_z \partial_z =    \begin{pmatrix}
      \partial_z&\partial_x -i\partial_y\\
      \partial_x+i\partial_y&-\partial_z
    \end{pmatrix} ~~~\leadsto \\
D^2 = \Delta_3 ~~{\mathbb 1}_3.
$$
And so on. The link provided will let you write down the Dirac gamma matrices in all dimensions, rectify your Minkowski metric by multiplying the spacelike ones with i, and dot them with the d-dimensional gradient to factorize your Laplacian in all dimensions, as above. For instance, for d =4, $D= -\vec \partial \cdot ~ (\sigma_2\otimes\vec \sigma )+ \partial_w ~\sigma_1\otimes {\mathbb 1}_2$, hence, yet again, $D^2 = \Delta_4 ~~{\mathbb 1}_4$.

*

*Leon's "much more involved" refers to the tensor product structure of coordinate space with spinor space. This latter part you cannot eschew, and it is a misunderstanding to expect it to be absent from the picture.

I'm not clear what you imagine you could get for CFTs in all dimensions this way. (To me the gig looks useless, but I should not aim to discourage creative thinking...)
A: In general, sums of squares can be factored for 2 dimensions (using complex numbers), for 4 dimensions (using quaternions), and for 8 dimensions (using octonions). The two factors are thereby conjugate to each other. The foundation of this are the three algebraic sum-of-squares identities (there are only three of them):
The product of two sums of two squares is again a sum of two squares (Diophantes); the product of two sums of four squares is again a sum of four squares (Leonhard Euler); and the product of two sums of eight squares is again a sum of eight squares (Ferdinand Degen).
These sums-of-squares identities are the foundations of the complex numbers, the quaternions, and the octonions, respectively. The length of the product of complex numbers, quaternions and octonions is always equal to the product of their lengths, which makes differential calculus possible in these number systems.
Instead of using complex number systems, one can also write the result as vector equations, following the multiplication rule of the corresponding complex number systems.
Hence, factoring a 4-dimensional Laplace operator yields two systems of each four coupled first-order differential equations:
$$\Delta = (\partial ^2_0 + \partial ^2_1 + \partial ^2_2 + \partial ^2_3) = (\partial_0 + i\partial_1 + j\partial_2 + k\partial_3)(\partial_0 - i\partial_1 - j\partial_2 - k\partial_3)$$
We have thus (the associative law holds for quaternions; thus we can shift the parentheses):
$$\Delta\Psi = \{(\partial_0 + i\partial_1 + j\partial_2 + k\partial_3)(\partial_0 - i\partial_1 - j\partial_2 - k\partial_3)\}\Psi = 0 $$
$$\Delta\Psi = (\partial_0 + i\partial_1 + j\partial_2 + k\partial_3)\{(\partial_0 - i\partial_1 - j\partial_2 - k\partial_3)\Psi\} = 0 $$
$\Psi$ is a quaternion function, comprising four components.
Substituting now:
$$(\partial_0 - i\partial_1 - j\partial_2 - k\partial_3)\Psi\ = \Phi$$
wherein $\Phi$ is a quaternion function, too, we obtain:
$$(\partial_0 + i\partial_1 + j\partial_2 + k\partial_3)\Phi = 0  $$
These are two systems of each four coupled first-order differential equations, which can be written in vectorial form as:
$$ (\partial_0\Phi_0 - \partial_1\Phi_1 - \partial_2\Phi_2 - \partial_3\Phi_3) = 0  $$
$$ (\partial_1\Phi_0 + \partial_0\Phi_1 - \partial_3\Phi_2 + \partial_2\Phi_3) = 0  $$
$$ (\partial_2\Phi_0 + \partial_3\Phi_1 + \partial_0\Phi_2 - \partial_1\Phi_3) = 0  $$
$$ (\partial_3\Phi_0 - \partial_2\Phi_1 + \partial_1\Phi_2 + \partial_0\Phi_3) = 0  $$
And similarly for $\Psi$:
$$ (\ \ \ \partial_0\Psi_0 + \partial_1\Psi_1 + \partial_2\Psi_2 + \partial_3\Psi_3) = \Phi_0  $$
$$ (-\partial_1\Psi_0 + \partial_0\Psi_1 + \partial_3\Psi_2 - \partial_2\Psi_3) = \Phi_1  $$
$$ (-\partial_2\Psi_0 - \partial_3\Psi_1 + \partial_0\Psi_2 + \partial_1\Psi_3) = \Phi_2  $$
$$ (-\partial_3\Psi_0 + \partial_2\Psi_1 - \partial_1\Psi_2 + \partial_0\Psi_3) = \Phi_3  $$
Inserting the second equations into the first ones yields back the original Laplace equation:
$$ (\partial^2_0\Psi_0 + \partial^1_1\Psi_0 + \partial^2_2\Psi_0 + \partial^3_3\Psi_0) = 0  $$
$$ (\partial^2_0\Psi_1 + \partial^2_1\Psi_1 + \partial^2_2\Psi_1 + \partial^2_3\Psi_1) = 0  $$
$$ (\partial^2_0\Psi_2 + \partial^2_1\Psi_2 + \partial^2_2\Psi_2 + \partial^2_3\Psi_2) = 0  $$
$$ (\partial^2_0\Psi_3 + \partial^2_1\Psi_3 + \partial^2_2\Psi_3 + \partial^2_3\Psi_3) = 0  $$
The 3-dimensional Laplace equation is a particular case of the 4-dimensional one, whose factoring can be obtained from the 4-dimensional factoring by setting the partial derivatives of the missing dimension to zero. 
The 8-dimensional Laplace equation can also be factored in octonion algebra. Albeit the associative law does no longer generally hold for octonions, it still holds in any sub-algebra spanned by only two octonions, which is the case here ($ \partial \partial^* \Psi $ is associative). 
Factorings for the 5-, 6-, and 7-dimensional Laplace equation can be obtained from the 8-dimensional factoring, by setting the partial derivatives of the missing dimensions to zero.
The wave equation can be conformally mapped onto the Laplace equation through the change of variables to complex, quaternion, or octonion coordinates, respectively, and then factored accordingly.
QED
