Well, you could take that term from your first post:
F = \vec{\eta}\cdot\vec{\nabla}
where I let
\vec{\eta} = \vec{x}-\vec{x}_0
and multiply it by the Identity Matrix, I,
I = \sum_k{\left|k\right\rangle\left\langle k\right|}
from the left and from the right:
D \equiv IFI = \sum_{i,j}{\left|i\right\rangle\left\langle i\right| \vec{\eta}\cdot\vec{\nabla} \left|j\right\rangle\left\langle j\right|}
= \sum_{i,j}{\left|i\right\rangle\left\langle i\right| \eta_j\frac{\partial}{\partial x_j}}
Notice that, then, the vector function \vec{f}\left(\vec{x}\right) is written as:
\vec{f}\left(\vec{x}\right) = \sum_k{\left|k\right\rangle f_k\left(\vec{x}\right)}
so that
D\vec{f} = \sum_{i,j,k}{\left|i\right\rangle\left\langle i\right| \eta_j\frac{\partial f_k}{\partial x_j} \left|k\right\rangle} = \sum_{i,j}{\eta_j\frac{\partial f_i}{\partial x_j} \left|i\right\rangle}
Notice also that D^2 is, explicitly:
D^2 = \sum_i{\left(\vec{\eta}\cdot\vec{\nabla}\right)^2 \left|i\right\rangle\left\langle i\right|}
= \sum_i{\left(\vec{\eta}\cdot\vec{\nabla}\right) \left(\vec{\eta}\cdot\vec{\nabla}\right) \left|i\right\rangle\left\langle i\right|}
= \sum_{i,j,k}{\eta_j\partial _j\left(\eta_k\partial _k\right) \left|i\right\rangle\left\langle i\right|}
Assuming \partial _i=\frac{\partial}{\partial x_i} (then \partial _j\eta_i=\delta_{ij} -- Kronecker delta); thus:
D^2 = \sum_{i,j,k}{\eta_j\partial _j\left(\eta_k\partial _k\right) \left|i\right\rangle\left\langle i\right|}
= \sum_{i,j,k}{\eta_j\left(\partial _j\eta_k\partial _k + \eta_k\partial _{jk}\right) \left|i\right\rangle\left\langle i\right|}
= \sum_{i,j,k}{\eta_j\partial _j\eta_k\partial _k \left|i\right\rangle\left\langle i\right|} + \sum_{i,j,k}{\eta_j\eta_k\partial _{jk} \left|i\right\rangle\left\langle i\right|}
= \sum_{i,j}{\eta_j\partial _j \left|i\right\rangle\left\langle i\right|} + \sum_{i,j,k}{\eta_j\eta_k\partial _{jk} \left|i\right\rangle\left\langle i\right|}
= D + \sum_{i,j}{\eta_j\eta_k\partial _{jk} \left|i\right\rangle\left\langle i\right|}
And finally:
\sum_{i,j,k}{\eta_j\eta_k\partial _{jk} \left|i\right\rangle\left\langle i\right|} = \sum_{i,j,k}{\eta_j \eta_k \frac{\partial^2}{\partial x_j \partial x_k} \left|i\right\rangle\left\langle i\right|}
= D^2-D=D(D-1)\equiv H
Then you may define H: (which I'll call H because of the similarity with the
Hessian Matrix):
H \equiv D(D-1) = \sum_i{\left|i\right\rangle \vec{\eta}\cdot\vec{\nabla} \left(\vec{\eta}\cdot\vec{\nabla}-1\right)\left\langle i\right|} = \sum_{i,j,k}{\eta_j \eta_k \frac{\partial^2}{\partial x_j \partial x_k} \left|i\right\rangle\left\langle i\right|}
so that:
H\vec{f} = \sum_{i,j,k,l}{\eta_j \eta_k \frac{\partial^2 f_l}{\partial x_j \partial x_k} \left| i \right\rangle \left\langle i \right|\left. l \right\rangle} = \sum_{i,j,k}{\eta_j \eta_k \frac{\partial^2 f_i}{\partial x_j \partial x_k} \left| i \right\rangle}
because \left\langle i \right|\left. l \right\rangle=\delta_{il}.
Now, your expansion may be written as:
\vec{f}\left(\vec{x}_0+\vec{\eta}\right) = \vec{f}\left(\vec{x}_0\right) + D\left.\vec{f}\right|_{\vec{x}=\vec{x}_0} + \frac{1}{2}H\left.\vec{f}\right|_{\vec{x}=\vec{x}_0} + O\left(\vec{\eta}^3\right)