Multiple choices of gauge for a single many-body wavefunction in non-relativistic quantum mechanics It is known that in a single particle quantum mechanics problem with the Hamiltonian, $H = \frac{(\vec p-q\vec A)^2}{2m} + V(\vec r)$, one can perform the following gauge transformation:
$$\vec A \rightarrow \vec A' = \vec A + \vec \nabla \lambda(\vec r),$$
provided we also transform the wave-function, $\Psi(\vec r) \rightarrow \Psi'(\vec r) = e^{iq\lambda(\vec r)/\hbar} \Psi(\vec r) $, so that the Schr$\ddot o$dinger equation remains satisfied.
However, if we have multiple particles with the Hamiltonian being,
$$H = \sum \frac{(\vec p_n-q\vec A_n)^2}{2m} + V(\vec r_n),$$
is one allowed to choose the gauge for each term $\vec A_n$ differently for each particle?
Also, would the corresponding wave-function after gauge transformation look like $$\Psi'(\vec r_1, ..., \vec r_N) = e^{\frac{iq}{\hbar}\sum \lambda_n(\vec r_n) } \Psi(\vec r_1, ..., \vec r_N)~?$$
I tried to verify if this is true and it seems to be so. In that case, this is very strange since a single source of magnetic field will give rise to as many gauge choices as there are particles.
 A: In electromagnetism there is one vector potential, $\vec{A}$. So the correct way to write the Hamiltonian is
\begin{equation}
H = \sum_n \frac{(p_n-q_n \vec{A})^2}{2m_n} + V(\vec{r}_1, \vec{r}_2,...\vec{r}_N)
\end{equation}
This Hamiltonian is invariant under gauge transformations of the form (apologies if I get a sign wrong)
\begin{eqnarray}
\Psi_n &\rightarrow& e^{i q_n \lambda(\vec{r}_n)/\hbar} \Psi_n \\
\vec{A} &\rightarrow& \vec{A} + \nabla \lambda
\end{eqnarray}
ie, the vector potential transforms the same way each place it appears, and all fields transform with the same $\lambda(\vec{r}_n)$, with the exact phase factor determined by the charge $q_n$. It is not consistent to allow the same field ($\vec{A}$) to transform in different ways if it appears in multiple places in the Hamiltonian.
(Note: in an earlier version I wrote $x$ instead of $\vec{r}_n$ as the argument to $\lambda$, but changed this due to a comment)
Aside
Based on the comments, what follows from this point is not the situation the OP is interested in. But, for completeness, I'll note we can consider a situation with $N$ vector potentials (where $N$ is the number of particles). We can represent this situation with a Hamiltonian
\begin{equation}
H = \sum_n \frac{(p_n-q_n \vec{A}_n)^2}{2m_n} + V(\vec{r}_1, \vec{r}_2,...\vec{r}_N)
\end{equation}
which would be invariant under $N$ different gauge symmetries (the group is then $U(1)^N$). In this situation, we can transform each field separately. Let me use $j$ to label the fields, to try to make it obvious there is a difference with the situation above, in that here only two fields are transforming, $\Psi_j$ and $\vec{A}_j$, as opposed to $N+1$ in the case of electromagnetism.
\begin{eqnarray}
\Psi_j &\rightarrow& e^{i q_j \lambda_j(\vec{r}_j)/\hbar} \Psi_j \\
\vec{A}_j &\rightarrow& \vec{A}_j + \nabla \lambda_j
\end{eqnarray}
In this situation, since there are $N$ different vector potentials, it is consistent for them each to transform in a different way.
A: You have a single electromagnetic field described by single
potentials $\vec{A}(\vec{r},t)$ and $V(\vec{r},t)$.
These potentials extend over the whole space (all $\vec{r}$).
And you have $N$ particles
represented by the wave function $\Psi(\vec{r}_1,\dots,\vec{r}_N,t)$.
Then the Hamiltonian is
$$H=\sum_{n=1}^N\left
(\frac{1}{2m}(\vec{p}_n-q\vec{A}(\vec{r}_n,t))^2
+qV(\vec{r}_n,t)\right).$$
Notice that you still have only one electromagnetic field ($\vec{A}$ and $V$).
However, only the field values at positions $\vec{r}_n$ are relevant,
because the electromagnetic interaction takes place only where the
particles are.
Then Schrödinger's equation is
$$\sum_{n=1}^N\left
(\frac{1}{2m}(\vec{\nabla}_n-q\vec{A}(\vec{r}_n,t))^2
+qV(\vec{r}_n,t)\right)\Psi(\vec{r}_1,\dots,\vec{r}_N,t)
=i\hbar\frac{\partial}{\partial t}\Psi(\vec{r}_1,\dots,\vec{r}_N,t)$$
where $\vec{\nabla}_n$ is with respect to $\vec{r}_n$.
It is easy to verify that this Schrödinger equation is invariant
to the following gauge transformation:
$$\begin{align}
\vec{A}'(\vec{r},t)&=\vec{A}(\vec{r},t)+\vec{\nabla}\lambda(\vec{r},t) \\
V'(\vec{r},t)&=V(\vec{r},t)-\frac{\partial}{\partial t}\lambda(\vec{r},t) \\
\Psi'(\vec{r}_1,\dots,\vec{r}_N,t)&=
  e^{iq\lambda(\vec{r}_1,t)/\hbar}...e^{iq\lambda(\vec{r}_N,t)/\hbar}
  \Psi(\vec{r}_1,\dots,\vec{r}_N,t)
\end{align}$$
Notice that you have one common function $\lambda$.
Otherwise gauge invariance would not hold.
