Transformation of Laplace Operator for arbitrary Diffeomorphism (using Matrix Calculus)

For a diffeomorphism $T\colon \mathbb{R}^n \rightarrow \mathbb{R}^n$ and some function $\phi \in C^2(\mathbb{R}^n)$ I am trying to prove $$\Delta \phi \circ T = \frac{1}{\gamma} \operatorname{div}\left( A \nabla(\phi \circ T) \right)$$ where $\gamma\colon \mathbb{R}^n \rightarrow \mathbb{R}, x \mapsto \det(DT(x))$ and $A\colon \mathbb{R}^n \rightarrow \mathbb{R}, x\mapsto \gamma(x) DT(x)^{-1} DT(x)^{-T}$. However, I am kind of stuck in my proof.

For my proof I tried to avoid the usage of $D^2 T \in \mathbb{R}^{n\times n\times n}$ (actually I did use it in another version but at some steps I felt like it is too unclear what is happening, so I instead chose a more "elementary" approach). So far, I did the following:

Because of $\nabla (\phi \circ T) = DT^{T} (\nabla \phi\circ T)$ we have \begin{align*} A\nabla(\phi\circ T) & = \gamma DT^{-1} (\nabla \phi\circ T) \\ & = \gamma \left( \sum_{j=1}^n DT^{-1}_{ij} (\partial_j\phi \circ T)\right)_{i\in\{1,\dots,n\}}\text{.} \end{align*} Consequently, \begin{align*} \operatorname{div}\left( A \nabla(\phi \circ T) \right) & = \operatorname{div} \left( \gamma \left( \sum_{j=1}^n DT^{-1}_{ij} (\partial_j\phi \circ T) \right)_{i\in\{1,\dots,n\}} \right) \\ & = \sum_{i=1}^n \partial_i \left( \gamma \sum_{j=1}^n DT^{-1}_{ij} (\partial_j\phi \circ T) \right) \\ & = \sum_{i,j=1}^n (\partial_i \gamma) DT^{-1}_{ij} (\partial_j \phi \circ T) + \sum_{i,j=1}^n \gamma \partial_i \left( DT^{-1}_{ij} (\partial_j \phi \circ T) \right)\text{.} \end{align*} For $\partial_i \gamma$ we consider the formula for the derivative of the determinant and get \begin{align*} \partial_i \gamma & = \gamma \operatorname{Trace}(DT^{-1} \partial_i DT) \\ & = \sum_{l,m=1}^n \gamma DT^{-1}_{lm} \partial_i DT_{ml}\text{.} \end{align*} For $\partial_i \left( DT^{-1}_{ij} (\partial_j \phi \circ T) \right)$ we have \begin{align*} \partial_i \left( DT^{-1}_{ij} (\partial_j \phi \circ T) \right) & = (\partial_i DT^{-1}_{ij}) (\partial_j \phi \circ T) + \sum_{l=1}^n DT^{-1}_{ij} (\partial_{jl} \phi \circ T) \partial_i T_l \\ & = - \sum_{l,m=1}^n DT^{-1}_{il} \partial_i DT_{lm} DT^{-1}_{mj} (\partial_j \phi \circ T) + \sum_{l=1}^n DT^{-1}_{ij} (\partial_{jl} \phi \circ T) DT_{li} \end{align*} where we used the formula for the derivative of inverse matrices, here given by \begin{align*} \partial_i DT^{-1} = - DT^{-1} (\partial_i DT) DT^{-1}\text{.} \end{align*} Putting everything together, we obtain \begin{align*} \operatorname{div}\left( A \nabla(\phi \circ T) \right) & = \sum_{i,j,l,m=1}^n \gamma DT^{-1}_{lm} \partial_i DT_{ml} DT^{-1}_{ij} (\partial_j \phi \circ T) \\ & \quad - \sum_{i,j,l,m=1}^n \gamma DT^{-1}_{il} \partial_i DT_{lm} DT^{-1}_{mj} (\partial_j \phi \circ T) \\ & \quad + \sum_{i,j,l=1}^n \gamma DT^{-1}_{ij} (\partial_{jl} \phi \circ T) DT_{li} \\ & = \sum_{i,j,l,m=1}^n \gamma DT^{-1}_{ml} \partial_i DT_{lm} DT^{-1}_{ij} (\partial_j \phi \circ T) \\ & \quad - \sum_{i,j,l,m=1}^n \gamma DT^{-1}_{il} \partial_i DT_{lm} DT^{-1}_{mj} (\partial_j \phi \circ T) \\ & \quad + \gamma\Delta \phi \circ T \end{align*} where we used \begin{align*} \sum_{i,j,l=1}^n DT^{-1}_{ij} \partial_i DT_{lm} DT^{-1}_{mj} (\partial_j \phi \circ T) = \operatorname{Trace}(DT^{-1} (D^2\phi\circ T) DT) = \operatorname{Trace}(D^2\phi\circ T) = \Delta \phi \circ T\text{.} \end{align*} At this point I am stuck because I don't see why the terms containing the $\partial_i DT_{lm}$ would cancel out. Do you see how to solve this or where I maybe did apply indices in the wrong order? Or do you have a simpler proof?

Remark: In another attempt where I solely used matrix notations (so no index battles) I got stuck with showing \begin{align*} \operatorname{Trace} \left( \operatorname{Trace}\left(DT^{-1} D^2T\right)DT^{-1} D\phi\circ T\right) =\operatorname{Trace}\left( DT^{-1} D^2T DT^{-1} D\phi\circ T\right)\text{,} \end{align*} but also here I do not see why this would hold.

Result: Let $\Omega_0, \Omega_1 \subseteq \mathbb{R}^n$. Furthermore, let $T\colon \Omega_0 \rightarrow \Omega_1$ be a diffeomorphism with the property that \begin{align*} H_0^2(\Omega_0) & \longrightarrow H_0^2(\Omega_1) \\ v & \longmapsto v \circ T \end{align*} is an isomorphism. Then we have for $h\in H_0^2(\Omega)$ that \begin{align*} \Delta h \circ T = \frac{1}{\gamma} \operatorname{div}( A(t) \nabla(h\circ T) ) \end{align*}
Proof: (Some parts are only sketched. Feel free to tell me if something is unclear or incorrect. Some parts will definitely be incorrect because there are a few more assumptions that I have and that are implicitly assumed in the proof, like $\gamma >0$.)
Let $h \in H_0^2(\Omega)$. We have for all $\phi \in H^1(\Omega)$ \begin{align*} \int_{\Omega_1} \Delta h \phi & = -\int_{\Omega_1} \nabla h \cdot \nabla \phi + \int_{\partial\Omega_1} \partial_n h \phi \\ & = -\int_{\Omega_1} \nabla h \cdot \nabla \phi \\ & = -\int_{\Omega_0} \nabla(\phi \circ T) \cdot \left(\gamma DT^{-1} DT^{-T} \nabla (h\circ T) \right) \\ & = \int_{\Omega_0} (\phi \circ T) \operatorname{div}( A \nabla(h\circ T) ) - \int_{\Omega_0} \operatorname{div}( (\phi \circ T)A \nabla(h\circ T)) \\ & = \int_{\Omega_0} (\phi \circ T) \operatorname{div}( A \nabla(h\circ T) ) \end{align*} because $h \circ T \in H_0^2(\Omega_0)$ and thus \begin{align*} \int_{\Omega_0} \operatorname{div}((\phi \circ T)A \nabla(h\circ T)) = \int_{\partial\Omega_0} (\phi \circ T) A \nabla(h\circ T) \cdot n = 0 \end{align*} in the line second to the last. On the other hand, we have \begin{align*} \int_{\Omega_1} \Delta h \phi & = \int_{\Omega_0} (\Delta h \circ T) (\phi \circ T) \gamma\text{.} \end{align*} By using an appropriate version of the fundamental theorem of variational calculus we get the assertion.