I was reading this document on how to get some common operators when dealing with general orthogonal curvilinear coordinates. I am interested in particular in equation (12). It basically defines three orthogonal coordinates $u,v,w$ that depend on $x,y,z$ (the cartesian coordinates) with scale factors $h_u, h_v, h_w$. The document asks to consider a "curvilinear" box (a very small one) and computes the flux through that box to find an expression for the divergence. This computation is divided in the directions of $u,v,w$. In the $u$ direction the computation yields that that the flux of the field $\boldsymbol{B}$ is $$(B_uh_vh_w)_{u+du/2}dvdw - (B_uh_vh_w)_{u-du/2}dvdw.$$
Next they use Taylor Series and they only take the first derivative into account ($\dfrac{\partial(B_uh_vh_w)}{\partial u}dudvdw$). Moreover they equate this term to the flux, as if the other terms of the series were $0$. If the coordinates are for instance the spherical coordinates, then equation (15) of the document: $$Div\boldsymbol{B} = \frac{1}{h_rh_\theta h_\phi}\left[\frac{\partial B_r h_\theta h_\phi}{\partial r} + \frac{\partial B_\theta h_r h_\phi}{\partial \theta} + \frac{\partial B_\phi h_\theta h_r}{\partial \phi}\right]$$ is not an approximation (as one would expect because the superior order terms were ignored), but it is exactly the same expression as the expression that one gets when one writes $Div\boldsymbol B$ in cartesian coordinates and uses the transformation formulas to write $Div\boldsymbol B$ in spherical coordinates.
Another situation in which this happens is in infinitesimal rotations of coordinate systems, where second and superior order terms are ignored (see Goldstein Mechanics, Ch 4). I won't talk about this though...
My question is: Why are second and superior order terms of taylor series ignored when one is expanding a differential (or something infinitesimal)?
$\newcommand{\dv}{\mathop{\rm div}}$ From the divergence theorem $$\int_\Omega \dv B = \int_{\partial \Omega} B \cdot n$$ we see that if we want to define the divergence via flux, the correct quantity is (assuming the divergence is continuous)
$$ \dv B (x) = \lim_{r \to 0} \frac1{V(\Omega(x,r))}\int_{\Omega(x,r)} \dv B =\lim_{r \to 0} \frac1{V(\Omega(x,r))}\int_{\partial \Omega(x,r)}B \cdot n$$
where $\Omega(x,r)$ is a "box" around the point $x$ with side lengths scaled by $r$ and $V$ denotes volume. Since the volume of a 3-dimensional box $B(x,r)$ is proportional to $r^3$, the omitted terms in the Taylor series do not contribute to the divergence: since they will on the order of $r^4$ or higher, when we divide by $r^3$ they will vanish in the limit $r \to 0$.