Confused about Brownian motion on Riemannian manifold Should I think of Riemannian Brownian motion as a stochastic process which, in some coordinate system (but not necessarily all), has the Laplace-Beltrami operator as its infinitesimal generator?

Consider the Riemannian manifold $(\mathbb{R}^n, \langle\cdot,\cdot\rangle)$ where $\langle\cdot,\cdot\rangle$ represents the standard Euclidean metric. A (global) coordinate system for the manifold is $(\mathbb{R}^n, \mathrm{Id})$.
The definition of Brownian motion on an $m$-dimensional Riemannian manifold $(M, g)$ is a diffusion process that, in local coordinates has one-half the Laplace-Beltrami operator as its infinitesimal generator. The Laplace-Beltrami operator has a representation in local coordinates $(U, \phi)$ as
\begin{align}
    \Delta f  = \frac{1}{\sqrt{\mathrm{det}(g_{ij})(x)}} \sum_{j=1}^m \frac{\partial}{\partial x_j}\left(\sqrt{\mathrm{det}(g_{ij})(x)} ~ \sum_{i=1}^m g^{ij}(x) \frac{\partial}{\partial x_i} f(\phi^{-1}(x))\right)
\end{align}
where $f\in C^\infty(M)$.
In the case of the manifold $(\mathbb{R}^m, \langle\cdot,\cdot\rangle)$, the representation of the Laplace-Beltrami operator in the local coordinate system $(\mathbb{R}^m,\mathrm{Id})$ is,
\begin{align}
    \Delta f(x) = \sum_{i=1} \frac{\partial^2}{\partial x_i^2} f(x).
\end{align}
Let $\mathrm{d}X_t^i = \mu_i(X_t)~\mathrm{d}t + \sigma_{j}^{i}(X_t)~\mathrm{d}B^{j}_t$ be a diffusion in the local coordinates $(\mathbb{R}^m,\mathrm{Id})$. The infinitesimal generator of this diffusion is,
\begin{align}
    -\sum_{i=1}^m \mu_i(x)~\frac{\partial}{\partial x_i} f(x) + \frac{1}{2}\sum_{i=1}^m\sum_{j=i+1}^m \sum_{k=1}^m \sigma^{i}_{k}(x)\sigma^{j}_{k}(x) ~\frac{\partial^2}{\partial x_i\partial x_j} f(x).
\end{align}
Thus, it is clear that the diffusion with $\mu_i(x) = 0$ for $i=1,\ldots, m$ and $\sigma_{i}^{j}(x) = \delta_{ij}$ yields a diffusion that corresponds to Brownian motion in these local coordinates. The diffusion in question is just $\mathrm{d}X_t = \mathrm{d}B_t$ and so this is perfectly aligned with our intuition that, in the chosen coordinates, Brownian motion should just be Euclidean Brownian motion.
Let $\phi :\mathbb{R}^m\to\mathbb{R}^m$ be a diffeomorphism and consider the global coordinate system $(\mathbb{R}^m, \phi^{-1})$. In these coordinates the Laplace-Beltrami operator has the form,
\begin{align}
    \Delta f(x) &= \frac{1}{\left|\mathrm{det}(\nabla \phi(x))\right|} ~\nabla \cdot \left(\left|\mathrm{det}(\nabla \phi(x))\right| ~\left[(\nabla \phi(x))^\top (\nabla \phi(x))\right]^{-1} \nabla f(\phi(x))\right) \\
    &= \frac{1}{\left|\mathrm{det}(\nabla \phi(x))\right|} \sum_{j=1}^m\frac{\partial}{\partial x_j}\left(\left|\mathrm{det}(\nabla \phi(x))\right|~ \sum_{i=1}^m g^{ij}(x)\frac{\partial}{\partial x_i} f(\phi(x))\right)
\end{align}
where
\begin{align}
g^{ij} &= \left(\left[(\nabla \phi)^\top (\nabla \phi)\right]^{-1}\right)_{ij} \\
&= \left((\nabla \phi)^{-1} ((\nabla \phi)^{-1})^\top\right)_{ij}.
\end{align}
Therefore, a diffusion whose infinitesimal generator is one-half the Laplace-Beltrami operator in these coordinates is given by,
\begin{align}
    \mathrm{d}X_t^i = \frac{1}{2\left|\mathrm{det}(\nabla \phi(X_t))\right|} \sum_{j=1}^m \frac{\partial}{\partial x_j}\left(\left|\mathrm{det}(\nabla \phi(X_t))\right| g^{ij}(X_t)\right)~\mathrm{d}t + \sum_{j=1}^m ((\nabla \phi(X_t))^{-1})^i_j~\mathrm{d}B_t^j.
\end{align}
On the other hand, if I set $Y_t = \phi^{-1}(B_t)$ then by Ito's formula I obtain,
\begin{align}
    \mathrm{d}Y^i_t &= \frac{1}{2} \sum_{j=1}^m\sum_{k=1}^m \frac{\partial^2}{\partial x_j\partial x_k} (\phi^{-1}(B_t))^i~ \mathrm{d}B_t^j~\mathrm{d}B_t^k + \sum_{j=1}^m (\nabla \phi^{-1}(B_t))_j^i ~\mathrm{d}B_t^j \\
    &= \frac{1}{2} \sum_{j=1}^m \frac{\partial^2}{\partial x_j^2}(\phi^{-1}(B_t))^i ~\mathrm{d}t + \sum_{j=1}^m ((\nabla \phi(B_t))^{-1})_j^i ~\mathrm{d}B_t^j
\end{align}
I think that $X_t$ and $Y_t$ do not represent the same diffusion. What surprises me about this computation is that the transformation of Brownian motion in the coordinate system $(\mathbb{R}^m, \mathrm{Id})$ via the diffeomorphism $\phi^{-1}$ is not Brownian motion in the coordinate system $(\mathbb{R}^m, \phi^{-1})$. Given this, I'm wondering how I should be conceptualizing Brownian motion on a manifold.
Should I think of it as a stochastic process which, in some coordinate system (but not necessarily all), has the Laplace-Beltrami operator as its infinitesimal generator?
 A: Let $(M,g)$ be a closed (compact, no boundary) Riemannian manifold with Laplace–Beltrami operator $\Delta_g$. The operator $(\Delta_g,\mathcal C^\infty(M))$ is an unbounded operator operator on $L^2(\mathrm{vol}_g)$, closable, with closure $(\Delta_g, \mathrm{dom}(\Delta_g))$. By Hille–Yoshida's Theorem, $(\Delta_g, \mathrm{dom}(\Delta_g))$ uniquely determines a strongly continuous contraction semigroup $P_t\colon L^2(\mathrm{vol}_g)\to L^2(\mathrm{vol}_g)$ defined by $P_t:= e^{t\Delta_g}$. It turns out that $P_t$ is an integral operator with kernel $p_t(x,y)$.
Definition (Brownian motion) The Brownian motion on $(M,g)$ is the unique Markov process $B$ with state space $M$ and transition probabilities $p_t(x,A)$ given by $p_t(x,A):=\int_A p_t(x,y)\,\mathrm{dvol}_g(y)$.
Note For simplicity, here we are using the analytic convention that $B$ is generated by the Laplace–Beltrami operator. The probabilistic convention is that $B$ is generated by $\frac{1}{2}\Delta_g$.
There is no substantial difference between the two.
Let me start by pointing out some problems in the understanding of the OP. The definition above makes apparent that being the generator of a Markov process is not a local property, but rather a global and intrinsic one. That is to say, the following sentence in the op

[the] diffusion process that, in local coordinates has one-half the Laplace-Beltrami operator as its infinitesimal generator

is a bit misguided.
Now to the computations:
\begin{align}
    \Delta_g f  = \frac{1}{\sqrt{\mathrm{det}(g_{ij})(x)}}\sum_{j=1}^m \frac{\partial}{\partial x_j}\left(\sqrt{\mathrm{det}(g_{ij})(x)} ~ \sum_{i=1}^m g^{ij}(x) \frac{\partial}{\partial x_i} f(x)\right)
\end{align}
is the expression for $\Delta_g$ in coordinates. (If $x\in M$, then it does not make sense to compute $f$ at $\phi^{-1}(x)$, since $f$ is defined on $M$. Here $(U,\phi)$ is a local chart with $\phi\colon U\to \mathbb R^n$ and $U\subset M$).
On the other hand, if $\phi\colon M\to M$ is a diffeomorphism with inverse $\psi:=\phi^{-1}$,

\begin{align}
  & \frac{1}{\left|\mathrm{det}(\nabla \phi(x))\right|} ~\nabla \cdot \left(\left|\mathrm{det}(\nabla \phi(x))\right| ~\left[(\nabla \phi(x))^\top (\nabla \phi(x))\right]^{-1} \nabla f(\phi(x))\right) \\
    &= \frac{1}{\left|\mathrm{det}(\nabla \phi(x))\right|} \sum_{j=1}^m\frac{\partial}{\partial x_j}\left(\left|\mathrm{det}(\nabla \phi(x))\right|~ \sum_{i=1}^m g^{ij}(x)\frac{\partial}{\partial x_i} f(\phi(x))\right)
\end{align}
where
\begin{align}
g^{ij} &= \left(\left[(\nabla \phi)^\top (\nabla \phi)\right]^{-1}\right)_{ij} \\
&= \left((\nabla \phi)^{-1} ((\nabla \phi)^{-1})^\top\right)_{ij}.
\end{align}

looks more like (up to some missing square roots) an expression in coordinates for $\Delta_{\psi^*g}$, the Laplace–Beltrami operator of the metric $\psi^*g$, i.e. the pullback metric of the standard metric $g$ on $\mathbb R^n$ via $\psi$.
Since the correspondence between Markov processes (regarded up to equivalence in law) and non-negative self-adjoint operators is one-to-one, the corresponding processes $X$ and $Y$ are the same if and only if $\Delta_g=\Delta_{\psi^*g}$. By e.g. this question, we have $\Delta_g=\Delta_{\psi^*g}$ if and only if $\psi$ is an isometry.
On $\mathbb R^n$, this means that $X=Y$ if and only if $\psi$ is a (composition of) a rotation, symmetry, or translation, in which case it is not difficult to verify that Y (when started at $0$) is a standard Brownian motion (e.g. by the multi-dimensional version of Lévy's characterization of BM).
Side note. It is certainly possible to get a local understanding of the BM on a manifold as the solution to a certain SDE on each chart of $M$. This is a done e.g. in [A] (see below). A rigorous proof of this understanding is however quite tough. One main issue appears on manifolds with more than one chart (this is not the case in the op, but I point it out for completeness): we start the process, then stop it when we hit the boundary of one chart, take another chart around the point we stopped the process at, and repeat. We need to ensure that the sequence of stopping times so constructed diverges (if not, the process explodes in finite time). Unfortunately, on general manifolds $(M,g)$ it can happen that the BM explodes in finite time, depending on $g$. The interplay between the properties of the BM and the metric
$g$ is quite subtle.
[A] N. Ikeda and S. Watanabe, Stochastic Differential Equations and Diffusion Processes
