Let $\Omega \subset \mathbb{R}^m$ be some compactly contained domain in $\mathbb{R}^m.$ Let $Q_t(x) = x+t\zeta$ where $\zeta\in C^{\infty}_c(\Omega, \mathbb{R}^m)$ and $t$ is small enough so that $Q_t$ is a diffeomorphism on $\Omega.$ My goal is to compute, $$\frac{d}{dt}\left(\int_{\Omega} |D(u(x+t\zeta))|^2dx\right)_{|t=0} = \int_{\Omega}\left(\frac{|Du|^2}{2}\operatorname{div}(\zeta) - \partial_{\alpha}u^i \partial_{\beta} u^i \partial_{\beta}\zeta^\alpha\right) dx$$ where $u\in H^{1}(\Omega, \mathbb{R}^n),$ and we implicitly sum over indices $1\leq \alpha, \beta \leq m$ and $1\leq i \leq n.$
I am reading a proof where the author considers the map, $u_t(x)=u(Q_t^{-1}(x))$ and argues that, $$\int_{\Omega} |D(u(x+t\zeta))|^2dx =\int_{\Omega} |Du_t(x)|^2dx.$$ I am not sure how to prove this, I tried change of variable, by setting $y=Q_t(x)$ in the first integral, to get $x=Q_t^{-1}(y)\implies dx = |\det DQ_t^{-1}(y)| dy$ and so, $$\int_{\Omega} |D(u(x+t\zeta))|^2dx = \int_{\Omega} |Du(x+t\zeta) DQ_t(x)|^2 dx \\ =\int_{\Omega} |Du(y)|^2 |DQ_t(Q_t^{-1}(y))|^2 |\det DQ_t^{-1}(y)| dy\\ =\int_{\Omega} |Du(y)|^2 |DQ_t(Q_t^{-1}(y))|^2 |\det DQ_t^{-1}(y)| dy.$$
On the other hand, the author claims that, $$\int_{\Omega} |Du_t(x)|^2dx = \int_{\Omega} |Du(x)|^2 DQ_t^{-1}(Q_t(x))|^2 |\det DQ_t(x)| dx.$$
How do I show that the two expressions are the same?
The equality $$\int_{\Omega} |D(u(x+t\zeta))|^2\,\mathrm{d}x =\int_{\Omega} |Du_t(x)|^2\,\mathrm{d}x$$ follows immediately from the equality $D(u(x+t\zeta)) = Du_t(x)$ by virtue that $u_t(x) = u(x+t\zeta(x)).$ This is largely a matter of notation, and there is no change of variable involved (and this is not actually needed to show your desired identity).
Your change of variables $y = Q_t(x)$ is also not quite correct, as you should have $$ D(u(x+t\zeta(x)) = DQ_t(x)Du(Q_t(x)) = Du(Q_t(x)) + tD\zeta^i(x)D_iu(Q_t(x)) $$ using the same summation convention. Here we see that $DQ_t(x)$ is a linear operator acting on $Du(Q_t(x)),$ so $$|DQ_t(x)Du(Q_t(x))|^2 \neq |DQ_t(x)|^2|Du(Q_t(x))|^2$$ in general. Note that it can be instructive at this point to write everything in coordinates, as it may clarify what's going on.
Hence a change of variables should give $$ \int_{\Omega} |D(u(x+t\zeta))|^2\,\mathrm{d}x =\int_{\Omega} |DQ_t(Q_t^{-1}(y))Du(y)|^2 |\det DQ_t^{-1}(y)| \,\mathrm{d}y \\ =\int_{\Omega} |Du(y) + tD\zeta^i(Q_t^{-1}(y))D_iu(y)|^2 |\det DQ_t^{-1}(y)| \,\mathrm{d}y.$$ This is not quite what you write at the end, and I don't believe the equality you quote holds true.
However the above does allow you to derive the identity you seek to prove. Note that $y$ is just a integration variable now, so you are free to write it as $y$ or $x;$ in particular it is independent of $t$ (which is important in what follows). Now observe that we have $$ \frac{\mathrm{d}}{\mathrm{d}t}\left.\left( tD\zeta^i(Q_t^{-1}(x))D_iu(x)\right)\right|_{t=0} = D\zeta^i(x)D_iu(x),$$ since if the derivative lands on the $\zeta$ term what you get vanishes as you set $t=0.$ Hence differentiating under the integral sign we have $$ \frac{\mathrm{d}}{\mathrm{d}t}\left.\left(\int_{\Omega} |D(u(x+t\zeta))|^2\,\mathrm{d}x\right)\right|_{t=0}=\int_{\Omega} \frac{\mathrm{d}}{\mathrm{d}t}\left.\left(|Du(x) + tD\zeta^i(Q_t^{-1}(x))D_iu(x))|^2\right)\right|_{t=0} \,\mathrm{d}x + \int_{\Omega} |Du(x)|^2\frac{\mathrm{d}}{\mathrm{d}t}\left.\left(|\det DQ_t^{-1}(x)|\right)\right|_{t=0} \,\mathrm{d}x \\ = \int_{\Omega} 2 D_ju^kD_j\zeta^iD_iu^k + |Du|^2 \operatorname{div} \zeta \,\mathrm{d}x,$$ as required (up to a factor of $\frac12$). I've used Jacobi's formula to differentiate the second term, which I presume your text also explains.
Also for reference to anyone else reading this, the above computes the inner variations for harmonic maps in Euclidean space.