How to get the n-order Taylor expansion of bivariate function?

Question

I want to get the n-order Taylor expansion of a bivariate function at point (x0,y0):

f[x_, y_]: = E^(x + y);
    
{x0,y0}={0,0};

The result calculated by hand is:

$\begin{aligned} \mathrm{e}^{x+y}=& 1+(x+y)+\frac{1}{2 !}\left(x^{2}+2 x y+y^{2}\right)+\frac{1}{3 !}\left(x^{3}+3 x^{2} y+3 x y^{2}+y^{3}\right) \\ &+\cdots+\frac{1}{n !}(x+y)^{n}+R_{n}=\sum_{k=0}^{n} \frac{(x+y)^{k}}{k !}+R_{n} \end{aligned}$

$R_{n}$ is Langrange Remainder. How to get this result by MMA?

ps:

1. There are some answers for reference, but they are all about univariate functions:

Taylor series representation as an infinite sum

Series with specific notation

For example,

E^x=

    Clear["Global`*"]
    series[expr_, x_, x0_] := 
     Defer[expr = Sum[#, {n, 0, \[Infinity]}]] &[
      FullSimplify@
        SeriesCoefficient[expr, {x, x0, n}, 
         Assumptions -> {n >= 0}] (x - x0)^n]
    series[E^x, x, 0]

($e^{x}=\sum_{n=0}^{\infty} \frac{x^{n}}{n !}$)

Similarly, how to get the n-order Taylor expansion of "E^(x + y)"?

2. Taylor series representation as an finite sum for bivariate function：

Taylor Expansion of and Exponential function but multivariable

https://community.wolfram.com/groups/-/m/t/2353053?p_p_auth=8FpuUkxs

For example,

3-order Taylor expansion for E^(x + y):

Clear["Global`*"]
f[x_, y_] = E^(x + y);

ef[x_, y_, x0_, y0_, n_Integer] := 
 Normal[Series[f[(x - x0)*t + x0, (y - y0)*t + y0], {t, 0, n}]] /. 
  t -> 1

ef[x, y, 0, 0, 3]

(*1 + x + y + 1/2 (x + y)^2 + 1/6 (x + y)^3*)

Similarly, how to get the "n"-order Taylor expansion of E^(x + y)?

In the comment area, the results can be obtained by the substitution method. Thanks for the comment from @Michael E2. But I want to ask, is there any method suitable for the general situation? Because most bivariate functions cannot be converted into univariate functions by the method of substitution.

For example, f[x_, y_] = E^x*Log[1 + y].

ps. The n-order Taylor expansion formula of bivariate function f(x,y) at (x0,y0).

(x0 + h, y0 + k) is any point in the neighborhood of point (x0, y0).

$\begin{aligned} & f\left(x_{0}+h, y_{0}+k\right) \\=& f\left(x_{0}, y_{0}\right)+\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right) f\left(x_{0}, y_{0}\right)+\\ & \frac{1}{2 !}\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right)^{2} f\left(x_{0}, y_{0}\right)+\cdots+\frac{1}{n !}\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right)^{n} f\left(x_{0}, y_{0}\right) \end{aligned}$

$\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right) f\left(x_{0}, y_{0}\right)$ ---> $h f_{x}\left(x_{0}, y_{0}\right)+k f_{y}\left(x_{0}, y_{0}\right)$

$\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right)^{2} f\left(x_{0}, y_{0}\right)$ ---> $h^{2} f_{x x}\left(x_{0}, y_{0}\right)+2 h k f_{x y}\left(x_{0}, y_{0}\right)+k^{2} f_{y y}\left(x_{0}, y_{0}\right)$

$\left(h \frac{\partial}{\partial x}+k \frac{\partial}{\partial y}\right)^{m} f\left(x_{0}, y_{0}\right)$ ---> $\left.\sum_{p=0}^{m} C_{m}^{p} h^{p} k^{m-p} \frac{\partial^{m} f}{\partial x^{p} \partial y^{m-p}}\right|_{\left(x_{0}, y_{0}\right)}$