why does Lie bracket of two coordinate vector fields always vanish?
This is really puzzling me. Say we are dealing with a Riemannian manifold $(M,g)$. Suppose $\nabla$ is the unique torsion free connection on $M$ that is compatible with $g$. Suppose we are in a neighbourhood $U$ with coordinate map $(x^1,\cdots, x^m )$. Since the connection is torsion free, $$[\frac{\partial}{\partial x_i},\frac{\partial}{\partial x_j}]=\nabla_{\frac{\partial}{\partial x_i}}{\frac{\partial}{\partial x_j}}-\nabla_{\frac{\partial}{\partial x_j}}{\frac{\partial}{\partial x_i}}.$$
And since the $\Gamma_{i,j}^k$ is symmetric on $i,j$, the right hand side of the above equation will vanish. So the Lie bracket will be 0. Now here is my confusion. If I start out with $m$ linearly independent vector fields $Y_1, \cdots, Y_m$, then I can find a coordinate system $(y_1,\cdots, y_m)$ such that $Y_i = \frac{\partial }{\partial y_i}$ (Correct me if I am wrong, because I am not sure about this). Then arguing as above, I can show that the Lie bracket of $Y_i$ and $Y_j$ vanishes. I know Lie bracket shouldn't vanish on any two random vector fields I pick. So there must be something wrong with my argument here. Thank you in advance!
Ultimately, the vanishing is due the fact that partial derivatives commute. The reason we cannot realize a given set of vector fields as coordinate derivations is ultimately due to non-trivial curvature on the manifold.
On the other hand, if a given set of nontrivial vector fields have vanishing brackets, or more generally Lie brackets which close on the span of the vector fields then we say such a set of vector fields is involutive. This means there exists a submanifold of the given manifold which takes the given set of vector fields as tangents. See Frobenius Theorem or this related MSE question for example.
I should also mention, it is always possible to take one nontrivial vector field and make it a coordinate derivation see straightening theorem (which is a trivial case of the more general result of Frobenius).
What you write before "Correct me if I'm wrong" is indeed wrong.
In $\mathbb{R}^n$ we know that for any smooth function $f$ we have$$\frac{\partial^2f}{\partial x_1\partial x_2}=\frac{\partial^2f}{\partial x_2\partial x_1}.$$In other words:$$\left[\frac{\partial}{\partial x_1},\frac{\partial}{\partial x_2}\right]=0.$$Hence, whenever $X_1,\ldots,X_n$, are local vector fields on a manifold induced by a parametrization, we necessarily have $[X_i,X_j]=0$ for $i,j=1,\ldots,n$.