How to calculate matrix divergence in polar coordinates

$defms{boldsymbol{sigma}} defs{sigma} defr{rho} deff{varphi} defo{cdot} defra{rightarrow} defve{{bf e}} defvv{{bf v}} defMM{{bf M}} defid{mathbb{I}} newcommanddv[1]{nablacdot #1} newcommandmt[4]{left[begin{array}{cc} #1 & #2 \ #3 & #4 end{array}right]} newcommandcvc[2]{left[begin{array}{cc} #1 & #2 end{array}right]} newcommandvc[2]{left[begin{array}{c} #1 \ #2 end{array}right]}$Consider $ms$ in the natural basis, $$ms = mt{s_{xx}(x,y)}{s_{xy}(x,y)}{s_{yx}(x,y)}{s_{yy}(x,y)}.$$ We can write this as $$ms = s_{ij}ve_i ve_j^T,$$ where $ve_x=cvc{1}{0}^T$ and $ve_y=cvc{0}{1}^T$. We can relate the components of $ms$ in different orthonormal bases in the following way. We have begin{align*} ms=s_{ij}ve_i ve_j^T = s_{i'j'}ve_{i'}ve_{j'}^T, tag{1} end{align*} and so $$s_{ij} = ve_i^T ve_{i'} s_{i'j'} ve_{j'}^T ve_j$$ or $$ms = MM^Tms' MM,$$ where $ms'=[s_{i'j'}]$ is the matrix of components in the primed basis and where $$M_{i'i} = ve_{i'}^T ve_i.$$ By assumption the bases are orthonormal, so $MM^T = MM^{-1}$. (Note that, by (1) and the fact that the Kronecker delta is unchanged by coordinate transformation, $id = ve_i ve_i^T = ve_{i'}ve_{i'}^T$. Thus, $[MM^TMM]_{ij} = ve_{i}^T ve_{i'} ve_{i'}^Tve_j = ve_i^T ve_j = [id]_{ij}.$)

A similar argument can be made to show that the components of a vector $vv$ in the different bases are related by begin{align*} vv' &= MM vv, tag{2} end{align*} where $vv$ and $vv'$ are the vectors whose components are in the natural and primed bases, respectively. (We assume that $MM=MM(x',y')$.)

We write the divergence in the primed basis as $[dvms]'$. Note that this quantity is a vector, and so transforms as indicated in (2). Thus, begin{align*} [dvms]' &= MM[dvms]_{(x,y)ra(x',y')} \ &= MM[dv(MM^Tms' MM)_{(x',y')ra(x,y)}]_{(x,y)ra(x',y')}. end{align*}

For this problem we have begin{align*} [dvms]' &= MM[dv(MM^Tms' MM)_{(r,f)ra(x,y)}]_{(x,y)ra(r,f)}.tag{3} end{align*} Note that begin{align*} ve_r &= cosf ,ve_x + sinf ,ve_y \ ve_f &= -sinf ,ve_x + cosf ,ve_y. end{align*} (This basis is orthonormal by inspection.) This implies, for example, that $$M_{r x} = ve_r^T ve_x = cvc{cosf}{sinf} vc{1}{0} = cosf.$$ Calculating the other components, one finds $$MM = mt{cosf}{sinf}{-sinf}{cosf}$$ or $$MM_{(r,f)ra(x,y)} = frac{1}{sqrt{x^2+y^2}} mt{x}{y}{-y}{x}.$$ Note also that $$[ms'(r,f)]_{(r,f)ra(x,y)} = ms'(sqrt{x^2+y^2},arctan y/x).$$ It is now a straightforward, if tedious, task to work out the correct form for $[dvms]'$.

Addendum

(* mma code to check (3) *) 

fs[foo_] := FullSimplify[foo, {r > 0, -Pi < f < Pi}];
rpc = {r -> Sqrt[x^2 + y^2], f -> ArcTan[x, y]};
rcp = {x -> r Cos[f], y -> r Sin[f]};
Apolar = {{srr[r, f], srf[r, f]}, {sfr[r, f], sff[r, f]}};
M = {{Cos[f], Sin[f]}, {-Sin[f], Cos[f]}};
Mt = Transpose[M];

M //. rpc // fs
Apolar[[1]][[1]] //. rpc
(* out: [a check on some relations above] *)

ansmma = Div[Apolar, {r, f}, "Polar"] // fs
(* out: [divergence according to mma] *)

ans = M.Div[Mt.Apolar.M //. rpc, {x, y}, "Cartesian"] //. rcp // fs
(* out: [divergence using (3)] *)

ans == ansmma // fs
(* out: True *)