Assuming $\mathbf{S}$ is a matrix and $\mathbf{u}$ is a vector, then $\operatorname{div}\mathbf{Su}$ is a scalar. Hence, we have
\begin{align}
\operatorname{div}\mathbf{Su} = \frac{\partial}{\partial x_i}(S_{ij}u_j) = \frac{\partial S_{ij}}{\partial x_i}u_j+S_{ij}\frac{\partial u_j}{\partial x_i}.
\end{align}
Note that
\begin{align}
\frac{\partial S_{ij}}{\partial x_i}u_j=\frac{\partial S_{ij}}{\partial x_i}\mathbf{e}_j\cdot u_j\mathbf{e}_j=\mathbf{u}\cdot \operatorname{div}\mathbf{S}.
\end{align}
Next, observe
\begin{align}
S_{ij}\frac{\partial u_j}{\partial x_i}= \sum_j \sum_i S_{ij}\frac{\partial u_j}{\partial x_i} = \sum_j \sum_i S_{ji}^Tu_{j, i} = \mathbf{S}^T\cdot \operatorname{grad}\mathbf{u}.
\end{align}
Here we used the fact that
\begin{align}
\text{d}\mathbf{u}= (\operatorname{grad}\mathbf{u})^T
\end{align}
otherwise, we would have the expression
\begin{align}
\text{d}\mathbf{u}\cdot \mathbf{S}.
\end{align}
Let $\boldsymbol{F} = \left(F_1, F_2, F_3\right)$ and $\boldsymbol{G} = \left(G_1, G_2, G_3\right)$ be two vector fields.
Then, their vector product is defined as
$$
\boldsymbol{F}\times \boldsymbol{G} = (F_2G_3-F_3G_2, F_3G_1-F_1G_3, F_1G_2-F_2G_1) \Rightarrow.
$$
$$
\begin{aligned}
\text{div}\boldsymbol{F}\times \boldsymbol{G} &= \frac{\partial}{\partial x}(F_2G_3-F_3G_2) + \frac{\partial}{\partial y}(F_3G_1-F_1G_3) + \frac{\partial}{\partial z}(F_1G_2-F_2G_1) = \\
&= G_3\frac{\partial}{\partial x}F_2 + F_2\frac{\partial}{\partial x}G_3 - G_2\frac{\partial}{\partial x}F_3 - F_3\frac{\partial}{\partial x}G_2 + \\
&+ G_1\frac{\partial}{\partial y}F_3 + F_3\frac{\partial}{\partial y}G_1 - G_3\frac{\partial}{\partial y}F_1 - F_1\frac{\partial}{\partial y}G_3 + \\
&+ G_2\frac{\partial}{\partial z}F_1 + F_1\frac{\partial}{\partial z}G_2 - G_1\frac{\partial}{\partial z}F_2 - F_2\frac{\partial}{\partial z}G_1 = \\
&= G_1\left(\frac{\partial}{\partial y}F_3-\frac{\partial}{\partial z}F_2\right)+
G_2\left(\frac{\partial}{\partial z}F_1-\frac{\partial}{\partial x}F_3\right)+
G_3\left(\frac{\partial}{\partial x}F_2-\frac{\partial}{\partial y}F_1\right)- \\
&-F_1\left(\frac{\partial}{\partial y}G_3-\frac{\partial}{\partial z}G_2\right)+
F_2\left(\frac{\partial}{\partial z}G_1-\frac{\partial}{\partial x}G_3\right)+
F_3\left(\frac{\partial}{\partial x}G_2-\frac{\partial}{\partial y}G_1\right) = \\
&= \boldsymbol{G}\cdot\text{curl}\boldsymbol{F}-\boldsymbol{F}\cdot\text{curl}\boldsymbol{G},
\end{aligned}
$$
where $\text{curl}\boldsymbol{F}$ is the the curl of the vector field $\boldsymbol{F}$, and it is defined as
$$
\text{curl}\boldsymbol{F} = \left(\frac{\partial}{\partial y}F_3-\frac{\partial}{\partial z}F_2, \frac{\partial}{\partial z}F_1-\frac{\partial}{\partial x}F_3, \frac{\partial}{\partial x}F_2-\frac{\partial}{\partial y}F_1\right).
$$
Now, we have
$$
\text{div}\nabla{f}\times\nabla{g} = \nabla{g}\cdot\text{curl}(\nabla{f})-\nabla{f}\cdot\text{curl}(\nabla{g}).
$$
Further, for any scalar function $f$, by definition of the $\nabla$ and $\text{curl}$,
$$
\begin{aligned}
\text{curl}(\nabla{f}) &= \left(\frac{\partial}{\partial y}\left(\frac{\partial}{\partial z}f\right)-\frac{\partial}{\partial z}\left(\frac{\partial}{\partial y}f\right), \frac{\partial}{\partial z}\left(\frac{\partial}{\partial x}f\right)-\frac{\partial}{\partial x}\left(\frac{\partial}{\partial z}f\right), \frac{\partial}{\partial x}\left(\frac{\partial}{\partial y}f\right)-\frac{\partial}{\partial y}\left(\frac{\partial}{\partial x}f\right)\right) = (0, 0, 0) = \boldsymbol{0}.
\end{aligned}
$$
So, it is a zero vector.
Finally,
$$
\text{div}\nabla{f}\times\nabla{g} = \nabla{g}\cdot\text{curl}(\nabla{f})-\nabla{f}\cdot\text{curl}(\nabla{g}) = \nabla{g}\cdot\boldsymbol{0} - \nabla{f}\cdot\boldsymbol{0} = 0.
$$
Best Answer
Note that we have
$$\begin{align} \nabla \cdot (\phi \mathbf S)\cdot \hat x_i&=\partial_j(\phi S_{ji})\\\\ &=\phi \partial_j(S_{ji})+(\partial_j\phi)S_{ji}\\\\ &=\phi \partial_j(S_{ji})+(\partial_j\phi)S^T_{ij}\\\\ &=\phi \nabla \cdot (\mathbf S)\cdot \hat x_i+\mathbf S^T\cdot \nabla(\phi)\cdot \hat x_i \tag 1 \end{align}$$
Since $(1)$ is valid for all $i$, then
$$\nabla \cdot (\phi \mathbf S)=\phi \nabla \cdot (\mathbf S)+\mathbf S^T\cdot \nabla(\phi)$$
as was to be shown!