For the purposes of this question, assume that ordered pairs, vectors, and column matrices are all equivalent.
$$\newcommand{r}{\mathbf{r}}\newcommand{b}{\mathbf{b}}\newcommand{H}{\boldsymbol{H}}$$
Consider the following table, which summarizes the $n$th terms of the Taylor expansions of the univariate function $g(x)=\sum u_n(x)$ about the point $x=a$ and the bivariate function $F(x,y)=F(\r)=\sum Z_n(\r)$ about the point $(x,y)=(h,k)=\b$:
$$\begin{array}{c|cc} n & u_n(x) & Z_n(\r) \\ \hline 0 & g(a) & F(\b) \\ 1 & g'(a)\,(x-a) & \nabla F(\b)\,(\r-\b) \\ 2 & g''(a)/2\,(x-a)^2 & \H F(\b)\,(\r-\b)^2 \\ 3 & g'''(a)/3!\,(x-a)^3 & ? \\ \vdots & \vdots & \vdots \\ t & g^{(t)}(a)/t!\,(x-a)^t & ? \\ \end{array}$$
where $\H F(\b)$ is the Hessian matrix of $F$ evaluated at $\b$,
$$\H F(\b) = \begin{pmatrix} F_{xx}(h,k) & F_{xy}(h,k) \\ F_{yx}(h,k) & F_{yy}(h,k) \\ \end{pmatrix}$$
How is the definition of $Z_n(\r)$ extended for values $n\ge3$? My gut tells me that it involves $n$th partial derivatives of $F$, that it is not pretty, and that there is no compact notation for it. Also, I am pretty sure that the order of my matrix terms is out of whack for multiplication to be defined, so that might also be messing with my understanding.
The $n$-th order terms of the Taylor polynomial of a sufficiently differentiable bivariate function $f(x,y)$ are determined by the $2^n$ mixed partial derivatives of order $n$.
You can think of these $2^n$ derivatives as encoded in a
$$\underbrace{2\times\cdots\times2}_{n\,\text{times}}$$
array. When $n=1$ this is a vector with $2$ components. When $n=2$ this is a $2\times 2$ matrix, the Hessian, with $4$ components. When $n=3$ we have a $2\times2\times2$ "cube" of third partials, with $8$ components. When $n=4$ we have a $2\times2\times2\times2$ four-dimensional hypercube of fourth partials, with $16$ components. And so on. (This is a general result: for a function of $m$ independent variables instead of $2$, just replace $2$ with $m$.) If you insist on taking this perspective (of $n$-dimensional arrays), you must talk about tensors of order $n$. It's a little more down-to-earth, if less succinct, to think about the direct expansion in terms of the mixed partial derivatives, instead of in terms of tensors.
By Clairaut's theorem on the equality of mixed partials, many of these $2^n$ mixed partials are equivalent. For example, $f_{xy}=f_{yx}$, and $f_{xxy}=f_{xyx}=f_{yxx}$, and $f_{xyy}=f_{yxy}=f_{yyx}$. (I'm glossing over some technicalities about the domain on which these partial derivatives are defined.) So the arrays of mixed partials have a great deal of symmetry to them. In general, there will only be $n+1$ different mixed partials of order $n$. When $n=2$, there's $3$ distinct mixed partials: $f_{xx}$, $f_{xy}$, and $f_{yy}$. When $n=3$, there's $4$ distinct mixed partials: $f_{xxx}$, $f_{xxy}$, $f_{xyy}$, and $f_{yyy}$. This symmetry greatly simplifies the Taylor expansion. The symmetry is similar, but harder to notate, when you have $m$ independent variables instead of $2$.
Now let's work out the bivariate case explicitly.
A specific mixed partial of order $n$, say
$$f_{\underbrace{x\cdots x}_{i}\underbrace{y\cdots y}_{j}}=\frac{\partial f}{\partial x^i\,y^j}:=f_{i,j}$$
where $i+j=n$, contributes a term
$$\frac{1}{n!}f_{i,j}(x-x_0)^i(y-y_0)^j$$
to the Taylor expansion. But ${n\choose i}$ of these terms are the same, by the symmetry mentioned above. (We can choose any of the $i$ places for $x$.)
So the $n$-th order terms of the Taylor expansion (there are $n+1$ one of them) look like
$$\frac{1}{n!}f_{n,0}(x-x_0)^n+\frac{1}{n!}{n\choose 1}f_{n-1,1}(x-x_0)^{n-1}(y-y_0)+\cdots+\frac{1}{n!}{n\choose n-1}f_{1,n-1}(x-x_0)(y-y_0)^{n-1}+\frac{1}{n!}f_{0,n}(y-y_0)^n$$
or, more succinctly,
$$\frac{1}{n!}\sum_{i=0}^n{n\choose i}f_{n-i, i}(x-x_0)^{n-i}(y-y_0)^i$$
Note that when I write $f_{n-i,i}$, I mean the mixed partial evaluated at the point $(x_0, y_0)$.