EDIT: I believe I have found the right answer, but would like confirmation. In the application of the mean value theorem, I overlooked that $f(\cdots)-f(\cdots)$ can be considered as one continuous function, whence the mean value should be the same for both; $\vartheta_{z_1}=\vartheta_{z_2}$, at which point my problem is solved trivially. I leave the workings below for posterity's sake, with the corrections in red, but as mentioned I'd like confirmation that this is rigorous. The issue of the existence of the limits still stands!
Op:
$\newcommand{\div}{\operatorname{div}}\newcommand{\curl}{\operatorname{curl}}\newcommand{\d}{\mathrm{d}}$I begin with the following definitions, where $f:\Bbb R^3\to\Bbb R^3$ is a vector field, in the case of curl, and where $f:\Bbb R^n\to\Bbb R^n$ is a generic vector field in $n\ge2$ dimensions in the case of divergence. In both cases $f$ is assumed to be $C^1$.
$$\div_f(p_0)=\lim_{\mu(V)\to0}\frac{1}{\mu(V)}\oint_{\partial V}f\cdot\vec n\,\d S$$ Where $V$ is chosen from the family of all closed neighbourhoods of $p_0$, $\mu$ is Lebesgue measure, and $\vec n$ is the vector normal to (and outward from) $V$ at each point on its surface.
$$\curl_f(x_0,y_0,z_0)\cdot\vec n=\lim_{\mu(A)\to0}\frac{1}{\mu(A)}\oint_{\partial A}f\cdot\d\vec r$$ Where $A$ is any area normal everywhere to $\vec n$, and containing $(x_0,y_0,z_0)$, and its boundary is oriented according to the right hand rule.
I feel like these definitions could use some work, namely stipulations on the differentiability and measurability of the boundaries and surfaces, but this was what I was able to recover from Wikipedia. Certainly these are intuitively the divergence and the curl, but I want to go from these definitions and arrive at the familiar $\nabla\cdot$, $\nabla\times$ formulae.
In both cases I am concerned with whether or not the limits actually exist, assuming $f$ is continuously differentiable. I believe Lebesgue's differentiation theorem can help us out with the existence of the divergence limit, but as $V$ is not necessarily a ball I am not sure. As for the curl, it is not at all clear to me that the limit should always exist, even for $C^1$ vector fields. This is the first sub-question that I have. This is now my only question, really.
Assuming now that the limits do exist, then any family of areas/volumes which tend to zero in Lebesgue measure will suffice. I then will choose cubes in $\Bbb R^n$ for divergence, and squares for curl.
Beginning with curl:
We can construct $\vec R=\curl_f(x_0,y_0,z_0)$ from its components in the direction of $\vec n=\vec i,\vec j,\vec k$. As all of these have the same derivation, I will focus arbitrarily on $\vec n=\vec i$, the $x$-component of the curl. Choose $h\gt0$, and $A$ to be the square with vertices: $$\{(x_0,y_0-h,z_0-h),(x_0,y_0+h,z_0-h),(x_0,y_0+h,z_0+h),(x_0,y_0-h,z_0+h)\}$$And I choose $\partial A$ to be oriented in precisely that order. The boundary has four segments where only the $y$ or $z$ coordinate varies, and thus $f\cdot\d\vec r$ will be the $y$ or $z$ component of $f$, denoted $f_y,f_z$, and so $\vec R\cdot\vec i$ can be found as the limit: $$\lim_{h\to0}\frac{1}{4h^2}\int_{-h}^h\left\{f_y(x_0,y_0+t,z_0-h)-f_y(x_0,y_0+t,z_0+h)\\\hspace{80pt}+f_z(x_0,y_0+h,z_0+t)-f_z(x_0,y_0-h,z_0+t)\right\}\,\d t$$We want to show that this results in the familiar $\partial_y f_z-\partial_z f_y$. The result follows if I can show: $$\lim_{h\to0}\frac{1}{4h^2}\int_{-h}^hf_z(x_0,y_0+h,z_0+t)-f_z(x_0,y_0-h,z_0+t)\,\d t=\frac{\partial f_z}{\partial y}\biggr |_{(x_0,y_0,z_0)}$$Since the case for $f_y$ is the same with the signs reversed.
As $f$ is assumed continuously differentiable, so too is the component $f_z$. $f_z(x_0,y_0+h,z_0+t)$ can be considered a $C^1$ function from $\Bbb R\to\Bbb R$, in the variable $t$, and I can employ the mean value theorem for integration, applied similarly to the other term, to find: $$\lim_{h\to0}\frac{1}{4h^2}[(2h)f_z(x_0,y_0+h,z_0+\vartheta_{z_1})-(2h)f_z(x_0,y_0-h,z_0+\vartheta_{z_2}]\\\hspace{10pt}=\lim_{h\to0}\frac{f_z(x_0,y_0+h,z_0+\vartheta_{z_1})-f_z(x_0,y_0-h,z_0+\vartheta_{z_2})}{2h}$$
For $\vartheta_{z_1,z_2}\in(-h,h)$. By differentiability, for any $\epsilon\gt0$ I can select $h$ so that: $$|f_z(x_0,y_0+h,z_0+\vartheta_{z_1})-\partial_yf_z(x_0,y_0,z_0)\cdot h-\partial_z f_z(x_0,y_0,z_0)\cdot\vartheta_{z_1}|\lt\epsilon/2$$And it follows, by doing the same for the other term, that the expression inside the limit is within $\epsilon$ of: $$\partial_yf_z(x_0,y_0,z_0)+\partial_zf_z(x_0,y_0,z_0)\cdot\frac{\vartheta_{z_1}-\vartheta_{z_2}}{2h}\color{red}{=\partial_yf_z(x_0,y_0,z_0)\,\because\vartheta_{z_1}=\vartheta_{z_2}}$$
Since $-h\lt\vartheta_{z_1,z_2}\lt h$, the quotient on the right is less than $1$ but I do not know how to show it tends to zero. I am almost there but not quite!
I have similar formal concerns with the divergence:
Assuming existence of the limit, I can pick volumes $V$ which are cubes with faces oriented to the coordinate axes. A very similar derivation occurs, but is harder to write down as it occurs in arbitrary $n$ dimensions. I will switch now to the notation of $x_1,x_2,x_3,\cdots,x_n$ as the components, and $e_1,e_2,\cdots,e_n$ as the coordinate bases. The classic divergence formula follows if I can show that: $$\lim_{h\to0}\frac{1}{(2h)^n}\underset{\text{$n-1$ times}}{\int\cdots\int_{-h}^h}f_{x_1}(p_0+he_1+t_2e_2+\cdots+t_ne_n)\hspace{60pt}\\\hspace{100pt}-f_{x_1}(p_0-he_1+t_2e_2+\cdots+t_ne_n)\,\d t_2\,\d t_3\,\cdots\,\d t_{n}$$Evaluates to $\frac{\partial f_{x_1}}{\partial x_1}$.
By Fubini's theorem, assuming "nice" $f$ as we usually do in vector calculus, I can split this into $n-1$ single integrals, and I believe I can apply the mean value theorem $n-1$ times, by considering the variable $t_k$ free and all other $t$ fixed, with $f_{x_1}(\cdots+t_ke_k+\cdots)$ a $C^1$ function from $\Bbb R\to\Bbb R$, obtaining mean values $\vartheta_{t^1_2,t^2_2,t^1_3,t^2_3,\cdots,t^2_n}\in(-h,h)$, since there are two instances of $f_{x_1}$ featuring each $t$.
Using the exact same differentiability argument, noting that I get a factor of $(2h)^{n-1}$ in the iterated application of MVT, I get that the expression inside the limit is within $\epsilon$, for any $\epsilon\gt0$, of: $$\partial_{x_1}f_{x_1}(p_0)+\sum_{k=2}^n\frac{f_{x_1}(p_0+\vartheta_{t^1_k})-f_{x_1}(p_0+\vartheta_{t^2_k})}{2h}\color{red}{=\partial_{x_1}f_{x_1}(p_0)\,\because\vartheta_{t^1_k}=\vartheta_{t^2_k}}$$
I here have the same difficulty. It is not clear to me how to prove that the right hand sum should go to zero, but it must go to zero in order for the divergence formula to be correct.
Any assistance with this mean value issue would be greatly appreciated, and justifying/proving the existence of the limits in the first place - otherwise the use of only cubes as opposed to generic volumes cannot be justified. Alternatively, if the reader is familiar with another (rigorous!) derivation that'd also be appreciated. I do not have access to any vector calculus textbooks, and all these workings are my own ramblings, frustratingly close to the right answer.
