Why does a force field leave the momentum operator unchanged in the Schrödinger equation? The reasoning leading to the Schrödinger equation goes as follows:
A plane wave in empty space has the following form:
$$\psi = e^{i(kx-\omega t)}$$
Einstein had previously explained the photoelectric effect, i.e. the emission of electrons from a metal surface through light, by suggesting that light is made from photons containing the momentum:
$$p_\text{Photon}=\hbar k$$
and energy:
$$E_\text{Photon}=\hbar\omega$$
This was proven to be correct by experiment.
De Broglie then suggested that the same relations hold for electrons, so that, given the momentum and energy of an electron, one could find the wavelength and frequency of an electron's plane wave:
$$\begin{gather}
k_\text{Electron}=\frac{p_\text{Electron}}{\hbar}\\
\omega_\text{Electron}=\frac{E_\text{Electron}}{\hbar}
\end{gather}$$
Given the plane-wave function of an electron, one can then use certain operators to extract momentum and energy from it:
$$\begin{gather}
\hat{p}\psi = -i\hbar\partial_x\psi = \hbar k \psi\\
\hat{E}\psi = i\hbar\partial_t\psi = \hbar\omega\psi
\end{gather}$$
But for a free electron outside of a force field, the relationship between energy and momentum is given by:
$$E = \frac{p^2}{2m}$$
The Schrödinger Equation for a free electron can be derived from this equation by replacing energy and momentum by the extraction operators:
$$\hat{E}\psi = i\hbar\partial_t\psi  = \frac{\hat{p}^2}{2m}\psi=-\frac{
\hbar^2\partial^2_x}{2m}\psi$$
Schrödinger now argued that placing the free electron in a Potential would modify the equation simply through the addition of the Potential Energy:
$$i\hbar\partial_t\psi  =\biggl(-\frac{\hbar^2\partial^2_x}{2m}+V(x)\biggr)\psi$$
Now my question is: The introduction of a potential is going to mess up the plane waves very seriously, transforming them into something quite different, so the extraction operators for Energy and Momentum - which only work for plane waves, for which they were designed - are no longer going to work! How can the Schrödinger equation still hold up?
 A: While plane waves aren't eigenstates of Hamiltonians with nonconstant potentials, they're still perfectly valid states. They also form a complete set, so any state can be expressed as a sum of plane waves. The Schrödinger equation is linear, so if it applies to every plane wave, it applies to every sum of plane waves, and thus every state.
A: Think of classical mechanics. For a free particle we have a Lagrangian:
$$L=\frac{m}2 v^2.$$
Trajectories of such a particle would be just straight lines, and momentum would conserve. The momentum is defined as
$$p=\partial_{v}L=mv.$$
Energy of this particle is:
$$E=v\;\partial_vL-L=\frac{mv^2}2$$
Now, similarly, we can subtract potential energy to get Lagrangian for particle in non-uniform potential:
$$L'=\frac{m}2v^2-U.$$
Now as you say, introduction of the potential is going to mess up the trajectories very seriously, transforming them into something other. Does this mean that equation for momentum should change? In no way: $\partial_vU=0$, thus
$$p'=mv.$$
Should expression for energy change? Surely it does:
$$E'=v\;\partial_vL-L=\frac{mv^2}2+U$$
Now what is different is that $p$ is no longer a conserved quantity, so defining a single value of momentum doesn't describe motion of our particle in potential.
Quantum mechanics is quite similar. Momentum operator by definition is an operator, whose eigenstates are states with definite momentum, and corresponding eigenvalues are those definite values of momentum. Hamiltonian is by definition an operator, whose eigenstates are states with definite energy, and corresponding eigenvalues are those definite values of energy.
Now, as inhomogeneous potential makes momentum not conserve, we no longer can say that eigenstates of momentum operator remain eigenstates of energy operator. Moreover, as in classical mechanics, we can't specify single value of momentum to describe eigenstate of Hamiltonian: in state with definite total energy momentum can be different.
Another way to think of this: position operator is defined as an operator, whose eigenvalues are positions and corresponding eigenstates are states with definite position, i.e. in position representation the states are $\delta(x-x_0)$. Now, this operator would commute with energy operator only for some very weird unphysical potentials (I'm not even sure how to construct them if it's possible at all). Do we change this operator for normal potentials? Of course not, position as a physical quantity doesn't depend on potential. Yet its operator in momentum representation has the same form as momentum operator in position representation:
$$\left\langle p\right|\hat r\left| p'\right\rangle=i\hbar\nabla\delta(p-p').$$
Compare with momentum operator in position representation:
$$\left\langle r\right|\hat p\left| r'\right\rangle=-i\hbar\nabla\delta(r-r').$$
A: I found a proof that the momentum operator does not change when the wave it operates upon is not a plane wave:
Let $\phi(k) = \frac{1}{2\pi}\int\limits_{-\infty}^{\infty}e^{-ikx}\psi(x)dx$
be the Fourier transform of an (arbitrary) wave function in position space.
The expectation value of k - using the wave function $\psi(x)$ in position space - is:
$\overline{k}=\int\limits_{-\infty}^{\infty}\psi^*(x)k\psi(x)dx$
Rewriting the above expression with the Fourier Transform of $\psi(x)$ we get:
$\overline{k}=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\psi^*(x)ke^{ikx}\phi(k) dxdk = \frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\psi^*(x)ke^{ikx}e^{-ikx'}\psi(x')dxdkdx'$
We can now write $ke^{-ikx'}$ in the above expression as $i\frac{\partial}{\partial x'} e^{ikx'}$:
$\overline{k}=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\psi^*(x)e^{ikx}[i\frac{\partial}{\partial x'} e^{ikx'}]\psi(x')dxdkdx'$
Let us integrate the two last factors containing x' by parts, using the fact that the wave equation $\psi(x')$ will vanish at infinity:
$\overline{k}= \frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\psi^*(x)\frac{1}{i}\frac{\partial\psi(x')}{\partial x'}e^{ik(x-x')}dx dkdx'$
The Integral over k of the last exponential together with the factor $\frac{1}{2\pi}$is the Dirac Delta Function, so when we integrate over x, we end up with:
$\overline{k}=\int\limits_{-\infty}^{\infty}\psi^*(x')\frac{1}{i}\frac{\partial\psi(x')}{\partial x'}dx'$
Using the de Broglie equation p = $\hbar k$ and using x instead of x' for clarity we get:
$\overline{p}=\int\limits_{-\infty}^{\infty}\psi^*(x)\frac{\hbar}{i}\frac{\partial\psi(x)}{\partial x}dx$
where we see that the momentum operator $\hat{p}=-i\hbar\partial_x$
Exactly the same demonstration can be made with $\omega$ and $t$ instead of $k$ and $x$, and the result is that the Energy operator is $\hat{E}= i\hbar\partial_t$.
However, I'm still not completely convinced by this argument, because I don't really know what physical interpretation to give to the  waves $e^{ikx}$ and $e^{-i\omega t}$ in which we have Fourier-analysed the wave function $\psi (x,t)$ under consideration. We have presumed that they are plane waves with Momentum $\hbar k $ and Energy $\hbar \omega$, respectively, but $e^{ikx}$ is time independent and $e^{-i\omega t}$ is space independent - Something like this is completely unphysical, much more so than a plane wave $e^{i(kx-\omega t)}$ which might not exist for practical purposes, but could at least theoretically exist and has definite momentum AND Energy.
But I believe this is is as good as it gets.
A: The problem lies in your statement "The introduction of a potential is going to mess up the plane waves very seriously, transforming them into something quite different". The Hamiltonian operator, which is just $\hat H=\frac{\hat p^2}{2m}+V(x)$, is an operator which 'extracts' the energy.
For each operator there are states called 'eigenstates', which are very important. When operators act on these kinds of states they act as to multiply with a number. For example $e^{ipx/\hbar}$ is an eigenstate of the momentum operator: when acted on by $\hat p$, it multiplies the state by the number $p$. We say these states have definite momentum and we call this number $p$ the eigenvalue of this state. When we add two states with different eigenvalues you get a superposition. For example $e^{ip_1x/\hbar}+1/2e^{ip_2x/\hbar}$ is in a superposition. It has both momentum $p_1$ and $p_2$ at the same time. Because $p_1$ is multiplied by a higher number, there is a higher probability of observing $p_1$.
The Hamiltonian also has eigenstates, called energy eigenstates. These states have definite energy. You are correct in your suspicion that momentum eigenstates are generally not energy eigenstates and vice versa. In general, when we have an energy eigenstate it is in an infinite superposition of momentum eigenstates. There is some probability of observing $p=0.1$, some of observing $p=1.5$, some of observing $p=1000$ etc. Only when the potential is zero do these operators share eigenstates, i.e. it is possible to have definite momentum and energy at the same time.
