Xrfgtjtk

While researching multivariable Taylor series I found this expression for a Taylor polynomial in two variables:
$$T_n(x,y)=sum_{i=1}^nleft((x-x_0)frac{partial}{partial x}+(y-y_0)frac{partial}{partial y}right)^ifrac{f(x_0,y_0)}{i!}$$
As an example, the second-degree Taylor polynomial of $f(x,y)=frac{1}{xy}$ at $(1,1)$ is
$$f(1,1)+(x-1)frac{partial f}{partial x}(1,1)+(y-1)frac{partial f}{partial y}(1,1)+frac{1}{2}left((x-1)^2frac{partial^2f}{partial x^2}(1,1)+2(x-1)(y-1)frac{partial^2f}{partial x partial y}(1,1)+(y-1)^2frac{partial^2f}{partial y^2}(1,1)right)$$
$$=1-(x-1)-(y-1)+(x-1)^2+(x-1)(y-1)+(y-1)^2$$
However this notation seems sort of "cheaty" to me. Can we use the differential operator in this way? Namely, is raising the differential operator to a power and then "multiplying" it against the function at the desired point in this manner legitimate?

This is completely legitimate.
The multivariable Taylor formula can be written as
$$
T_m(a+Delta x)= sum_{k=0}^{m} frac{d^{k}f(a)}{k!}+dfrac{d^{m+1}f(a+theta Delta x)}{(m+1)!},quad
thetain (0,1).
$$
We have
$$
df(x)= frac{partial f(x)}{partial x_1} Delta x_1+ldots+frac{partial f(x)}{partial x_n}Delta x_n,
$$
where $Delta x=(Delta x_1,ldots,Delta x_n)$ is a known constant vector and $x$ is a point where the differential is calculated.

The second order differential is defined as
$$
d^2f(x)= d(df(x))= frac{partial df(x)}{partial x_1}Delta x_1+ldots+frac{partial df(x)}{partial x_n}Delta x_n
$$
$$
=frac{partial left( frac{partial f}{partial x_1}Delta x_1+ldots+frac{partial f}{partial x_n}Delta x_nright) }{partial x_1}Delta x_1+ldots+frac{partial left( frac{partial f}{partial x_1}Delta x_1+ldots+frac{partial f}{partial x_n}Delta x_nright) }{partial x_n}Delta x_n
$$
$$
=sum_{i=1}^{n}sum_{j=1}^{n} frac{partial^2 f}{partial x_i partial x_j}Delta x_iDelta x_j.
$$
This shows that if we consider the differential as a differential operator
$$
df(x)= left( frac{partial}{partial x_1}Delta x_1+ldots+frac{partial}{partial x_n}Delta x_n right)f(x),
$$
then we can write
$$
d^2 f(x)=left( frac{partial}{partial x_1}Delta x_1+ldots+frac{partial}{partial x_n}Delta x_n right)^2f(x).
$$
Similarly, one can obtain formulas for differentials of higher orders:
$$
d^3f(x)= d(d^{2}f(x))= left( frac{partial}{partial x_1}Delta x_1+ldots+frac{partial}{partial x_n}Delta x_n right)^{3}f(x)
$$
$$
ldots
$$
$$
d^mf(x)= d(d^{m-1}f(x))= left( frac{partial}{partial x_1}Delta x_1+ldots+frac{partial}{partial x_n}Delta x_n right)^{m}f(x)
$$
In the particular case of two variables, we get the
formula from Wikipedia
$$
d^mf = sum_{k=0}^m binom{m}{k}frac{partial^m f}{partial x^k partial y^{m-k}}(dx)^k(dy)^{m-k}.
$$

Partial derivatives, and differential operators in general, are linear operators; raising them to powers or multiplying makes perfect sense, as composition.

There's only one problem here - that multiplication by $(x-x_0)$ and $(y-y_0)$. Those don't commute with the partial derivatives. For example, taking it all literally as composition of linear operators, I get
begin{align*}left((x-x_0)frac{partial}{partial x} + (y-y_0)frac{partial}{partial y}right)^2 &= (x-x_0)^2frac{partial^2}{partial x^2} +(x-x_0)(y-y_0)left(frac{partial^2}{partial xpartial y}+frac{partial^2}{partial ypartial x}right)cdots\ &quad +(x-x_0)frac{partial}{partial x}+(y-y_0)frac{partial}{partial y}end{align*}
That bit with the mixed partials isn't too bad; there's a standard theorem that for reasonably nice functions they're the same. The problem is the terms on the second line, which come from differentiating the $(x-x_0)$ and $(y-y_0)$ when we apply the product rule. Those terms aren't in the Taylor expansion.

To see where this goes, apply this formula to finding the Taylor series of a linear function $f(x)=ax+by$ at the origin. We have $left(xfrac{partial}{partial x}+yfrac{partial}{partial y}right)f(0,0) = ax+by=f(x,y)$ (that $f(x_0,y_0)$ in the formula actually means to evaluate the derivatives there), so $left(xfrac{partial}{partial x}+yfrac{partial}{partial y}right)^if(0,0) = ax+by$ for all $ige 1$. Add up the series, and we get a series $0+(ax+by)sum_{i=1}^{infty}frac{1}{i!} = (e-1)(ax+by)$. Absurd.

As noted in the comments below - there is a fix. The $frac{partial}{partial x}$ operator is actually "differentiate with respect to the first variable", and the $x$ in $(x-x_0)$ is a new variable. Time to rename things. Call the point we're evaluating the Taylor series at $(u,v)$, so the operator we're raising to various powers is $(u-x_0)frac{partial}{partial x}+(v-y_0)frac{partial}{partial y}$, or the equivalent with more subscripts in higher dimensions. Now the multipliers are constant with respect to the derivatives, and everything works out right with no extra terms. Another one for the hazards of using variables out of scope.