Denote by $\mathcal{D}_m$ the Lorentz covariant derivative,
$$\mathcal{D}_m=\partial_m-\frac{1}{2}\omega_m{}^{ab}M_{ab} \tag{1}$$
where indices $m,n,p,\dots$ are world indices, indices $a,b,c,\dots$ are tangent indices, $\omega_{m}{}^{ab}$ is the Lorentz connection and $M_{ab}$ are the Lorentz generators (in an arbitrary representation). I would like to compute the covariant derivative of the determinant of the vielbein, $e:=\text{det}(e_m{}^{a})$, in both the torsion and torsionless cases. Lets consider $D=3+1$ for concreteness.
In the torsionless case, it is known that the Lorentz covariant derivative $\mathcal{D}_m$ coincides with the regular covariant derivative $\nabla_m$ of general relativity (in tosion free case) involving the Christoffel symbols. Therefore, since we know that $\nabla_m\sqrt{-g}=0$ and $\sqrt{-g}=\sqrt{-\text{det}(g_{mn})}=$ $=\sqrt{-\text{det}(e_m{}^{a}e_n{}^{b}\eta_{ab})}=\sqrt{e^2}=|e|$, we can conclude that $\mathcal{D}_me=0$.
How would I go about computing $\mathcal{D}_me$ without using the knowledge $\nabla_m\sqrt{-g}=0$?
In the torsionless case, the spin connection could be expressed in terms of the vielbein via the torsion free condition, which would eliminate $\omega_{m}{}^{ab}$ in $(1)$ from the calculation in favour of the vielbein. Further more the torsion free condition is equivalent to the statement that the vielbein is covariantly constant, $\hat{\nabla}_me_{n}{}^{a}=0$ where $\hat{\nabla}_a=\nabla_m-\frac{1}{2}\omega_{m}{}^{bc}M_{bc}$, that might also come in handy. But these are just ideas of what might be useful in the calculation, i'm not sure how to put them to use. To start with, i'm not even sure how to interpret the expression $M_{bc}e$. Normally we would have $M_{bc}\phi=0$ for any scalar field $\phi$, but I suspect that this is not the case here (given that I anticipate $\mathcal{D}_me=0$ and $\partial_m e\neq 0$). I know that under an infinitesimal local Lorentz transformation parameterized by $K^{ab}=-K^{ba}$, $e$ is invariant, but I don't think I can apply that here.
Once I see how the torsionless case is done, I'll have a crack at the case with torsion myself.
Edit 1
Might someone confirm whether or not the following is valid:
\begin{align}
\frac{1}{2}\omega_n{}^{bc}M_{bc}\text{det}(e_{m}{}^{a})&=\frac{1}{2}\omega_n{}^{bc}M_{bc}\big(\frac{1}{4!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}e_{m_1}{}^{a_1}e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)\\
&=\frac{4}{4!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}\big(\frac{1}{2}\omega_{nbc}M^{bc}e_{m_1}{}^{a_1}\big)e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\\
&=\frac{1}{3!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}\big(\omega_{nbc}\eta^{a_1b}e_{m_1}{}^{c}\big)e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\\
&=-\frac{1}{3!}\omega_{nc}{}^{a_1}\varepsilon_{a_1a_2a_3a_4}\big(\varepsilon^{m_1m_2m_3m_4}e_{m_1}{}^{c}e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)\\
&=-\frac{1}{3!}\omega_{nc}{}^{a_1}\varepsilon_{a_1a_2a_3a_4}\big(\text{det}(e_{m}^{a})\varepsilon^{ca_2a_3a_4}\big)\\
&=-\frac{1}{3!}\cdot 3!e\omega_{nc}{}^{a_1}\delta^{c}_{a_1}\\
&=0
\end{align}
where in the last line I have used the antisymmetry $\omega_{n}{}^{bc}=-\omega_{n}{}^{cb}$. This would imply that $\mathcal{D}_me=\partial_me$ in both the torsion and tosionless case, which I am doubtful of (but coincidentally, it is exactly what I need for a separate calculation I am doing).
Edit 2
Similar to the above I can show that $\nabla_m e=0$,
\begin{align}
\nabla_m e &=\nabla_m\big(\frac{1}{4!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}e_{m_1}{}^{a_1}e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)\\
&=\partial_m\big(\frac{1}{4!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}e_{m_1}{}^{a_1}e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)~~~~~~~~~~~~~~~~~~~-4\big(\frac{1}{4!} \varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}(\Gamma^p_{mm_1}e_{p}{}^{a_1})e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)\\
&=\partial_me-\frac{1}{3!}\varepsilon^{m_1m_2m_3m_4}\Gamma^p_{mm_1} e\varepsilon_{pm_2m_3m_4}\\
&=\partial_me-e\Gamma^p_{mp} \\
&=\partial_me-e(\frac{1}{2}g^{-1}\partial_mg)\\
&=\partial_me-\partial_me\\
&=0
\end{align}
These two results (in edit 1 and 2) are consistent with the fact that the Vielbein is covariantly constant (true for torsionless case) since on the one hand we have (using $\hat{\nabla}_me_{n}{}^{a}=0$)
\begin{align}
\hat{\nabla}_me=\hat{\nabla}_m\big(\frac{1}{4!}\varepsilon^{m_1m_2m_3m_4}\varepsilon_{a_1a_2a_3a_4}e_{m_1}{}^{a_1}e_{m_2}{}^{a_2}e_{m_3}{}^{a_3}e_{m_4}{}^{a_4}\big)=0
\end{align}
and the the other hand we have
\begin{align}
\hat{\nabla}_me=(\nabla_m-\frac{1}{2}\omega_m{}^{bc}M_{bc})e=0.
\end{align}
However there is an inconsistency with the claim that $\mathcal{D}_m$ coincides with $\nabla_m$ in the torsionless case (since according to the above we have $\mathcal{D}_me=\partial_me$ and $\nabla_me=0$). Does anyone see where the problem is?
Best Answer
Notation:
Here in this answer I will break with some conventions present in the OP's question. In particular, greek $\mu,\nu...$ indices will denote components taken with respect to a holonomic (coordinate-) frame (also known as "world indices" or "curved indices") and latin $a,b...$ indices will denote local Lorentz indices (also known as "flat indices"). I will work in $n$ dimensions.
In addition to their question, OP is also confused/interested about the relationship between torsion and the covariant constancy of the volume element. For this reason we will use a connection $\nabla_\mu$ which is metric compatible, but not torsionless.
Preliminaries:
On a vector field with world indices, the connection acts as $$ \nabla_\mu X^\nu=\partial_\mu X^\nu + \Gamma^\nu_{\mu\rho}X^\rho, $$ where the coefficients $\Gamma^\nu_{\mu\rho}$ are not the usual Christoffel symbols, but the Christoffel symbols extended with the contorsion tensor.
On a vector field with local Lorentz indices, the connection acts as $$ \nabla_\mu X^a=\partial_\mu X^a+\omega_{\mu}{}^{a}{}_{b}X^b, $$ where the $\omega_{\mu}{}^{a}{}_{b}$ coefficients are not the usual Ricci rotation coefficients, but the Ricci rotation coefficients extended with terms involving torsion. Note that I will not use the Lorentz algebra generators directly, the $\omega_{\mu}{}^{a}{}_{b}$ coefficients are Lorentz algebra valued, and any other representations will be derived directly from these.
It should also be noted that technically these are the same connection, in two different representations.
The vielbein is $e^\mu_a$ and the covielbein is $\theta^a_{\mu}$. These are inverses of one another as matrices.
I will also note that the vielbein (and the covielbein) is covariantly constant, even if the connection is torsionful. In fact, they are covariantly constant even if the connection is not even metric-compatible.
Confusion can arise from the fact that the covariant derivative of the vielbein depends on which indices of the vielbein one considers as "tensor indices". If both indices are considered as tensor indices, then $\theta^a_\mu$ is just the identity tensor $\delta^\mu_\nu$ expressed with mixed indices. It's derivative vanishes as a consistency postulate: $$ \nabla_\mu\theta^a_\nu=\partial_\mu\theta^a_\nu+\omega_{\mu}{}^{a}{}_{b}\theta^b_\nu-\Gamma^\rho_{\mu\nu}\theta^a_\rho=0. $$
On the other hand, if a so-called covariant exterior derivative is calculated as $$ \nabla_\mu\theta^a_{(\nu)}-\nabla_\nu\theta^a_{(\mu)}=T^a_{\mu\nu}, $$ where the indices in parentheses are excempt from the covariant derivative (not considered as tensor indices), then the right hand side of this expression is the torsion 2-form.
The connection is assumed to be metric compatible, however, so we have $$ 0=\nabla_\mu g_{ab}=\partial_\mu g_{ab}-\omega_{\mu}{}^{c}{}_{a}g_{cb}-\omega_{\mu}{}^{c}{}_{b}g_{ac}=-(\omega_{\mu ba}+\omega_{\mu ab}), $$ where we have used that in a local Lorentz frame, $g_{ab}$ has constant components. Ergo, the connection is skew-symmetric in $ab$ even if it is torsionful, only metric compatibility is required.
Densities in vielbein formalism - covariant derivative of densities:
For the sake of simplicity, I will assume that the manifold we are working on is orientable, and a positive orientation has been chosen. I will express everything in either positive coordinates, or positive frames. This way, we can get rid of the absolute value in the transformation rule for densities.
Let $\pi_{\mu_1...\mu_n}$ denote the Levi-Civita symbol. It is not assumed to be a geometric object at all. A totally antisymmetric tensor $\rho_{\mu_1...\mu_n}$ has only one independent component, and it can be written as $$ \rho_{\mu_1...\mu_n}=\rho\cdot\pi_{\mu_1...\mu_n}, $$ where $\rho=\rho_{12...n}$.
It is not difficult to check that $\rho$ transforms as a scalar density of weight 1 (without the absolute value sign, but we got rid of it by choice of orientation).
Hence, we obtain a duality between scalar densities of weight 1, and totally antisymmetric covariant tensors of order $n$ ($n$-forms).
The canonical volume element on a pseudo-Riemannian manifold is defined as $$ \mu_{\mu_1...\mu_n}=\sqrt{|\det g|}\pi_{\mu_1...\mu_n}=|\det\theta|\pi_{\mu_1...\mu_n}. $$
Note the greek indices - this form of the volume tensor is taken with respect to a holonomic frame.
The form of the volume tensor in a local Lorentz frame can be calculated as $$ \mu_{a_1...a_n}=\mu_{\mu_1...\mu_n}e^{\mu_1}_{a_1}...e^{\mu_n}_{a_n}=|\det\theta|\pi_{\mu_1...\mu_n}e^{\mu_1}_{a_1}...e^{\mu_n}_{a_n}=\det\theta\det e \pi_{a_1...a_n}=\pi_{a_1...a_n}. $$ Therefore, if we define the single, independent component of the volume tensor as the function $f$ multiplying the Levi-Civita symbol $\pi_{a...}$, then we arrive at the important result that in a local Lorentz frame, the component of the volume density/tensor is 1.
The covariant derivative preserves symmetries, so if $\rho_{a_1...a_n}$ is a completely antisymmetric tensor field expressed in a local Lorentz frame, we can calculate its covariant derivative as $$ \nabla_\mu\rho_{a_1...a_n}=(\nabla_\mu\rho_{1...n})\pi_{a_1...a_n}, $$ so let us calculate $\nabla_\mu\rho_{1...n}$ for a generic $\rho$.
We have $$ \nabla_\mu\rho_{1...n}=\partial_\mu\rho_{1...n}-\omega_\mu{}^b{}_{1}\rho_{b2...n}-...-\omega_\mu{}^b{}_{n}\rho_{12...b} \\ =\partial_\mu\rho-\rho(\omega_\mu{}^b{}_{1}\pi_{b2...n}+\omega_\mu{}^b{}_{n}\pi_{12...b}) \\ =\partial_\mu\rho-\rho\omega_{\mu}{}^b{}_b=\partial_\mu\rho.$$
Due to the antisymmetry of $\pi$, in each contraction in the parentheses, only one term survives - the one where $b$ has the same value as the second lower index of $\omega$, so the entire parenthetical expression reduces to a trace, and we have utilized the fact that the connection coefficients have vanishing trace due to antisymmetry (which itself stems from metric-compatibility, as stated in the beginning).
We have obtained that when expressed in local Lorentz frames, densities can be partially differentiated.
However we have also obtained that $$ \mu_{a_1...a_n}=1\cdot\pi_{a_1...a_n}, $$ and so $$ \nabla_\mu \mu=\partial_\mu 1=0, $$ and we obtained the desired result.
Also note that we have only assumed compatibility, but not torsionlessness on the connection $\nabla$.