Let $\Omega$ be a bounded set in $\mathbb{R}^3$ and $f,g:\mathbb{R}^3\to\mathbb{H}$ are two quaternion-valued function \begin{align*} f(x) &= f_0(x) + if_1(x) + jf_2(x) + kf_3(x) \\ g(x) &= g_0(x) + ig_1(x) + jg_2(x) + kg_3(x) \end{align*} where $i,j,k$ are basic quaternions and for all $m,n\in\{1,2,3\}$, $f_m, g_n:\mathbb{R}^3\to\mathbb{R}$ have continuous partial derivatives with respect to $x_1,x_2,x_3$. Let $D$ be Dirac operator defined by \begin{equation*} D := i\dfrac{\partial}{\partial x_1} + j\dfrac{\partial}{\partial x_2} + k\dfrac{\partial}{\partial x_3}. \end{equation*} I want to show that \begin{equation*} \iint_{\partial\Omega} f(x)\mathbf{n}g(x)dS = \iiint_{\Omega} (fD\cdot g + f\cdot Dg)dx_1 dx_2 dx_3 \end{equation*} where $n = in_1 + jn_2 + kn_3$ is an unit outer normal vector of $\partial\Omega$.
My intuition is to expand $f\mathbf{n}g$ in L.H.S and $fD\cdot g + f\cdot Dg$ in R.H.S then use Divergence theorem to show that the equation is true. But this process is very tedious and easily to make mistake. I also search for some property of Dirac operator (maybe $fD\cdot g + f\cdot Dg = D(fg)$?), but all I found is some other types of Dirac operator that I can't understand. Is there any elegant solution for this problem that don't require any high level knowlegde?
$ \newcommand\R{\mathbb R} \newcommand\Cl{\mathrm{Cl}} \newcommand\PD[2]{\frac{\partial#1}{\partial#2}} \newcommand\dd{\mathrm d} \newcommand\bound\partial \newcommand\grade[1]{\langle#1\rangle} $
Let $\{e_1, e_2, e_3\} \subset \R^3$ be an orthonormal basis of vectors in the Euclidean Clifford algebra $\Cl_3(\R)$, meaning $e_ie_j + e_je_i = 2\delta_{ij}$, and let $I = e_1e_2e_3$. Note that $I$ commutes with all of $\Cl_3(\R)$ and that $I^2 = -1$. The even subalgebra $\{x_0 - xI \;:\; x_0 \in \R,\: x \in \R^3\}$ of $\Cl_3(\R)$ is isomorphic to the quaternions, and we define $$ i := -e_1I = e_3e_2,\quad j := -e_2I = e_1e_3,\quad k := -e_3I = e_2e_1. $$ The vector derivative is $$ \nabla_x = e_1\PD{}{x_1} + e_2\PD{}{x_2} + e_3\PD{}{x_3} $$ where $x = x_1e_2 + x_2e_2 + x_3e_3$. Your operator $$ D_x = i\PD{}{x_1} + j\PD{}{x_2} + k\PD{}{x_3} $$ is then related to $\nabla_x$ by $$ D_x = -\nabla_xI,\quad \nabla_x = D_xI. $$
The fundamental theorem of geometric calculus states that $$ \int_{\bound\Omega}f(x)\,\dd^2x\,g(x) = \int_\Omega f(x)\,\dd^3x\nabla_x\,g(x), \tag{$*$} $$ where $\dd^2x$ is the bivector-valued surface measure and $\dd^3x$ the pseudoscalar-valued volume measure, and $\nabla_x$ acts on both $f$ and $g$. We may write $$ \dd^2x = |\dd^2x|In(x),\quad \dd^3x = |\dd^3x|I, $$ where $|\dd^2x|, |\dd^3x|$ are the usual scalar surface area and volume measures and $n(x)$ is the outward-pointing unit normal of $\bound\Omega$ at $x$. This is related to your quaternion $\mathbf n$ by $$ \mathbf n(x) = -n(x)I,\quad n(x) = \mathbf n(x)I. $$ Putting this all together in ($*$), we see $$ \int_{\bound\Omega}|\dd^2x|f(x)I\mathbf n(x)Ig(x) = \int_\Omega|\dd^3x|f(x)ID_xIg(x) $$$$ \implies \int_{\bound\Omega}|\dd^2x|f(x)\mathbf n(x)g(x) = \int_\Omega|\dd^3x|f(x)D_xg(x). $$ I've realized that I don't know what you mean by e.g. $f(x)\cdot D$. If you simply mean quaternion multiplication with $D$ acting on $f$, then this last equation is the result you want, since by the product rule we get $$ f(x)D_xg(x) = \dot f(x)\dot D_xg(x) + f(x)\dot D_x\dot g(x), $$ where the dots indicate what is being differentiated. You can verify this sort of identity by expanding in coordinates and using the product rule for the partial derivatives $\PD{}{x_i}$.
If you mean something like $f(x)\cdot Dg(x) = \mathrm{Re}[f(x)D]g(x)$ with $D$ acting on both $f$ and $g$, however, then the result you want is false (and I don't believe it's true if $D$ only acts on the function its dotted with either). We may write $$ f(x)D_xg(x) = (f(x)\cdot D_x + f(x)\times D_x)g(x) = f(x)(D_x\cdot g(x) + D_x\times g(x)) $$ where $A\cdot B = \mathrm{Re}[AB] = \grade{AB}_0$ is the scalar part of $AB$ and $A\times B = \frac12(AB - BA)$ is the commutator product. Hence $$ 2\int_{\bound\Omega}|\dd^2x|f(x)\mathbf n(x)g(x) = \int_\Omega|\dd^3x|\Bigl[f(x)\cdot D_x\:g(x) + f(x)\:D_x\cdot g(x)\Bigr] + \int_\Omega|\dd^3x|\Bigl[f(x)\times D_x\:g(x) + f(x)\:D_x\times g(x)\Bigr]. $$ In order to have the result you want, we would need $$ \int_\Omega|\dd^3x|\Bigl[f(x)\cdot D_x\:g(x) + f(x)\:D_x\cdot g(x)\Bigr] = \int_\Omega|\dd^3x|\Bigl[f(x)\times D_x\:g(x) + f(x)\:D_x\times g(x)\Bigr], $$ But choosing $f(x) = x_1i$ and $g(x) = 1$ gives $-\int_\Omega|\dd^3x| = 0$, which is false for arbitrary $\Omega$.