Meaning of cross terms in multivariable Taylor expansion The cross terms in the Taylor expansion of $f = f(x,y)$ in $(x_0,y_0)$ 
$$
f(x,y) = f(x_0,y_0) + \ldots + \frac{1}{2!}\bigg( \frac{\partial^2 f}{\partial x^2} (\Delta x)^2 + \color{green}{2\frac{\partial^2 f}{\partial x \partial y} \Delta x \Delta y} + \frac{\partial^2 f}{\partial y^2}(\Delta y)^2 \bigg) + \ldots \tag{1}
$$ can be seen as to have arosen from the cross terms of 
$$
f(x,y) = \sum_{n=0}^\infty \frac{1}{n!} \bigg[ \bigg( \Delta x \frac{\partial}{\partial x} + \Delta y \frac{\partial}{\partial y}\bigg)^n \ f(x,y)\bigg]_{x_0,y_0} , \tag{2}
$$ 
where we consider the partial derivatives to only operate on $f$.
I'm missing the interpretation of these cross terms. 
$(2)$ somewhat explains where the (algebraic structure of) these cross terms comes from, but it doesn't give me any insight as to why we need to consider the product of two changes. Moreover, I also don't grasp how my text arrives at $(2)$. Why are the cross terms included?
I understand the operation of mixed derivatives in terms of calculations, but I want to know how I can interpret them.
Furthermore, I also understand the concept of stationary points of multi-variable functions. I don't think the concept of stationary points are relevant here.
I'm not necessarily looking for a geometric interpretation; as long as the relevance of the mixed partial derivative terms is made evident.
 A: One good way to visualize these sorts of formulas is to probe them with smooth curves so that you can invoke your 1-dimensional intuition.
Starting at $(x_0, y_0)$, suppose we move along the velocity vector $(v, w)$.  How do the values of $f$ change?  To answer this question we can look at the values of $f$ along the parametric curve $r(t) = (x_0 + vt, y_0 + wt)$.  By the chain rule we get:
$$\frac{d}{dt} f(r(t)) = f_x(r(t)) v + f_y(r(t))w$$
and:    
$$\frac{d^2}{dt^2} f(r(t)) = \left(f_{xx}(r(t)) v^2 + f_{xy}(r(t)) vw \right) + \left(f_{yx}(r(t)) wv + f_{yy}(r(t)) w^2 \right)$$
$$= f_{xx}(r(t))v^2 + 2f_{xy}(r(t))vw + f_{yy}(r(t))w^2$$
Hopefully this explains the cross terms: they represent the fact that rate of change of $f$ in the $x$ direction - represented by $f_x$ - is itself changing in both the $x$ and $y$ directions, and similarly for $f_y$.  The quadratic behavior in $v$ and $w$ just says that faster motion away from $(x_0, y_0)$ (corresponding to larger displacements $\Delta x$ and $\Delta y$) causes more dramatic changes in the values of $f$.

A separate remark: if you really don't like the cross terms then you can always get rid of them by changing coordinates (assuming $f$ is sufficiently regular).  To see this, write:
$$f_{xx} \Delta x^2 + 2f_{xy} \Delta x \Delta y + f_{yy} \Delta y^2 = (\Delta x, \Delta y) \left(\begin{array}{cc} f_{xx} & f_{xy} \\ f_{xy} & f_{yy} \end{array}\right) \left(\begin{array}{c} \Delta x \\ \Delta y \end{array}\right)\ \ \ $$
The $2 \times 2$ matrix in the middle - called the Hessian matrix of $f$ - is symmetric and hence it can be diagonalized by an orthonormal coordinate system $(u, v)$.  In this coordinate system the second derivative takes the form $f_{uu} \Delta u^2 + f_{vv} \Delta v^2$ - informally, this means that the rate of change in the $u$ direction does not itself change in the $v$ direction, and vice-versa.
