When the gradient and the divergence of a solenoidal function in $\mathbb{R}^d$ commute? Notation: Suppose that we are working with vector valued functions in $\mathbb{R}^d$, $d\geq1$, and the divergence is defined as $\text{div}\,\textbf{u}\equiv\partial_i\textbf{u}_i$. Moreover, divergence of a matrix is defined row-wise, i.e.,  $(\text{div}\,\textbf{U})_i=\text{div}\,\textbf{U}_{i*}$. Therefore, $(\text{div}\,\nabla\textbf{u})_i=(\Delta\textbf{u})_i=\Delta\textbf{u}_i$ is just the vector laplacian.
Question: Is there a non-singular matrix $A$ such that
$$\text{div}\,\nabla(A\textbf{v})=\nabla\text{div}\,(A\textbf{v})$$
for all $\textbf{v}$ satisfying
$$\text{div}\,\textbf{v}=0?$$
I found out that if $d=2$, then one can take $A=\left(\matrix{0&-1\\1&0}\right)$. Therefore, I need help with $d\geq3$. Should the answer be negative, I would like to know the "deep" reason for this. Thanks.
 A: Let $\mathbf w=A\mathbf v$ and $f = \operatorname{div} w$. The condition $\Delta \mathbf w = \nabla f$ requires $\Delta \mathbf w$ to be conservative. More precisely, there is a function $g$ such that $\Delta g=f$ (a solution of Poisson's equation), hence $\Delta(\mathbf w - \nabla g) = 0$. This means that 
$\mathbf w = \nabla g + \mathbf h$  where $\Delta \mathbf h=0$. Moreover, plugging such $\mathbf w$ into $\Delta \mathbf w = \nabla \operatorname{div} \mathbf w$ yields 
$$
\Delta \nabla g = \nabla \Delta g + \nabla \operatorname{div}\mathbf h
$$
where the first two terms cancel, and so $\operatorname{div} h$ must be constant.  Conclusion: the necessary and sufficient condition for $\mathbf w$ is that
$$ \mathbf w = \nabla g + \mathbf h \quad \text{where }\Delta \mathbf h=0 \text{ and } \operatorname{div}\mathbf h =  \operatorname{const}
\tag1$$
(Actually, this is equivalent to requiring $\operatorname{div}\mathbf h = 0$ because one can subtract a field of the form $\mathbf u = (cx_1, 0, \dots , 0)$ from $\mathbf h$ while  adding $cx_1^2/2$ to function $g$; this preserves equality.)
In two dimensions, rotating a divergence-free field by $90$ degrees, that is applying the matrix $A=\begin{pmatrix}0 & -1 \\ 1 & 0\end{pmatrix}$, makes it conservative. There is no such way to move from divergence-free to conservative in dimensions $d\ge 3$, and the reason is: 


*

*Divergence-free fields have $d-1$ degrees of freedom ($d$ scalar functions subject to $1$ scalar equation)

*Conservative fields have $1$ degree of freedom (it's the gradient of a scalar function)


We only have $d-1 = 1$ when $d=2$.   
And the term $\mathbf h$ in $(1)$ cannot supply the missing degrees of freedom, because the space of harmonic functions is relatively tiny (compared to all smooth functions). To be precise, recall that a harmonic function is determined by its values on any open set. So, if we have a decomposition $(1)$ for two smooth fields $\mathbf w$, and then join those functions using a smooth partition of unity, the conclusion is that the same $\mathbf h$ must work for both functions. So there are no extra degrees of freedom coming from $\mathbf h$. 
Compare with the  Helmholtz decomposition where we get, for a given field ($d$ arbitrary functions) two vector fields, one conservative ($1$ degree of freedom) and one solenoidal ($d-1$ degrees of freedom). There can also be a harmonic component there, but it's a minor part.
