In the context of a 4D (pseudo-) Riemannian manifold, I've been trying to show that the Lie derivative of the volume form $\Omega =\sqrt{-g}\,dx^{1}\wedge dx^{2}\wedge dx^{3}\wedge dx^{4}:=\sqrt{-g}\,d^{4}x$ is $\mathcal{L}_{X}\Omega=\nabla_{\mu}X^{\mu}\Omega$, however, I'm a bit confused about how $\sqrt{-g}$ behaves under lie and exterior differentiation.
As I understand it, $\sqrt{-g}$ is technically a scalar density, i.e. it transforms as $\sqrt{-g}\rightarrow \lvert J\rvert^{-1}\sqrt{-g}$ under coordinate transformations (where $J$ is the Jacobian determinant of the coordinate transformation). This leaves me confused as to whether $\sqrt{-g}$ behaves as a scalar or something different under Lie or exterior differentiation.
What has further confused me, is that if I treat $\sqrt{-g}$ as a scalar, then I can get the Lie derivative of it in two different ways and end up with different results. The first is by brute force: $$\mathcal{L}_{X}\sqrt{-g}=\frac{\partial\sqrt{-g}}{\partial g_{\mu\nu}}\,\mathcal{L}_{X}g_{\mu\nu}=\sqrt{-g}\,g^{\mu\nu}\nabla_{\mu}X_{\nu}=\sqrt{-g}\,\nabla_{\mu}X^{\mu}$$ and the second using Cartan's identity: $$\mathcal{L}_{X}\sqrt{-g}=\iota_{X}\mathrm{d}(\sqrt{-g})=X^{\mu}\partial_{\mu}(\sqrt{-g})=\frac{1}{2}X^{\mu}\sqrt{-g}g^{\lambda\sigma}\partial_{\mu}g_{\lambda\sigma}=\sqrt{-g}\,X^{\mu}\Gamma^{\nu}_{\mu\nu}$$ I know I must be doing something wrong, but I can't see the wood for the trees. If someone could explain this for me I'd much appreciate it.
You can think about $\sqrt{|g|}$ as a density if you wish, but I don't think it will help you with this calculation. To keep everything concrete, let's treat it as a scalar function, which requires us to fix a single coordinate system $x^i$.
Since the Lie derivative of a scalar is just the ordinary (directional) derivative, your "brute force" calculation should start off with $\mathcal L_X \sqrt{-g} = X^\mu \partial_\mu \sqrt{-g},$ at which point it is clearly the same as the expression you get from Cartan's identity.
You instead seem to be assuming some kind of chain rule for Lie derivatives, which (as this example shows) is incorrect. I think the main point of confusion here is the distinction between the metric tensor itself and its coordinate components - it is the latter which $\sqrt{- g}$ is written as a function of, since we really mean the determinant of the matrix of $g$ in our fixed coordinates $x^i$. Your formula is correct if you change $\mathcal L_X g_{\mu \nu}$ (which we interpret as $(\mathcal L_X g)_{\mu \nu}$, i.e. the Lie derivative of a 2-tensor) to $\mathcal L_X(g_{\mu \nu})$, which is just the Lie derivative of the coordinate components; but of course this is just $X^\alpha \partial_\alpha (g_{\mu\nu}).$