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$.
Best Answer
I am not sure if this is what you are looking for, but nevertheless, let me say the following:
In full generality, a "vielbein" on a d-dimensional spacetime is defined to be a vector bundle isomorphism
$$e:T\mathcal{M}\to\mathcal{T},$$
where $T\mathcal{M}=\coprod_{p\in\mathcal{M}}T_{p}\mathcal{M}$ denotes the tangent bundle of the spacetime manifold $(\mathcal{M},g)$ and where $\mathcal{T}$ is some other vector bundle, called "internal space". Note that equivalently one can view $e$ as a $\mathcal{T}$-valued $1$-form, i.e. $e\in\Omega^{1}(\mathcal{M},\mathcal{T})$. Now, the bundle $\mathcal{T}$ is usually chosen to be the associated vector bundle $\mathcal{T}:=F_{\mathcal{Ort}}(\mathcal{M})\times_{\rho}\mathbb{M}^{d}$, where $F_{\mathcal{Ort}}(\mathcal{M})$ is the "orthonormal frame bundle", which is a principal $\mathrm{SO}(1,d-1)$-bundle, whose fibres are just the orthonormal bases of the tangent spaces $T_{p}\mathcal{M}$, where $\mathbb{M}^{d}=\left(\mathbb{R}^{d},\eta\right)$ is the Minkowski space and where $\rho:\mathrm{SO}(1,d-1)\to\mathrm{Aut}(\mathbb{M}^{d})$ is the fundamental representation.
Now, if the manifold $\mathcal{M}$ is parallelizable, then you can identify $e\in\Omega^{1}(\mathcal{M},\mathcal{T})$ with an element $e\in\Omega^{1}(\mathcal{M},\mathbb{M}^{d})$ and hence, you can write that $e=e^{a}v_{a}$, where $v_{a}$ is an orthonormal basis of $\mathbb{M}^{d}$ and where $e^{a}$ are real-valued coordinate forms. Furthermore, you can write the real-valued forms $e^{a}$ as $e^{a}=e^{a}_{\mu}\mathrm{d}x^{\mu}$, as usual. Using orthonormality, we then have that $\eta_{ab}e^{a}_{\mu}e^{b}_{\nu}=g_{\mu\nu}$, which in physics is often used as definition of the vielbein fields.
Of course, strictly speaking, the map $e$ discussed so far are the "co-vielbein-fields" and its inverse $e^{\mu}_{a}$ are the "vielbein fields".
This answer is mathematically quite technical, but since you asked for a mathematical definition, I thought I just write it here. If something is not clear, then do not hesitate to write a comment below. :-)
Answers to the additional questions of OP's Edit:
Answer to the question in the comments:
Yes, the vector bundle $F_{\mathrm{Ort}}(\mathcal{M})\times_{\rho}\mathbb{M}^{d}$ is defined using an equivalence relation. This is a particular example of so-called "associated vector bundles". Consider an arbitrary principal $G$-bundle $P$ over $\mathcal{M}$, which is a fibre bundle $G\to P\to\mathcal{M}$ with some extra properties. Now, take a finite-dimensional real or complex representation $(\rho,V)$ of $G$. Then one defines a vector bundle denoted by $P\times_{\rho}V=:E$ to be the vector bundle with fibres given by $E_{x}:=(P_{x}\times V)/G$, where the equivalence relation is defined by $$(p,v)\sim (p\cdot g,\rho(g)^{-1}v).$$ The symbol "$\cdot$" here denotes the group action of the principal bundle. The vector space structure of the fibres $E_{x}$ is the obvious one, i.e. $$[p,v]+\lambda [p,w]:=[p,v+\lambda w].$$ Let me stress that the notion of associated vector bundles is also very important in physics. You might know that "fields" are in mathematical physics usual defined to be sections of bundles. Now, if you have a gauge theory described by some principal bundle, then matter fields coupled to gauge fields (like in Klein Gordon theory with gauge fields, Higgs field, standard model,...) are sections of the associated vector bundles of the principal bundle with the representation in which the transform. (At least for the bosonic case. For fermions, you need the notion of twisted spinor bundles, which is more complicated.)
Literature:
I am also not an expert on vielbein fields, but I needed the mathematical definition myself some time ago. Hence, I know that it is quite hard to find a mathematical discussion in the literature, since most of the ressources just define them locally. However, after some time I found some ressources, which I can share here with you. These references are not precisely about vielbein fields in general, but include some nice discussion of their mathematical definition and structure: