Symbol of differential operator and change of coordinates $\DeclareMathOperator{I changed index to upper one}{but edit has to take al least 6 elements}$
Some time ago I posted the question about the change of coordinates in differential operator. Here is the relevant discussion
Symbol of differential operator transforms like a cotangent vector
The answer which was given to this question explains why we can invariantly define the principial symbol of differential operator. However I would like to understand straightforward what happens with the principial symbol when we introduce new coordinates. Let me recall some formulas for transformation laws for tensors: if $(x^i)$ are old coordinates and $(y^i)$ are new coordinates then we have following formulas for canonical basis in tangent and cotangent spaces (everything takes place over manifold $M$):
$$\frac{\partial}{\partial y^i}=\sum_{k=1}^n \frac{\partial x^k}{\partial y^i} \frac{\partial}{\partial x^k}, \qquad \frac{\partial}{\partial x^i}=\sum_{k=1}^n \frac{\partial y^k}{\partial x^i} \frac{\partial}{\partial y^k},$$
$$dy^i=\sum_{k=1}^n\frac{\partial y^i}{\partial x^k} dx^k, \qquad dx^i=\sum_{k=1}^n\frac{\partial x^i}{\partial y^k} dy^k $$
If we have vector field $X=\sum_{k=1}^nX^k \frac{\partial}{\partial x^k}=\sum_{k=1}^n Y^k \frac{\partial}{\partial y^k}$ and differential one-form $\omega=\sum_{k=1}^n \alpha_k dx^k=\sum_{k=1}^n\beta_k dy^k$ expressed in these two system of coordinates then we have formulas:
$$Y^i=\sum_{k=1}^n\frac{\partial y^i}{\partial x^k}X^k, \qquad X^i=\sum_{k=1}^n\frac{\partial x^i}{\partial y^k}Y^k,$$
$$\beta_i=\sum_{k=1}^n\frac{\partial x^k}{\partial y^i}\alpha_k, \qquad \alpha_i=\sum_{k=1}^n\frac{\partial y^k}{\partial x^i}\beta_k.$$
For general tensors of type $(r,s)$ if we have in two differenet coordinate systems $$t=\sum_{i_1,...,i_r,j_1,...,j_s}t^{i_1,...,i_r}_{j_1,...,j_s}\frac{\partial}{\partial x^{i_1}} \otimes ... \otimes \frac{\partial}{\partial x^{i_r}}\otimes dx^{j_1} \otimes ... \otimes dx^{j_s}=\sum_{i_1,...,i_r,j_1,...,j_s}u^{i_1,...,i_r}_{j_1,...,j_s}\frac{\partial}{\partial y^{i_1}} \otimes ... \otimes \frac{\partial}{\partial y^{i_r}}\otimes dy^{j_1} \otimes ... \otimes dy^{j_s}$$
then we have 
$$\frac{\partial}{\partial y^{i_1}} \otimes ... \otimes \frac{\partial}{\partial y^{i_r}}\otimes dy^{j_1} \otimes ... \otimes dy^{j_s}=\sum_{k_1,...,k_r,l_1,...,l_s}\frac{\partial x^{k_1}}{\partial y^{i_1}} ... \frac{\partial x^{k_r}}{\partial y^{i_r}}\frac{\partial y^{j_1}}{\partial x^{l_1}} ... \frac{\partial y^{j_s}}{\partial x^{l_s}}  \frac{\partial}{\partial x^{k_1}} \otimes ... \otimes \frac{\partial}{\partial x^{k_r}}\otimes dx^{l_1} \otimes ... \otimes dx^{l_s},$$
$$\frac{\partial}{\partial x^{i_1}} \otimes ... \otimes \frac{\partial}{\partial x^{i_r}}\otimes dx^{j_1} \otimes ... \otimes dx^{j_s}=\sum_{k_1,...,k_r,l_1,...,l_s}\frac{\partial y^{k_1}}{\partial x^{i_1}} ... \frac{\partial y^{k_r}}{\partial x^{i_r}}\frac{\partial x^{j_1}}{\partial y^{l_1}} ... \frac{\partial x^{j_s}}{\partial y^{l_s}}  \frac{\partial}{\partial y^{k_1}} \otimes ... \otimes \frac{\partial}{\partial y^{k_r}}\otimes dy^{l_1} \otimes ... \otimes dy^{l_s}$$and also 
$$u^{i_1,...,i_r}_{j_1,...,j_s}=\sum_{k_1,...,k_r,l_1,...,l_s}\frac{\partial y^{i_1}}{\partial x^{k_1}} ... \frac{\partial y^{i_r}}{\partial x^{k_r}} \frac{\partial x^{l_1}}{\partial y^{j_1}} ... \frac{\partial x^{l_s}}{\partial y^{j_s}}t^{k_1,...,k_r}_{l_1,...,l_s},$$
$$t^{i_1,...,i_r}_{j_1,...,j_s}=\sum_{k_1,...,k_r,l_1,...,l_s}\frac{\partial x^{i_1}}{\partial y^{k_1}} ... \frac{\partial x^{i_r}}{\partial y^{k_r}} \frac{\partial y^{l_1}}{\partial x^{j_1}} ... \frac{\partial y^{l_s}}{\partial x^{j_s}}u^{k_1,...,k_r}_{l_1,...,l_s}.$$
However general differential operator contains higher order differentials for which I haven't seen the appriopriate formulas. For general differential order $m$ operator of the form $D=\sum_{|\alpha| \leq m}a_{\alpha}(x)\frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} ... \partial x_n^{\alpha_n}}$ (where $\alpha=(\alpha_1,...,\alpha_n)$ a multiindex and $a_{\alpha}(x)$ are matrices) the principial synmbol is defined as the expression $a(x,\xi)=\sum_{|\alpha|=m}a_{\alpha}(x)\xi_1^{\alpha_1} ... \xi_n^{\alpha_n}$ where $\xi=(\xi_1,...,\xi_n)$ is a vector: what we do, we replace each $\frac{\partial}{\partial x_k}$ by $\xi_k$. So far we have some formal expression $a$ which takes a point $x$ of a manifold and vector $\xi \in \mathbb{R}^n$ and produces some matrix. And here are my questions:

  
*
  
*What transformation law the function $a$ must obey in order to correctly define a mapping from the cotangent bundle $T^*M$? Or maybe we look for some transformation laws for the variables $\xi_k$?    
  
*How to show that such transformation law holds?  
  

I repeat once again: I already understood how to define symbol globally as a map from cotangent space and that this definition gives exactly the local expression described above-however I would like to understand the reason for such definition which at the moment is pulled out of hat.Concerning the second question I would like to understand whether there are some general formulas for transformation laws for an arbitrary differential operators: as far as I remeber, when I once saw the computations for spherical Laplacian, it was not clear whether such method will work in general arbitrary example.
 A: In your tensor transoformation formulas you gracefully use lower and upper indices, but when you start talking about linear differential operators (LDO) you write it strangely as $a_\alpha\frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} ... \partial x_n^{\alpha_n}}.$ 
Notice that any vector field $X$ is a LDO of order 1 and locally can be described, as you wrote above, as $X^i\frac{\partial}{\partial x^i}.$ Hence LDOs are generalizations (almost compositions) of vector fields.
So let us use the convention to write LDOs in map as $a^\alpha\frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} ... \partial x_n^{\alpha_n}}.$ Now it looks like vector field case, i.e $a^\alpha=X^i.$
Notice that we write $\frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} ... \partial x_n^{\alpha_n}}$ with indices on bottom. However we still apply functions defined on $M$ to this operator. So it rather should be denoted as
$$\frac{\partial^{|\alpha|}}{\partial (x^1)^{\alpha_1} ... \partial (x^n)^{\alpha_n}}$$
But convention is different, I guess due to lack of space above $x$'s.
For simplicity let's say that we have LDO $\mathcal{D}:C^\infty(M)\to C^\infty(M).$ By this constrain we get rid of matrices $a^\alpha$ and simply we have functions instead. Assume that $\mathcal{D}$ is of order $m.$ Hence locally $\mathcal{D}$ takes a form
$$D=\sum_{\alpha\leqslant m}a^\alpha\frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} ... \partial x_n^{\alpha_n}}$$
At this point you should complain that we still do not know how does $a^\alpha$'s transform. 
However $a^\alpha$'s transformation rule depends on the order $|\alpha|.$ To make long story short: transformation of $k$ order part influences parts of $k,k-1,\dots, 1$ orders. Leibniz rule is guilty for such behavior.
For example consider 3rd order LDO $\mathcal{D},$ which in some map has only $3$ and $0$ order parts, i.e $\mathcal{D}$ locally takes form
$$Df=\sum_{i_i,i_2i_3}a^{i^1i^2i^3}\frac{\partial^3f}{\partial x^{i_1}\partial x^{i_2}\partial x^{i_3}}+fh.$$
Recall that $0$ part is just a multiplication by other function (in our case $h$). 
Now we change coordinates from $(x^1,\dots,x^n)$ to $(y^1,\dots,y^n).$ We use transformation rules that you wrote above and we see that in new coordinates $\mathcal{D}$ takes form
$$Df=\sum_{i_i,i_2i_3}a^{i_1i_2i_3}\sum_{j_1}\frac{\partial y^{j_1}}{\partial x^{i_1}}\frac{\partial}{\partial y^{j_1}}\left(\sum_{j_2}\frac{\partial y^{j_2}}{\partial x^{i_2}}\frac{\partial}{\partial y^{j_2}}\left( \sum_{j_3}\frac{\partial y^{j_3}}{\partial x^{i_3}}\frac{\partial f}{\partial y^{j_3}} \right)\right)+fh.$$
As I mentioned now we will use only Leibniz rule. Hence 
$$Df=\sum_{j_1,j_2,j_3}\sum_{i_i,i_2i_3}a^{i_1i_2i_3}\frac{\partial y^{j_1}}{\partial x^{i_1}}\frac{\partial}{\partial y^{j_1}}\left(\frac{\partial y^{j_2}}{\partial x^{i_2}}\left(\overbrace{\frac{\partial}{\partial y^{j_2}}\left(\frac{\partial y^{j_3}}{\partial x^{i_3}}\right) \frac{\partial f}{\partial y^{j_3}}}^{(1)} +\overbrace{\frac{\partial y^{j_3}}{\partial x^{i_3}}\frac{\partial^{i_2i_3}f}{\partial y^{i_2}y^{i_3}}}^{(2)}\right)\right)+fh.$$
Again we have to use Leibniz rule, but now to 3 factors. But I guess you can see that $(1)$ would produce $2$nd and $1$st order parts and $(2)$ would produce $3$rd and $2$nd order parts.
So in different coordinates we produced orders of $2$ and $1$ as well, but haven't touched $0$ order part.
It is sufficient to see what have happened to the order $3$ part. Keeping track of $3$rd order part in the above, we see that it takes form
$$\sum_{j_1,j_2,j_3}\sum_{i_i,i_2i_3}a^{i_1i_2i_3}\frac{\partial y^{j_1}}{\partial x^{i_1}}\frac{\partial y^{j_2}}{\partial x^{i_2}}\frac{\partial y^{j_3}}{\partial x^{i_3}}\frac{\partial^{i_1i_2i_3}f}{\partial y^{j_1}y^{j_2}y^{j_3}}.$$
Hence we get that $$a^{j_1j_2j_3}=\sum_{i_i,i_2i_3}a^{i_1i_2i_3}\frac{\partial y^{j_1}}{\partial x^{i_1}}\frac{\partial y^{j_2}}{\partial x^{i_2}}\frac{\partial y^{j_3}}{\partial x^{i_3}}$$
In order to make principal part $C^\infty$ linear we have to kill $\frac{\partial y^{j_1}}{\partial x^{i_1}}\frac{\partial y^{j_2}}{\partial x^{i_2}}\frac{\partial y^{j_3}}{\partial x^{i_3}}$ somehow. So for covector $\alpha$ we define $\sigma(x,\alpha):C^\infty(M) \to C^\infty(M)$ in $(x^1,\dots,x^n)$ coordinates as 
$$\sigma(x,\xi_i dx^i)f=\xi_{i_1}\xi_{i_2}\xi_{i_3}a^{i_1i_2i_3}f.$$
(From now I will skip sums and use Einstein notation) By your formula above we know that $\xi_j=\xi_i\frac{\partial x^i}{\partial y^j}$ hence principal symbol in $(y^1,\dots,y^n)$ coordinates takes form
$$\sigma(x,\xi_jdy^j)f=\xi_{j_1}\xi_{j_2}\xi_{j_3}a^{j_1j_2j_3}f=\xi_{i_1}\xi_{i_2}\xi_{i_3}\frac{\partial x^{i_1}}{\partial y^{j_1}}\frac{\partial x^{i_2}}{\partial y^{j_2}}\frac{\partial x^{i_3}}{\partial y^{j_3}}a^{i'_1i'_2i'_3}\frac{\partial y^{j_1}}{\partial x^{i'_1}}\frac{\partial y^{j_2}}{\partial x^{i'_2}}\frac{\partial y^{j_3}}{\partial x^{i'_3}}=\\ \xi_{i_1}\xi_{i_2}\xi_{i_3}a^{i_1i_2i_3}f=\sigma(x,\xi_i dx^i)f.$$
Hence $\sigma(x,\alpha)$ is well defined.
BTW. To see how general transformation rules are sick, let's think about LDO as a composition of vector fields $\mathcal{D} = X_1\circ\dots\circ X_m.$ Now try compute $X_1\circ\dots\circ gX_k\circ\dots\circ X_m.$ Notice that it will produce a lot of summands. While changing coordinates you have to multiply by such $g_k$ each $X_k.$ This is why it is hard.
