I'm going to explain my thought process. Question and counter-arguments at the end.
Divergence Theorem:
$$\int_{\mathscr{D}_i\subset \mathbb{R}^3} \nabla \cdot \vec{f} dV_i=\int_{\partial\mathscr{D}_i\subset \mathbb{R}^2} \vec{f} \cdot \hat{n}_i dS_i,\quad i\in\mathbb{N}\quad (1)$$
The volume integral over some domain, ${\mathscr{D}}$, of the divergence of a vector, $\vec{f}$, equals the surface integral of the inner product of, $\vec{f}$, and the unit vector, $\hat{n}$, normal to the boundary of that volume (denoted ${\partial\mathscr{D}}$). The reason I call it the $i^{th}$ domain will be clear in a moment.
Vanishing Theorem:
$$\int_{\forall\mathscr{D}_i\subset \mathbb{R}^n} g \ dx= 0 \rightarrow g\equiv0,\quad i\in\mathbb{N} \quad (2)$$
If the integral of some function, g, over some n-dimensional domain is zero for all domains in $\mathbb{R}^n$, then g can only be zero.
We are generally taught to use (1) on a single volume at a time. If we were interested in two non-intersecting volumes, we would simply solve the system:
$$\int_{\mathscr{D}_1\subset \mathbb{R}^3} \nabla \cdot \vec{f} dV_1=\int_{\partial\mathscr{D}_1\subset \mathbb{R}^2} \vec{f} \cdot \hat{n}_1 dS_1\quad (3a)$$
$$\int_{\mathscr{D}_2\subset \mathbb{R}^3} \nabla \cdot \vec{f} dV_2=\int_{\partial\mathscr{D}_2\subset \mathbb{R}^2} \vec{f} \cdot \hat{n}_2 dS_2\quad (3b)$$
Suppose $\vec{f}$ were divergence free ($\nabla \cdot \vec{f}=0$). (3a-b) are now,
$$\int_{\partial\mathscr{D}_1\subset \mathbb{R}^2} \vec{f} \cdot \hat{n}_1 dS_1=\int_{\partial\mathscr{D}_2\subset \mathbb{R}^2} \vec{f} \cdot \hat{n}_2 dS_2=0 \quad (4)$$
The nature of the solution depends on many things. I want to draw your attention to $\hat{n}_1$, and $\hat{n}_2$. Above all else, their relative orientation determines the existence and uniqueness of the solution… Now onto my question… Suppose I have an equation like the continuity equation from physics. A generic form of this equation could be:
$$\int_{\forall\mathscr{D}_i\subset \mathbb{R}^3} a+\nabla \cdot \vec{b} \ dV= 0\quad (5)$$
Normally authors don't write "$\forall\mathscr{D}_i$". Rather they say that the equation holds for an "arbitrary volume." To me, the concepts are equivalent but I welcome pedantic feedback if you see room for error. If $\nabla \cdot \vec{b}=0$, then by linearity of integration I can simply write,
$$\int_{\forall\mathscr{D}_i\subset \mathbb{R}^3} \nabla \cdot \vec{b} \ dV_i= 0\quad (6)$$
By (1) it then follows that,
$$\int_{\forall\partial\mathscr{D}_i\subset \mathbb{R}^2} \vec{b} \cdot \hat{n}_i dS_i= 0\quad (7)$$
Here's where my problems start. By (2) it should follow that,
$$\int_{\forall\partial\mathscr{D}_i\subset \mathbb{R}^2} \vec{b} \cdot \hat{n}_i dS_i= 0 \, \rightarrow \vec{b}\cdot \hat{n}_i=0 \, \forall \; i\in\mathbb{N} \quad (8)$$
Since "$\forall\partial\mathscr{D}_i$" is an infinitely large set, then from (8) I should be able to conclude that,
$$\vec{b}\cdot \hat{n}_i=0\rightarrow\vec{b}\equiv 0\quad (9)$$
Thus if (2) holds, the only divergence-free vector permitted by (5) is the zero vector. Question: Is (9) a valid special case, or a fallacy? I presented this to two math PhD's and got the following counter arguments:
Counter #1: If $\vec{b}$ is constant it is divergence free. My response to that: If $\vec{b}$ is constant, it cannot satisfy (7) and (7) should follow directly from (1), therefore, a constant non-zero vector should not be a valid solution even though it's divergence free.
Counter #2: The divergence theorem (1) does not give us "new" information, therefore, we cannot infer anything about a vector given only the fact that it is divergence free. My response to that: Since (1) is for a single volume, then I agree that for a single volume the inference is impossible. However, I claim that the "new information" comes from (2). By combining (1) and (2) I am saying (1) must hold over all possible volumes, and that "new" information should make (9) valid.
A quick note: all feedback is appreciated but I respectfully request that the responses be pure math responses. Responses based on physics would broaden the scope of the discussion too much. I want to get to the heart of the mathematics here.
Best Answer
Step 8 is incorrect. You know that for any appropriate domain $\mathscr{D}$, $\int_{\partial\mathscr{D}}\vec{b} \cdot \hat{n} dS=0$, where $\hat{n}$ is the normal vector to $\partial\mathscr{D}$. But you cannot conclude from this that the function $\vec{b}\cdot\hat{n}$ is identically zero, because $\hat{n}$ is different for each different choice of $\mathscr{D}$. There is no single function $g$ such that $g=\vec{b}\cdot\hat{n}$ and the integral of $g$ on any domain is $0$. Rather, for each domain, you have a function $\vec{b}\cdot\hat{n}$ depending on that domain (since $\hat{n}$ depends on $\mathscr{D}$) whose integral is $0$. You can't deduce from this that any function vanishes identically.
As your "Counter #1" stated, the conclusion you reach in (9) is in fact not true, since a constant vector field $\vec{b}$ is a counterexample. A constant vector field actually does satisfy (7) (and this is no problem since your step 8 does not follow).
(There are also some more minor issues. For instance, it is not correct to write that $\partial\mathscr{D}\subset\mathbb{R}^2$: the set $\partial\mathscr{D}$ is a surface in $\mathbb{R}^3$, not an open subset of the plane. So you can't apply the vanishing theorem in the form you stated: instead you would need a version of it where you have a function on $\mathbb{R}^3$ and you are integrating it over all possible closed surfaces you can get as boundaries. But this is not an essential point--the key error in your reasoning is in Step 8.)