The quadratic multivariable Taylor approximation of a function $f(x, y)$ around a point $(a, b)$ is given by $f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b) + \frac{1}{2}f_{xx}(a, b)(x - a)^2 + f_{xy}(a, b)(x - a)(y - b) + \frac{1}{2}f_{yy}(a, b)(y - b)^2$ There is a vectorized version of this expression which works for functions and points with an arbitrary number of variables, $f(\vec{a})\ + ∇f(\vec{a})\cdot(\vec{x} - \vec{a}) + \frac{1}{2}(\vec{x} - \vec{a})^TH(\vec{a})(\vec{x} - \vec{a})$, where the elements of $\vec{a}$ are the coordinates of the point around which we're approximating, the elements of $\vec{x}$ are the inputs of $f$, and $H$ is the Hessian matrix of $f$.
Following the pattern of single variable Taylor approximations, the cubic multivariable Taylor approximation should be $f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b) + \frac{1}{2}f_{xx}(a, b)(x - a)^2 + f_{xy}(a, b)(x - a)(y - b) + \frac{1}{2}f_{yy}(a, b)(y - b)^2 + \frac{1}{3!}f_{xxx}(a, b)(x - a)^3 + \frac{3}{3!}f_{xxy}(a, b)(x - a)^2(y - b) + \frac{3}{3!}f_{xyy}(a, b)(x - a)(y - b)^2 + \frac{1}{3!}f_{yyy}(a, b)(y - b)^3$
The amount of simplification in my coefficients is inconsistent, and I'm intentionally not generalizing the approximation order by introducing additional variables or sigma notation, as this expression already requires more than enough mental unpacking in my opinion, and it can only get worse in vector notation. Instead I'll use an image and words to describe the pattern I have in mind.
Every unique partial derivative up to the order of the approximation constitutes a term (third order approximation includes all zeroth, first, second, and third partial derivatives), and the numerator of the fraction coefficient for each term is the number of equivalent permutations of that partial derivative ($f_{xyy} = f_{yxy} = f_{yyx}$, so the numberator is $3$).
Is this the correct extension of the quadratic multivariable Taylor approximation to cubic approximations, and a correct pattern for writing quartic and higher approximations, and how can we put this into vector notation? An optimal answer, in order to avoid layering on too much abstraction, will likely include the expression for a cubic approximation in vector notation (possibly also quartic, if helpful) and an explanation in words of how to continue the pattern into higher order approximations (i.e., do we require an increasingly higher order "tensor," a "Hyperhessian" of sorts, as we continue to add more terms?).
Your expression is correct.
If you want to write it in "vector notation" we simply use the usual way of writing derivatives. Denote by $f^{(p)}$ the $p$th derivative of $f$ (which, by the way, is a $p$-linear continuous function) and write $h^{(p)}$ to mean the vector $(h, \ldots, h)$ ($h$ appearing $p$ times). Then, the Taylor polynomial of $f$ centred at $a$ of degree $n$ is $$\mathbf{T}_nf(a) \cdot h = f(a) + f'(a) \cdot h + \dfrac{f''(a) \cdot h^{(2)}}{2!} + \ldots + \dfrac{f^{(n)}(a) \cdot h^{(n)}}{n!}.$$