Big-O Notation for remainder terms in Taylor expansion The Big-O notation is commonly used in Taylor expansions of the form
$$f(x+\epsilon)=f(x)+\epsilon f'(x)+O(\epsilon^2)$$ to say that the remainder term grows at least quadratic around $\epsilon=0$.
Considering the analogue case of a function $f:\mathbb{R}^2\to\mathbb{R}$:
$$f(x+\epsilon_1,y+\epsilon_2)=f(x,y)+\epsilon_1 f_x(x,y)+\epsilon_2 f_y(x,y)+\epsilon_1\epsilon_2 f_{x,y}(x,y)+O(?)$$
my question is: What would be a suitable notation for the remainder term? 
The notation should say precisely that the term contains only contributions of order $\geq2$ in at least one of the variables $\epsilon_1,\epsilon_2$ around $(\epsilon_1,\epsilon_2)=(0,0)$.
I am also not sure if Big-O notation is the right thing to use here, as i read here its generalisation to 2D is quite nontrivial. However, the difficulties of this issue might be unimportant in my context, because they arise when describing the asymptotic behavior $\epsilon\to\infty$, whereas my question is on describing the local behaviour of the function for $\epsilon\to0$.
 A: Consider $f=f(x_1,x_2,\ldots,x_N)$ be a function of $N$ variables in a open convex set $D$ of $\textbf{R}^{N}$. Then the directional derivative of $f$ allong $\overline{h}$ is
$$
\nabla_{\overline{h}}f=\left\langle \nabla  f,\overline{h}\right\rangle=\sum^{N}_{k=1}\frac{\partial f}{\partial x_k}h_k,
$$
where
$$
\nabla=\overline{e}_1\frac{\partial}{\partial x_1}+\overline{e}_2\frac{\partial}{\partial x_2}+\ldots+\overline{e}_{N}\frac{\partial}{\partial x_{N}}
$$
Also symbolical we write $\left\langle \nabla f,\overline{h}\right\rangle=\left\langle\nabla,\overline{h}\right\rangle f$ and it holds
$$
\nabla^{\nu}_{\overline{h}}f=\left\langle \nabla,\overline{h}\right\rangle^{\nu} f
$$
For eample with $N=3$, then
$$
\nabla^2_{\overline{h}}f=\ldots=\sum^{3}_{i=1}h_i^2\frac{\partial^2 f}{\partial x_i^2}+2\sum^{3}_{i<j}h_ih_j\frac{\partial^2 f}{\partial x_i\partial x_j}=\left(\sum^{3}_{i=1}h_i\frac{\partial}{\partial x_i}\right)^2f
$$
and in general for $N=3$ we have
$$
\nabla^{\nu}_{\overline{h}}f=\left(\sum^{3}_{i=1}h_i\frac{\partial}{\partial x_i}\right)^{\nu}f=\sum_{\begin{array}{cc}
\nu_1+\nu_2+\nu_3=\nu\\
\nu_1,\nu_2,\nu_3\geq 0
\end{array}}\frac{\nu !}{\nu_1!\nu_2!\nu_3!}h_1^{\nu_1}h_2^{\nu_2}h_3^{\nu_3}\frac{\partial^{\nu} f}{\partial x^{\nu_1}_1\partial x^{\nu_2}_2\partial x^{\nu_3}_3}.
$$
Theorem.
If $f$ is differentiatable function of at least $n+1$ order derivatives in $D\subset \textbf{R}^N$. Then the Taylor expansion is
$$
f(P)=f(P_0)+\sum^{n}_{\nu=1}\frac{1}{\nu!}\left\langle \nabla,(P-P_0)\right\rangle^{\nu}f(P_0)+R_n(h),
$$
where
$$
R_n(h)=\frac{1}{(n+1)!}\left\langle \nabla,(P-P_0)\right\rangle^{n+1}f(P^*),
$$
where $P^*$ is a specific point of the $PP_0$ segment. Also if $M_{n+1}$ is the upper bound of the absolute value of the $n+1$ order derivatives in the "ball" $B(P_0,|h|)$, $|h|=|P-P_0|$, then
$$
|R_{n}(h)|\leq \frac{(h\sqrt{N})^{n+1}}{(n+1)!}M_{n+1}
$$
Example.
For $N=2$ we have $P_0=(a,b)$
$$
f(x,y)=f(P_0)+(x-a)\left(\frac{\partial f}{\partial x}\right)_{P_0}+(y-b)\left(\frac{\partial f}{\partial y}\right)_{P_0}+
$$
$$
+\frac{1}{2}\left[(x-a)^2\left(\frac{\partial^2 f}{\partial x^2}\right)_{P_0}+2(x-a)(y-b)\left(\frac{\partial^2 f}{\partial x\partial y}\right)_{P_0}+(y-b)^2\left(\frac{\partial^2 f}{\partial y^2}\right)_{P_0}\right]+R_2,
$$
where
$$
R_2=\frac{(|h|\sqrt{2})^3}{3!}M_{3},
$$
where
$$
M_3=\textrm{max}\left|\frac{\partial^3 }{\partial^i x\partial^j y}f(x,y)\right|\textrm{, }i+j=3\textrm{, }i,j\geq 0.
$$
Hence if $P_0=(x_0,y_0)$ and $x=x_0+\epsilon_1$, $y=y_0+\epsilon_2$, then
$$
f(x,y)=f(P_0)+f'_{x}(P_0)\epsilon_1+f'_{y}(P_0)\epsilon_2+
$$
$$
+\frac{1}{2}\left[f''_{xx}(P_0)\epsilon_1^2+2f''_{xy}(P_0)\epsilon_1\epsilon_2+f''_{yy}(P_0)\epsilon_2^2\right]+O\left(|h|^{3}\right),
$$
where $|h|=\sqrt{(x-x_0)^2+(y-y_0)^2}=\sqrt{\epsilon_1^2+\epsilon_2^2}$. If the 3rd order derivatives are bounded in $D$.
