Derivation of multivariable Taylor series

Let $\phi(\boldsymbol{r})$ be a scalar field, and $\boldsymbol{a} \cdot \nabla \phi$ gives the directional derivative of $\phi$ in the direction of $a$. That is,

$$\boldsymbol{a} \cdot \nabla \phi(\boldsymbol{r}) = \lim_{t\to 0} \frac{\phi(\boldsymbol{r} + \boldsymbol{a} t) - \phi(\boldsymbol{r})}{t}$$

Now let's consider $\Phi(t) = \phi(\boldsymbol{r}_0 + \boldsymbol{a}t)$ for some finite $t$. Now, let's expand this in powers of $t$. This is a one-dimensional Taylor series.

$$\Phi(t) = \Phi(0) + \Phi'(0)t + \frac{1}{2!} \Phi''(0) t^2 + \ldots$$

To substitute back in $\Phi(t) = \phi(\boldsymbol{r}_0+\boldsymbol{a}t)$, we must compute derivatives of $\Phi$ in terms of $\phi$. Again, we resort to the basic definition of the derivative.

$$\Phi'(0) = \lim_{t\to 0} \frac{\phi(\boldsymbol{r}_0+\boldsymbol{a}t) - \phi(\boldsymbol{r}_0)}{t} = \boldsymbol{a} \cdot \nabla \phi(\boldsymbol{r})\Big|_{\boldsymbol{r}=\boldsymbol{r}_0}$$

And similarly for higher derivatives. This enables us to write,

$$\phi(\boldsymbol{r}_0+\boldsymbol{a}t) = \phi(\boldsymbol{r}_0) + [\boldsymbol{a} \cdot \nabla \phi(\boldsymbol{r})] \Big|_{\boldsymbol{r}=\boldsymbol{r}_0} t + \frac{1}{2!} [\boldsymbol{a} \cdot \nabla][\boldsymbol{a} \cdot \nabla]\phi(\boldsymbol{r}) \Big|_{\boldsymbol{r}=\boldsymbol{r}_0} t^2 + \ldots$$

It is not difficult to show that this form reproduces the form of the original question. Take $t=1$ and let $\boldsymbol{a} = (x-x_0, y-y_0)$ and $\boldsymbol{r}_0 = (x_0, y_0)$. Thus, we have built multivariate Taylor series from the well-established case of a single variable, just by use of the directional derivative.

I think the easiest way to understand this is coming from the place of operators and linear transformations. A Taylor series in one dimension can be understood by exponentiating the derivative operator:

$$ f(x+a) = e^{a\frac{d}{dx}}f(x) = f(x) + af^\prime(x) + \frac{1}{2!}a^2f^{\prime\prime}(x)+... $$

You can see this in one way as follows. The infinitesimal (linear order) transformation $f(x+dx) = f(x) + dx f^\prime(x)$ is known, and we can build up the finite transformation by an infinite succession of infinitesimal transformations:

$$ f(x+a) = \lim_{N\rightarrow\infty} \left(1+\frac{a}{N}\frac{d}{dx}\right)^N f(x) = e^{a\frac{d}{dx}} f(x). $$

It is straightforward to extend this to multiple variables if we know the infinitesimal transformation (sometimes referred to as the generator), which you intuitively know as, $f(x+dx, y+dy) = f(x,y) + dx\frac{\partial}{\partial x}f(x,y) + dy\frac{\partial}{\partial y}f(x,y)$.

The finite transformation is then, $$ f(x+a,y+b) = e^{a\frac{\partial}{\partial x}+b\frac{\partial}{\partial y}} f(x,y)\\ = \left[1+a\frac{\partial}{\partial x} + b\frac{\partial}{\partial y} + \frac{1}{2!}\left(a^2\frac{\partial^2}{\partial x^2} + 2a\frac{\partial}{\partial x}b\frac{\partial}{\partial y}+ b^2\frac{\partial^2}{\partial y^2}\right) + ...\right]f(x,y). $$

Let $u \in \mathbb{R}^m, \, h \in \mathbb{R}^m, \, t \in \mathbb{R},$ and $F(t)=f(u+th).$ Suppose that $F$ can be expanded into Taylor's series $$F(t)=\sum\limits_{n=0}^{\infty}{\frac{1}{n!}}F^{(n)}(0)t^n.\tag{*}$$ Taylor's expansion for $f$ can be obtained from $({}^{*})$ by differentiating $f$ and then put $t=1$.

For the case $n=2$ $$f(u+h)=\sum\limits_{n=0}^{\infty}{{\frac{1}{n!}}d^{n}f(u)},$$ where $u=(x, \, y)\quad h=(dx,\, dy),$ $$d^{n}f(u)=\sum\limits_{k=0}^{n}{\binom{n}{k}}\frac{\partial^n{f}}{\partial{x}^k {}\partial{y}^{n-k}}dx^kdy^{n-k}.$$