dg.differential geometry - What is the symbol of a differential operator?

I find Wikipedia's discussion of symbols of differential operators a bit impenetrable, and Google doesn't seem to turn up useful links, so I'm hoping someone can point me to a more pedantic discussion.

Background

I think I understand the basic idea on $mathbb{R}^n$, so for readers who know as little as I do, I will provide some ideas. Any differential operator on $mathbb{R}^n$ is (uniquely) of the form $sum p_{i_1,dotsc,i_k}(x)frac{partial^k}{partial x_{i_1}dotspartial x_{i_k}}$, where $x_1,dotsc,x_n$ are the canonical coordinate functions on $mathbb{R}^n$, the $p_{i_1,dotsc,i_k}(x)$ are smooth functions, and the sum ranges over (finitely many) possible indexes (of varying length). Then the symbol of such an operator is $sum p_{i_1,dotsc,i_k}(x)xi^{i_1}dotsoxi^{i_k}$, where $xi^1,dotsc,xi^n$ are new variables; the symbol is a polynomial in the variables ${xi^1,dotsc,xi^n}$ with coefficients in the algebra of smooth functions on $mathbb{R}^n$.

Ok, great. So symbols are well-defined for $mathbb{R}^n$. But most spaces are not $mathbb{R}^n$ — most spaces are formed by gluing together copies of (open sets in) $mathbb{R}^n$ along smooth maps. So what happens to symbols under changes of coordinates? An affine change of coordinates is a map $y_j(x)=a_j+sum_jY_j^ix_i$ for some vector $(a_1,dotsc,a_n)$ and some invertible matrix $Y$. It's straightforward to describe how the differential operators change under such a transformation, and thus how their symbols transform. In fact, you can forget about the fact that indices range $1,dotsc,n$, and think of them as keeping track of tensor contraction; then everything transforms as tensors under affine coordinate changes, e.g. the variables $xi^i$ transform as coordinates on the cotangent bundle.

On the other hand, consider the operator $D = frac{partial^2}{partial x^2}$ on $mathbb{R}$, with symbol $xi^2$; and consider the change of coordinates $y = f(x)$. By the chain rule, the operator $D$ transforms to $(f'(y))^2frac{partial^2}{partial y^2} + f''(y) frac{partial}{partial y}$, with symbol $(f'(y))^2psi^2 + f''(y)psi$. In particular, the symbol did not transform as a function on the cotangent space. Which is to say that I don't actually understand where the symbol of a differential operator lives in a coordinate-free way.

Why I care

One reason I care is because I'm interested in quantum mechanics. If the symbol of a differential operator on a space $X$ were canonically a function on the cotangent space $T^ast X$, then the inverse of this Symbol map would determine a "quantization" of the functions on $T^ast X$, corresponding to the QP quantization of $mathbb{R}^n$.

But the main reason I was thinking about this is from Lie algebras. I'd like to understand the following proof of the PBW theorem:

Let $L$ be a Lie algebra over $mathbb{R}$ or $mathbb{C}$, $G$ a group integrating the Lie algebra, $mathrm{U}L$ the universal enveloping algebra of $L$ and $mathrm{S}L$ the symmetric algebra of the vector space $L$. Then $mathrm{U}L$ is naturally the space of left-invariant differential operators on $G$, and $mathrm{S}L$ is naturally the space of symbols of left-invariant differential operators on $G$. Thus the map Symbol defines a canonical vector-space (and in fact coalgebra) isomorphism $mathrm{U}Ltomathrm{S}L$.