Flux divergence formulation help

I’m trying to reformulate my PDE
\partial_L P(L,x,y,z) = \bigg(\sum_{i,j=1}^{3}a_{i,j}\frac{\partial^2}{\partial x_i\partial x_j} + \sum_{i=1}^{3}b_i\frac{\partial}{\partial x_i}\bigg)P(L,x,y,z), \quad P(L=0,x,y,z) = \delta(x)\delta(y)\delta(z-1)

where

a = \begin{bmatrix}\frac{A_3^2z^2}{2} & 0 & A_3^2xz \\ 0 & \frac12\bigg(A_1^2+A_2^2+A_3^2\frac{z^2}{x^2}\bigg) & 0 \\ A_3^2xz & 0 & \frac{A_3^2x^2}{2} \end{bmatrix}, \quad b = \begin{bmatrix}\bigg(\frac{A_3^2}{2}-A_4\bigg)x \\ -A_5 \\ \bigg(\frac{A_3^2}{2}+A_4\bigg)z\end{bmatrix}

I start by writing (which is think is correct)
so \partial_L P(L,\mathbf{x}) = (a\nabla+b)\cdot\nabla P(L,\mathbf{x})

How can I rewrite the (a\nabla+b)\cdot\nabla term in a way like \nabla\cdot Q + c\cdot \nabla where Q,c are new?