Gradient of the second order term of Newton's Method 
I know that Netwon's method can be pushed to the second order using the 1st Taylor expansion. However, how can I generalize Netwon's method to take x_0 as a vector and have the ability take the gradient?

The main point here is that $\nabla^2f(\mathbf x_0)$ is symmetric (because $\frac{\delta^2}{\delta x_i\delta x_k} = \frac{\delta^2}{\delta x_k\delta x_i}$).
Let's presume the space is $d$-dimensional and consider a Cartesian coordinate system, such that the gradient can be defined as $\nabla g := (\frac{\delta g}{\delta x_1}, \frac{\delta g}{\delta x_2}, \ldots, \frac{\delta g}{\delta x_d})^t$.
Below, I will abbreviate $\nabla^2f(\mathbf x_0)$ with the constant matrix $\mathbf F = (F_{ij})$ and I will use the Kronecker symbol $\delta_{ij}$:
$$
\begin{align}
\frac{\delta}{\delta x_k} \big[(\mathbf x - \mathbf x_0) \mathbf F (\mathbf x - \mathbf x_0)\big] & = \frac{\delta}{\delta x_k} \big[\sum_i \sum_j (x_i - x_{0i}) F_{ij} (x_j - x_{0j})\big]\\
    & = \sum_i \sum_j \delta_{ik} F_{ij} (x_j - x_{0j}) + \sum_i \sum_j (x_i - x_{0i}) F_{ij}\delta_{jk}) \\
    & = \sum_j F_{kj} (x_j - x_{0j}) + \sum_i (x_i - x_{0i}) F_{ik} \\
    & = \sum_j F_{kj} (x_j - x_{0j}) + \sum_i F_{ki} (x_i - x_{0i}) \\
    & = 2 \sum_j F_{kj} (x_j - x_{0j}),
\end{align}
$$
where I used the product rule of derivation in the second line and the symmetry of $\mathbf F$ in the fourth line.
This gives your required identity if you let $k$ run through all indices $k=1,\ldots,d$.
