How does the canonical momentum $p_i\equiv\frac{\partial L}{\partial\dot q_i}$ transform under a coordinates change $\mathbf q\to\mathbf Q$?

The canonical momentum is defined as

$$p_{i} = \frac {\partial L}{\partial \dot{q_{i}}},$$ where $$L$$ is the Lagrangian.

So actually how does $$p_{i}$$ transform in one coordinate system $$\textbf{q}$$ to another coordinate system $$\textbf{Q}$$?

https://en.wikipedia.org/wiki/Hamiltonian_mechanics#Charged_particle_in_an_electromagnetic_field

When dealing with the Hamiltonian of the electromagnetic field, the derivation of $$p_{j} = m \dot{x_j} + eA_j$$ on the above link is usually written as $$\textbf{p} = m \textbf{v} + e\textbf{A}$$

but the derivation is based on using Cartesian coordinates, does it mean that $$\textbf{p}$$ is really a vector? If we are using another general coordinates, say, spherical coordinates, can we still have $$\textbf{p} = m \textbf{v} + e\textbf{A}~ ?$$ If no, I think the form of Hamiltonian in electromagnetic field

$$H = \frac{(\textbf{p} - e\textbf{A})^2}{2m} + e\phi$$

will only be valid in Cartesian coordinates. In any other coordinates, $$H$$ carries a different form!

1) In curvilinear coordinates $x^i$, $i=1,2,3$, with metric $g_{ij}=g_{ij}(x)$, the non-relativistic Lagrangian for a charged particle in an electromagnetic field reads

$$L~=~T-U, \qquad T~=~\frac{m}{2}\dot{x}^i g_{ij} \dot{x}^j , \qquad U ~=~ q(\phi - \dot{x}^i A_i),$$

where $q$ is the charge. The canonical momentum is

$$p_i ~=~\frac{\partial L}{\partial \dot{x}^i} ~=~ mg_{ij} \dot{x}^j + q A_i,$$

and the Hamiltonian becomes

$$H ~=~ \frac{(p_i -q A_i) g^{ij} (p_j -q A_j)}{2m} + q\phi,$$

where $g^{ij}$ is the inverse metric.

2) These formulas behave covariantly under change of coordinates $x^i \to x^{\prime j}=x^{\prime j}(x)$ in position space. The generalized velocity

$$\dot{x}^{\prime j}~=~\frac{\partial x^{\prime j}}{\partial x^i}\dot{x}^i$$

behaves as a vector field, while the magnetic potential

$$A_i~=~A^{\prime}_j ~\frac{\partial x^{\prime j}}{\partial x^i},$$

and the canonical momentum

$$p_i~=~p^{\prime}_j ~\frac{\partial x^{\prime j}}{\partial x^i}$$

behave as covector fields (or equivalently, as one-forms). The Lagrangian $L$, the Hamiltonian $H$, and the electric potential $\phi$ are scalars (=invariant).

3) Finally, let us mention that the electric field

$$E_i~=-\frac{\partial \phi}{\partial x^i} -\frac{\partial A_i}{\partial t}$$

is a covector field, and the magnetic field

$$B^i~=\sqrt{g} ~\epsilon^{ijk}\frac{\partial A_k}{\partial x^j}, \qquad g~=~\det g_{ij}, \qquad \epsilon^{123}~=~+1,$$

is a vector field wrt. change of coordinates $x^i \to x^{\prime j}=x^{\prime j}(x)$ in position space.

The answer depends on how general a Lagrangian (or, equivalently, Hamiltonian) you want to allow. The most general Lagrangian surely allows arbitrary functions of $q,p$ and their mixing. If that's so, you shouldn't be thinking about the coordinate space for $q$ only. You have to think about the whole phase space and coordinate transformations on the whole space, i.e. reparametrizations of $(q_i,p_i)$. The phase space is a symplectic space i.e. one equipped with a non-singular antisymmetric $\omega_{ab}$ tensor, encoding the Poisson brackets of $\{X_a,X_b\}$ where $X=q$ or $p$.

In more restricted contexts, you may separate the phase space coordinates to the coordinates $q$ and momenta $p$. If you do so, $p_i$ live in a tangent space of the manifold spanned by $q_i$. The whole phase space may be interpreted as a fiber bundle over the $q_i$-generated phase space. The first identity you wrote, $p_i=\partial L / \partial \dot q_i$, gives you a relationship between $\dot q_i$ variables and $p_i$ that may also depend on $q_i$. Because the velocities $\dot q_i$ transform as tangent vectors (including the detailed values of the components) on the $q$ configuration space, you also know how $p_i$ transforms.

At any rate, I think it is a bit misleading to imagine that $p_i$ is a sort of a "field" defined in the $q$ space. In reality, $p_i$ are additional coordinates of the phase space. They may have a special relationship to the $q$'s but the transformations of the numerical values may be complex, act both on $q$ and $p$, and be affected by the Lagrangian or Hamiltonian which may be very nonlinear and complicated.