On two equivalent definitions of divergence

This is my answer to a question that was asked on math.SE, by someone who needed clarification on the limit-integral definition of the divergence operator. There is a lot to talk about, so I thought I would post a complete version of it here.

Let's put it this way. Suppose you have defined the divergence as a differential operator: that is, for a general function $\mathbf f : A \to \mathbb R^n$, where $A$ is a subset of $\mathbb R^n$, and for $\mathbf x_0 \in A \setminus \partial A$ at which $\mathbf f$ is differentiable, let the divergence of $\mathbf f$ at $\mathbf x_0$ be defined as the real number $$\operatorname{div} \mathbf f (\mathbf x_0) := \operatorname{tr} \mathsf D^1\mathbf f(\mathbf x_0) = \sum_{i=1}^n \frac{\partial f_i}{\partial x_i}(\mathbf x_0), \tag1$$ where $\mathsf D^1\mathbf f(\mathbf x_0)$ is the Jacobian matrix of $\mathbf f$ at $\mathbf x_0$ and $\operatorname{tr}$ indicates the trace operator.

Then the following theorem holds:

Theorem 1. Let $\Omega$ be an open subset of $\mathbb R^n$, and let $\mathbf f : \Omega \to \mathbb R^n$ be of class $\mathcal C^1$. Suppose furthermore that $\mathbf x_0 \in \Omega$, and $\{A_k\}_{k \in \mathbb N}$ is a sequence of subsets of $\Omega$ such that

1. For all $k$, $A_k$ is a regular open set (see below);
2. For all $k$, $A_k$ contains the point $\mathbf x_0$;
3. For all $\varepsilon>0$ there is an index $k \in \mathbb N$ such that $\operatorname{diam} A_k < \varepsilon$ — or, equivalently, $$\lim_{k\to\infty} \operatorname{diam}A_k = 0.$$

Then, if $\hat{\mathbf n}_k : \partial A_k \to \mathbb R^n$ is the function associating, to each point of $\partial A_k$, the unit normal vector pointing outward w.r.t. $A_k$, the limit $$\lim_{k\to \infty} \frac{1}{\operatorname{vol}_n A_k } \int_{\partial A_k} \mathbf f \cdot \hat{\mathbf n}_k \mathop{}\!da \tag2$$ exists in $\mathbb R$ and its value is $\operatorname{div} \mathbf f(\mathbf x_0)$.

By $\operatorname{diam}A_k$ we mean the diameter of the set $A_k$, i.e. the greatest possible distance between two points in the closure of $A_k$: $$\operatorname{diam} A_k := \sup_{\mathbf x, \mathbf y \in A_k} |\mathbf x-\mathbf y|;$$ also by $\operatorname{vol}_n$ we mean the $n$-dimensional standard measure defined on $\mathbb R^n$. For $n = 3$ we're dealing with the usual notion of volume.

By regular open set we mean that $A_k$ should have a "nice" border, i.e. that $\partial A_k$ should be a $(n-1)$-dimensional compact manifold of certain regularity (otherwise we will have problems defining the outward unit normal function $\hat{\mathbf n}_k$). More precisely, for an open set $U \subseteq \mathbb R^n$ to be regular, we require that it be bounded, and that for all $\mathbf x_0 \in \partial U$ there exist a neighborhood $I$ of $\mathbf x_0$, a neighborhood $J$ of the origin of $\mathbb R^n$, and a bijection $\mathbf r : J \to I$ of class at least $\mathcal C^1$ such that

1. $\mathbf r^{-1} : I \to J$ is of the same regularity as $\mathbf r$;
2. $I \cap  U = \mathbf r(J \cap \mathbb R^n_+)$;
3. $I \cap \partial U = \mathbf r(J \cap \partial \mathbb R^n_+)$;
4. $I \setminus \overline{U} = \mathbf r(J \setminus \overline{\mathbb R^n_+})$.

Of course, by $\mathbb R^n_+$ we mean the (open) region of $\mathbb R^n$ where $x^n > 0$ (or for that matter any other one of the coordinates in place of $x^n$); this is also called the $n$-dimensional upper half-space.

There is also a useful characterization of these regular open sets:

Proposition 2. A bounded open set $U \subset \mathbb R^n$ is regular if and only if there exists a function $g : \mathbb R^n \to \mathbb R$ that is at least of class $\mathcal C^1$, that verifies $$\begin{split} U &= \{\mathbf x \in \mathbb R^n\ |\ g(\mathbf x) < 0 \}, \\ \partial U &= \{ \mathbf x \in \mathbb R^n\ |\ g(\mathbf x) = 0\}, \\ \mathbb R^n \setminus \overline U &= \{\mathbf x \in \mathbb R^n\ |\ g(\mathbf x) > 0\}, \end{split}$$ and is such that $\nabla g(\mathbf x) \neq \mathbf 0$ for all $\mathbf x \in \partial U$.

Of course, that famous theorem which states that contour lines of regular functions can only be locally parametrized by maps of the same regularity prevents the border $\partial U$ from having any "sharp edges" or "corners" or other similar singular portions. That is why smooth balls are usually employed as the prototypical choice of $A_k$.¹

As an exercise: try to rewrite the statement of the theorem above by replacing the (discrete) sequence $\{A_k\}_{k\in\mathbb N}$ with the (continuous) family of open balls $\{B_{\varepsilon}(\mathbf x_0)\}_{\varepsilon > 0}$. (Spoiler: it becomes much simpler.)

---

Addendum. For completeness, I thought I'd add the proof to Theorem 1 presented above, in the case where $n=3$ and the ever smaller sets are ball-shaped: for $\mathbf x_0 \in \Omega$, we set $$A_k = B_{\varepsilon_k}(\mathbf x_0) = \{ \mathbf x \in \mathbb R^3\ |\ |\mathbf x - \mathbf x_0| < \varepsilon_k \}, \quad \text{with} \quad \varepsilon_k \xrightarrow{k\to\infty} 0^+.$$
Let us rewrite the statement of the theorem for this special case:

Theorem 1a. Let $\Omega$ be an open subset of $\mathbb R^3$, and let $\mathbf f : \Omega \to \mathbb R^3$ be of class $\mathcal C^1$. Suppose furthermore that $\mathbf x_0 \in \Omega$. Then the limit
$$\lim_{\varepsilon \to 0^+} \frac{1}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \int_{\partial B_\varepsilon (\mathbf x_0)} \mathbf f(\mathbf x) \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \mathop{}\!da \tag{2a}$$ exists in $\mathbb R$ and its value is $\operatorname{div}\mathbf f(\mathbf x_0)$.

Proof. As $\mathbf f$ is of class $\mathcal C^1$ over the open set $\Omega$, it is differentiable at all its points (due to a standard theorem in analysis). Because of this, we may surely expand $\mathbf f(\mathbf x)$ to first order: $$\mathbf f(\mathbf x) = \mathbf f(\mathbf x_0) + [\mathsf D^1\mathbf f(\mathbf x_0)](\mathbf x -\mathbf x_0) + |\mathbf x - \mathbf x_0|\mathbf q(\mathbf x - \mathbf x_0), \tag3$$ where $\mathbf q$ is a function such that $$\lim_{\mathbf h \to \mathbf 0} \mathbf q(\mathbf h) = \mathbf 0.$$ Let us set $$D(\varepsilon) := \frac{1}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \int_{\partial B_\varepsilon (\mathbf x_0)} \mathbf f(\mathbf x) \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \mathop{}\!da.$$ and let $\varepsilon > 0$. Substituting our expansion $(3)$ into $(2\text a)$ yields, when $B_\varepsilon(\mathbf x_0) \subseteq \Omega$, $$\begin{split} D(\varepsilon) &= \frac{1}{\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} \mathbf f(\mathbf x_0) \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \mathop{}\!da \\ &\quad + \frac{1}{\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} [\mathsf D^1\mathbf f(\mathbf x_0)](\mathbf x -\mathbf x_0)  \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \mathop{}\!da \\ &\quad + \frac{\varepsilon}{\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} \mathbf q(\mathbf x -\mathbf x_0)  \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \mathop{}\!da, \end{split} \tag4$$ where we have tacitly substituted $|\mathbf x-\mathbf x_0| = \varepsilon$ in the last integral, since it is performed on the surface of the ball. Let us now evaluate each integral, marking each of them as $J_k(\varepsilon)$, $k=1,2,3$.

Notice that $\hat{\mathbf n}_\varepsilon$ is an odd function, i.e. $$\hat{\mathbf n}_\varepsilon(-\mathbf x) = -\hat{\mathbf n}_\varepsilon(\mathbf x), \quad \forall \mathbf x \in \partial B_\varepsilon(\mathbf x_0);$$ hence the integrand $\mathbf f(\mathbf x_0) \cdot \hat{\mathbf n}_\varepsilon$ of $J_1(\varepsilon)$ is also odd, and since the domain $\partial B_\varepsilon(\mathbf x_0)$ is symmetric about the origin, we may conclude that $$J_1(\varepsilon) = 0.$$
Moreover, by Schwarz's inequality, $$\mathbf q(\mathbf x - \mathbf x_0) \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \leq |\mathbf q(\mathbf x - \mathbf x_0)| | \hat{\mathbf n}_\varepsilon(\mathbf x )|,$$ so that by applying a few common properties of integrals we obtain $$\begin{split} |J_3(\varepsilon)| &= \left| \frac{\varepsilon}{\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} \mathbf q(\mathbf x -\mathbf x_0)  \cdot \hat{\mathbf n}_\varepsilon(\mathbf x)\ \mathop{}\!da \right| \\ &\leq \frac{\varepsilon}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \int_{\partial B_\varepsilon(\mathbf x_0)} \Big| \mathbf q(\mathbf x -\mathbf x_0)  \cdot \hat{\mathbf n}_\varepsilon(\mathbf x) \Big| \mathop{}\!da \\ &\leq \frac{\varepsilon}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \int_{\partial B_\varepsilon(\mathbf x_0)} | \mathbf q(\mathbf x -\mathbf x_0)|  | \hat{\mathbf n}_\varepsilon(\mathbf x) | \mathop{}\!da \\ &= \frac{\varepsilon}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \int_{\partial B_\varepsilon(\mathbf x_0)} | \mathbf q(\mathbf x -\mathbf x_0)| \mathop{}\!da \\ &\leq \frac{\varepsilon \operatorname{area} \partial B_\varepsilon(\mathbf x_0)}{\operatorname{vol} B_\varepsilon(\mathbf x_0)} \sup_{\mathbf x \in \partial B_\varepsilon(\mathbf x_0)} | \mathbf q(\mathbf x -\mathbf x_0)|.  \\ \end{split}$$ Therefore, since $$\operatorname{vol} B_\varepsilon(\mathbf x_0) = \frac 4 3 \pi \varepsilon^3, \qquad \operatorname{area} \partial B_\varepsilon(\mathbf x_0) = 4\pi \varepsilon^2,$$ we may conclude that $$|J_3(\varepsilon)| \leq 3  \sup_{\mathbf x \in \partial B_\varepsilon(\mathbf x_0)} | \mathbf q(\mathbf x -\mathbf x_0)|,$$ and since $\mathbf q$ is infinitesimal as $\mathbf x \to \mathbf x_0$, this quantity will disappear as we let $\varepsilon \to 0$.

Turning to $J_2$, we see that for all $\mathbf x \in \partial B_\varepsilon(\mathbf x_0)$, $$\hat{\mathbf n}_\varepsilon(\mathbf x) =  \frac{\mathbf x - \mathbf x_0}{|\mathbf x - \mathbf x_0|} = \frac{\mathbf x - \mathbf x_0}{\varepsilon},$$ so that we may write $$\begin{split} J_2(\varepsilon) &= \frac{1}{\varepsilon\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} [\mathsf D^1\mathbf f(\mathbf x_0)](\mathbf x -\mathbf x_0)  \cdot (\mathbf x - \mathbf x_0) \mathop{}\!da  \\ &= \frac{1}{\varepsilon\operatorname{vol} B_\varepsilon(\mathbf x_0)}  \int_{\partial B_\varepsilon(\mathbf x_0)} (\mathbf x - \mathbf x_0)^\top[\mathsf D^1\mathbf f(\mathbf x_0)](\mathbf x -\mathbf x_0)  \mathop{}\!da.   \end{split}$$ With the substitution $\boldsymbol \xi = (\mathbf x - \mathbf x_0)$, the integrand may be fully expanded as a double sum $$(\mathbf x - \mathbf x_0)^\top[\mathsf D^1\mathbf f(\mathbf x_0)](\mathbf x -\mathbf x_0) = \sum_{i,j=1}^3 \frac{\partial f_i}{\partial x_j}(\mathbf x_0) \xi_i \xi_j,$$ so that $$J_2(\varepsilon) = \frac{1}{\varepsilon\operatorname{vol} B_\varepsilon(\mathbf x_0)} \sum_{i,j=1}^3 \frac{\partial f_i}{\partial x_j}  \int_{\partial B_\varepsilon(\mathbf 0)} \xi_i \xi_j  \mathop{}\!da.\tag5$$

The issue is now reduced to calculating this last integral. First, we notice that the integrand function $\varphi_{ij}$ satisfies this symmetry property: if $\mathbf s_i : \mathbb R^3 \to \mathbb R^3$ is the linear transformation sending $\xi_i \mapsto -\xi$, while $\xi_j \mapsto \xi_j$ for all $j \neq i$, then $$\varphi_{ij}(\mathbf s_i(\boldsymbol \xi)) = -(\xi_i\xi_j) = -\varphi_{ij}(\boldsymbol \xi).$$ Furthermore, the transformation $\mathbf s_i$ is a reflection through the plane normal to the $\xi_i$-axis, so the domain of integration $\partial B_\varepsilon(\mathbf 0)$ is left unchanged; therefore, $$\int_{\partial B_\varepsilon(\mathbf 0)} \xi_i\xi_j \mathop{}\!da = 0, \quad \forall i \neq j.$$ The more interesting case is when $i=j$: the integrals $$\int_{\partial B_\varepsilon(\mathbf 0)} \xi_1^2 \mathop{}\!da, \quad \int_{\partial B_\varepsilon(\mathbf 0)} \xi_2^2 \mathop{}\!da, \quad \int_{\partial B_\varepsilon(\mathbf 0)} \xi_1^3 \mathop{}\!da$$ all have the same value, because the sphere $\partial B_\varepsilon (\mathbf 0)$ is left unchanged by the transformations swapping $\xi_i$ with $\xi_j$ for $i\neq j$, which turn one of the above integrands into another. By the linearity property of integrals, we may write $$\begin{split}\int_{\partial B_\varepsilon(\mathbf 0)} \xi_1^2 \mathop{}\!da + \int_{\partial B_\varepsilon(\mathbf 0)} \xi_2^2 \mathop{}\!da + \int_{\partial B_\varepsilon(\mathbf 0)} \xi_3^2 \mathop{}\!da &= \int_{\partial B_\varepsilon(\mathbf 0)} (\xi_1 + \xi_2 + \xi_3)^2 \mathop{}\!da  \\ &= \int_{\partial B_\varepsilon(\mathbf 0)} \varepsilon^2 \mathop{}\!da = \varepsilon^2(4\pi\varepsilon^2),\end{split}$$ and thus $$\int_{\partial B_\varepsilon(\mathbf 0)} \xi_i^2 \mathop{}\!da = \frac 1 3 \varepsilon^2(4\pi\varepsilon^2) = \varepsilon \frac 4 3 \pi \varepsilon^3 = \varepsilon \cdot \operatorname{vol} B_\varepsilon(\mathbf 0), \quad i = 1,2,3.$$

All in all, $$\int_{\partial B_\varepsilon(\mathbf 0)} \xi_i \xi_j  \mathop{}\!da = \varepsilon \cdot \operatorname{vol} B_\varepsilon(\mathbf 0) \delta_{ij}, \quad \forall i,j = 1,2,3. \tag6$$ From $(6)$ we desume that the double sum in $(5)$ may be reduced to a single sum "over the diagonal": $$\begin{split}J_2(\varepsilon) &= \frac{1}{\varepsilon\operatorname{vol} B_\varepsilon(\mathbf x_0)} \sum_{i,j=1}^3 \frac{\partial f_i}{\partial x_j}  \varepsilon \cdot \operatorname{vol} B_\varepsilon(\mathbf 0) \delta_{ij} \\ &= \frac{1}{\varepsilon\operatorname{vol} B_\varepsilon(\mathbf x_0)} \sum_{i=1}^3 \frac{\partial f_i}{\partial x_i}  \varepsilon \cdot \operatorname{vol} B_\varepsilon(\mathbf 0),  \end{split}$$ and the equality $\operatorname{vol} B_\varepsilon(\mathbf 0)= \operatorname{vol} B_\varepsilon(\mathbf x_0)$ finally yields $$J_2(\varepsilon) = \sum_{i=1}^3 \frac{\partial f_i}{\partial x_i} = \operatorname{tr} \mathsf D^1\mathbf f(\mathbf x_0) = \operatorname{div} \mathbf f(\mathbf x_0).$$

Going back to $(4)$, we see that for our arbitrary choice of $\varepsilon > 0$ we have obtained
$$D(\varepsilon) = 0 + J_2(\varepsilon) + J_3(\varepsilon) =  \operatorname{div} \mathbf f(\mathbf x_0) + J_3(\varepsilon); \tag7$$ since we had established that $J_3(\varepsilon) \xrightarrow{\varepsilon\to 0^+} 0$, taking the limit in $(7)$ yields $$\lim_{\varepsilon \to 0^+} D(\varepsilon) = \operatorname{div} \mathbf f(\mathbf x_0),$$ as we wanted to show. $\blacksquare$

Final note: it is possible to justify the derivation of $(6)$ more precisely in terms of Darboux sums. We will do so in a future post.

---

¹ It is however possible to consider a slightly larger class of subsets $A_k$, for example by allowing $\partial A_k$ to be the union of a finite collection of compact $(n-1)$-manifolds, and the integral over $\partial A_k$ to be calculated as the sum of the integrals over these regular components. This enables us to define the divergence in terms of integrals over shrinking cubes or other polyhedra. ^↩