The key is to do the Taylor expansion properly (I change the notation slightly, using initially $\delta$ instead $d$. I also highlight in $\color{red}{red}$ the terms in the Taylor expansion that differ to yours):
First, let's define $\delta V$:
$$\delta V := V(S_1(t_0)+\delta S_1, S_2(t_0)+\delta S_2, t_0+\delta t))- V(S_1(t_0), S_2(t_0), t_0))$$
Now let's do a Taylor expansion around zero (i.e. $S_1(t_0), S_2(t_0), t_0)$ ):
$$\delta V=\frac{\partial V}{\partial t} \delta t + \frac{\partial V}{\partial S_1} \delta S_1 + \frac{\partial V}{\partial S_2} \delta S_2 + \frac{1}{2}\frac{\partial^2 V}{\partial S_1^2} \color{red}{\delta S_1^2} + \frac{\partial^2 V}{\partial S_1 \partial S_2} \color{red}{\delta S_1 \delta S_2} + \frac{1}{2} \frac{\partial^2 V}{\partial S_2^2} \color{red}{\delta S_2^2}$$
We ignore all higher order terms (i.e. $\delta t \delta S_1$, $\delta t \delta S_2$ , $\delta t ^2$ etc.), assuming they go to zero as $\delta \to 0$).
The key results I will be using are:
For any finite $\delta t > 0$ and a standard Brownian motion $X_t$:
\begin{align*} \tag{1} \delta X_t := X(\delta t) \stackrel{d}{=}\sqrt{\delta t}X \end{align*}
\begin{align*} \tag{2} \mathbb{E}[\delta X^2] = \mathbb{E}[\delta t X^2]=\delta t \end{align*}
And as $\delta t \to 0$ we can argue that:
\begin{align*} \tag{3} V(\delta X^2) = V(\delta t X^2]=\delta t^2 V(X^2) \to 0 \end{align*}
(the last results is true because $\delta t^2 \to 0$ and $V(X^2)$ = 3 using Moment Generating Function of a Normal distribution)
Now evaluating the terms:
$$\delta S_1^2=(a_1\delta t + b_1 \delta X_1)^2=a_1^2\delta t^2+2a_1b_1\delta t \delta X_1+b_1^2\delta X_1^2$$
NowUsing (1), (2) and (3) above, we can now argue that as $\delta t \to 0$ that:
- $a_1^2\delta t^2 \to 0$
- $2a_1b_1\delta t \delta X_1 = 2a_1b_1\delta t \sqrt{\delta t} X_1 = 2a_1b_1\delta t^{\frac{3}{2}}X_1\to 0$ (using the scaling property of Brownian motion, i.e. for any finite $\delta t > 0$, we have that $\delta X_1 := X_1(\delta t) \stackrel{d}{=}\sqrt{\delta t}X_1(1)$)
- $b_1^2\delta X_1^2 \to b_1^2 dt$ (using the scaling property again and computing: $\mathbb{E}[b_1^2\delta X_1^2] = b_1^2\delta t\mathbb{E}[X_1^2] = b_1^2\delta t $, whilst the variance converges to zero, i.e. $V(\delta b_1^2\delta X_1^2) = b_1^2\delta t^2 \mathbb{E}[X_1^2] \to 0$ because \delta t^2 \to 0using (3) above)
So using the above $b_1^2\delta S_1^2 \to b_1^2dt$$\delta S_1^2 \to b_1^2dt$
Using the same machinery, we get that $b_2^2 \delta t S_2^2 \to b_2^2 dt$$\delta S_2^2 \to b_2^2 dt$. And we can use the same machinery to compute $\delta S_1 \delta S_2$:
$$\delta S_1 \delta S_2 = a_1a_2 \delta t^2 + a_1b_2 \delta t \delta X_1+b_1a_2 \delta t \delta X_1+ b_1 b_2 \delta X_1 \delta X_2$$
Using the results already shown, everything goes to zero except for the last term:
$$b_1b_2\delta X_1 \delta X_2= b_1b_2 X_1(\delta t) X_2 (\delta t) = b_1b_2 \sqrt{\delta}X_1(1) \sqrt{\delta} X_2 (1) = b1b2 \delta t X_1X_2 $$
We know that the variances are (by definition of Brownian motion) $V(X_1(t))=\mathbb{E}[X_1(t)^2]= t$ and $V(X_2(t))=\mathbb{E}[X_1(t)^2]= t$ and the correlation between $X_1(t)$ and $X_2(t)$ is $\rho$, which means that:
$$Cov(X_1(t), X_2(t))=\rho \sqrt{V(X_1)V(X_2)}= \rho t$$
We also know that:
$$Cov(X_1(t), X_2(t)) = \mathbb{E}[X_1(t) X_2(t)]-\mathbb{E}[X_1(t)]\mathbb{E}[X_1(t)] = \mathbb{E}[X_1(t) X_2(t)]=t\mathbb{E}[X_1 X_2] $$
Using the last two equations, we can say that:
$$\mathbb{E}[X_1 X_2]=\rho $$
Therefore, the expected value of the cross term $\delta S_1 \delta S_2$ is $b_1b_2\delta t \rho$ which goes to $b_1b_2 \rho dt$ as $\delta t \to 0$. The variance of this term goes to $0$ (using $V( b1b2 \delta t X_1X_2)=b_1^2b_2^2\delta t^2 V(X_1X_1)$ which goes to zero due to $\delta t^2 \to 0$), so this term converges to $b_1b_2 \rho dt$.
So the formula is (after grouping all terms):
$$dV =\left( \frac{\partial V}{\partial t} + \frac{\partial V}{\partial S_1}a_1+ \frac{\partial V}{\partial S_2}a_2+\frac{1}{2} \frac{\partial^2 V}{\partial S_1^2}b_1^2+\frac{\partial^2 V}{\partial S_1\partial S_2}\rho b_1b_2+\frac{1}{2} \frac{\partial^2 V}{\partial S_2^2}b_2^2\right) dt + b_1\frac{\partial V}{\partial S_1}dX_1 +b_2\frac{\partial V}{\partial S_2}dX_2 $$