Rigorous multivariate differentiation of integral with moving boundaries (Leibniz integral rule) The Leibniz integral rule, in its multivariate form, deals with differentiation of the following sort:
$$ \frac{\partial}{\partial t} \int_{D(t)} F({\bf x}, t) \, d{\bf x} \, , \qquad D(t)\in \mathbb{R}^d \, .$$
I am looking for a fully rigorous formulation of this theorem, as well as a proper proof. So far, I could only find:

*

*The one-dimensional case (see e.g., Courant calculus book).

*Physics-flavored proofs, where the normal speed $v_n ({\bf x})$ is not well defined (see e.g., here). While insightful, this is not what I need.

*Differential-geometry and form-based proofs (see the above paper). While these are valid, they are far more general than what I need.

My question: Are you familiar with a multivariate-calculus version and proof of this theorem?
I searched most of the textbooks mentioned here and here, but still no luck. Any ideas? I believe it is not usually covered in graduate curses of advanced multivariate analysis, and so it is pretty difficult to find it in the textbooks.
(This is cross posted from math.se, since a week has passed and it recieved little-to-no-answer. Also, this is not to be confused with the Leibniz rule $(fg)'=f'g +fg'$, as in this post.)
 A: This is probably just a more streamlined version of the proof in the article you linked to, mentioning somewhat less differential forms. I don't know a more elementary proof, but hope this is of use anyway.
I'll write $D_t \subset \mathbb{R}^d$ for the domain of integration and I'll assume that the map $t\mapsto D_t$ is sufficiently regular in the sense that for sufficiently short times, the motion of the domain can be described by the flow of an ambient vector field $X$ on $\mathbb{R}^d$. I.e. for each $t\in \mathbb{R}$ there should exist a vector field $X$ on $\mathbb{R}^d$ with flow $\Phi_h:\mathbb{R}^d \to \mathbb{R}^d$ such that 
$$
D_{t+h}=\Phi_{h}(D_{t}) \qquad \forall h \in (-\epsilon,\epsilon).
$$
(The field $X$ is not unique, one may add any vector field that's tangent to the boundary of $D_t$.)
With such a vector field $X$ given we compute
\begin{align}
\frac{d}{dt} \int_{D_t}F(\mathbf{x},t)d\mathbf{x}
&= \lim_{h\to 0}\frac{1}{h}\left( \int_{D_{t+h}}F(\mathbf{x},t+h)d\mathbf{x} - \int_{D_{t}}F(\mathbf{x},t)d\mathbf{x}\right)\\
&=\lim_{h\to 0}\frac{1}{h}\left( \int_{\Phi_h(D_{t})}F(\mathbf{x},t+h)d\mathbf{x} - \int_{D_{t}}F(\mathbf{x},t)d\mathbf{x}\right)\\
&=\lim_{h\to 0}\frac{1}{h}\left( \int_{D_{t}}\Phi^*_h(F(\mathbf{x},t+h)d\mathbf{x}) - \int_{D_{t}}F(\mathbf{x},t)d\mathbf{x}\right)\\
\end{align}
(The previous equality is just the substitution rule for integrals in several variables. So $\Phi^*_h\left(F(\mathbf{x},t+h)d\mathbf{x}\right)$ denotes the pullback of the differential form $F(\mathbf{x},t+h)d\mathbf{x}$ along $\Phi_h$. It's what you get by replacing every $\mathbf{x}$ with $\Phi_h(\mathbf{x})$ in the expression $F(\mathbf{x},t+h)d\mathbf{x}$.  We continue by swapping limit and  integral)
\begin{align}
&=\int_{D_{t}}\lim_{h\to 0} \frac{\Phi^*_h(F(\mathbf{x},t+h)d\mathbf{x}) - F(\mathbf{x},t)d\mathbf{x}}{h}
\end{align}
(To go on use Taylor $F(\mathbf{x},t+h)=F(\mathbf{x},t)+\frac{\partial F(\mathbf{x},t)}{\partial t}\cdot h+O(h^2)$ as well as linearity of $\Phi^*$ and $\Phi^*_0=\text{id}_{\mathbb{R}^n}$.)
$$
=\int_{D_{t}} \left(\frac{\partial F(\mathbf{x},t)}{\partial t}d\mathbf{x} 
+\lim_{h\to 0}\frac{\Phi^*_h(F(\mathbf{x},t)d\mathbf{x}) - F(\mathbf{x},t)d\mathbf{x}}{h}\right).
$$
Now that limit in the integral is per definition $L_X(F(\mathbf{x},t)d\mathbf{x})$, the Lie derivative of the form 
$\omega=F(\mathbf{x},t)d\mathbf{x}$ with respect to $X$. Using Cartan's magic formula $L_x \omega = d(i_X \omega)+i_X(d\omega)$, which in our case simplifies to $L_x \omega = d(i_X \omega)$ since $F(\mathbf{x},t)d\mathbf{x}$ is of top degree, we get
$$
=\int_{D_{t}}\frac{\partial F(\mathbf{x},t)}{\partial t}d\mathbf{x} 
+\int_{D_{t}} d \left(i_X(F(\mathbf{x},t)d\mathbf{x})\right)
$$
Applying Stokes theorem to the second integral we finally arrive at
$$
\frac{d}{dt} \int_{D_t}F(\mathbf{x},t)d\mathbf{x}
=
\int_{D_{t}}\frac{\partial F(\mathbf{x},t)}{\partial t}d\mathbf{x} 
+\int_{\partial D_{t}} i_X(F(\mathbf{x},t)d\mathbf{x}).
$$
For those who want to get rid of that insertion in the second integral, a direct computation shows that $\int_\Gamma i_X( f(\mathbf{x})d\mathbf{x})=\int_\Gamma f(\mathbf{x}) X\cdot d\vec{A}$ where $d\vec{A}$ denotes the usual "surface element" with respect to the standard metric. So you can write the Leibniz integral rule as 
$$
\frac{d}{dt} \int_{D_t}F(\mathbf{x},t)d\mathbf{x}
=
\int_{D_{t}}\frac{\partial F(\mathbf{x},t)}{\partial t}d\mathbf{x} 
+\int_{\partial D_{t}} F(\mathbf{x},t)\, X\cdot d\vec{A}.
$$
