The Momentum Operator in QM I've seen the 'derivation' as to why momentum is an operator, but I still don't buy it. Momentum has always been just a product $m{\bf v}$. Why should it now be an operator. Why can't we just multiply the wave function by $\hbar{\bf k}$? Why should momentum be a derivative of a wave function?
 A: The simplest possible solution to the Schrodinger equation is that of a plane wave: $$ \psi = e^{i(\boldsymbol{\mathbf{k}}\cdot \boldsymbol{\mathbf{r}}- \omega t )} $$ In this scenario, the position of the system is indeterminate, but the momentum is known.  Next, take the gradient of this wavefunction. $$ \nabla \psi = ik \psi$$ But we know that $p = \hbar k $, from De Broglie.  Substituting this into the equation and rearranging  yields $$ \frac{\hbar}{i} \nabla \psi = p \psi $$ This is an example of an eigenvalue equation - an operator acts on a function or vector to produce a constant multiple of the original function.   So, we see that the eigenvalue is the measured momentum.  Therefore, the operator producing it, $ \frac{\hbar}{i} \nabla $, is called  the momentum operator.
A: Consider a theory with certain dynamical objects of interest. This can be the particle trajectory $x(t)$ in classical mechanics, or the wave fucntion $\psi(x,t)$ in the quantum theory.
By Noether's theorem, if certain conditions on the dynamics of these objects hold, then an infinitesimal pertubation which is characterized by derivatives $Q$ of the object involved in the transformation, gives you a conserved quantity $I$. The quantity which is conserved because of time invariance is what you call "energy", that which is conserved because of translation invariance is what you call "momentum" and so on. To see why the spatial derivative pops up when talking about the momentum is to see how fields come together with translations. 
I don't know your math background, but essentially "$I\propto Q$". Still roughly speaking, but in more technical terms, you obtain $I$ by multiplying by the conjugate momentum, which is different for every theory and then integrate over. This section of the wikipedia page on Noether's theorem suggests how in classical mechanics, where the dynamical object of interest is the trajectory $x(t)$, the different conserved observables are computed from the respective transformations. Now the example of translation in classical mechanics isn't particularly illuminating for our purposes here, because the momentum happens to coincide with the conjugate momentum. But in the following case, you can for example see how the cross product in the angular momentum arises from the cross product involved in the infinitesimal change if you do a rotation.
Now if you work with a field $f$, then the Taylor series at the point $x$ is given by 
$f(y) = \sum_{n=0}^{\infty} \frac{1}{n!}\frac{\partial^n f(x)}{\partial x^n } (y-x)^n.$
Set $y=x-\lambda$ and you get
$f(x-\lambda)=\sum_{n=0}^{\infty} \frac{1}{n!} \frac{\partial^n f(x)}{\partial x^n }  ((x-\lambda)-x)^k=\sum_{n=0}^{\infty} \frac{1}{n!}\left(-\lambda \frac{\partial }{\partial x }\right)^n f(x),$
which expresses translation of the field as an operator. The Taylor expansion uses the derivatives at a point to reconstruct the values of the function at other points and this is why $\frac{\partial }{\partial x }$'s suffice to describe translations. Here a little demonstration on the polynomial $p(x):=2x-bx^3$, the second and the last line giving the same result: http://i.imgur.com/VuWbVg3.png
Physicists also like to write the above formula as 
$f(x-\lambda)=\text{e}^{-i \lambda\ (-i\frac{\partial}{\partial x})}f(x),$
because $-i\frac{\partial}{\partial x}$ is hermitean, i.e. mathematically symmetric in the space where physical states live in. 
The quantity which is conserved because of translation invariance (which we like to call momentum), is given by the infinitesimal version of this, so "$I\propto Q\propto \frac{\partial}{\partial x}f(x)$". 
The above is also true for a classical field theory. If you want to apply this to QM now, consider the Schrödinger theory with $f=\psi$. The canonical momentum is given as $-i\psi$, giving the conserved quantity
$I=\int (-i\psi(x))\cdot \frac{\partial}{\partial x}\psi(x),$
which you can read as the expectation value $\int \psi(x)\left(-i\frac{\partial}{\partial x}\right)\psi(x)$.
A: As opposed to classical mechanics, where a particle is represented by a momentum and position vector that have determined directions and magnitudes, in Q.M. the particle is in a "superposition" of different positions and momentums. So it has a whole bunch of positions and momentums, all at the same time. As Slaviks pointed out, multiplying by a specific momentum value would no longer work, as you would have to choose one of the several (generally very many, mostly infinite) possible momenta. In the  cases where a particle is in a pure momentum state, where it has only one momentum, you can indeed just as well multiply by the momentum instead of applying the momentum operator.
The momentum operator is useful for finding the average momentum should you repeatedly measure the momentum of particles that each have an identical wave function:
$$\overline{P}=<\psi|-i\hslash\partial_x|\psi>= \\ \sum\limits_{All \space momentum \space states \space \psi_i \\ with \space momentum \space P_i \space that \space \psi \space is\space composed\space of}<\psi_i|P_i|\psi_i> = \sum\limits_{All\space momenta \space P_i \space contained \space in \space \psi }(Probability\space of \space finding \space P_i)\cdotp(P_i)$$
Each individual measurement will not give the average value, but one specific momentum value. That's because measurement has caused the "collapse of the wave function" culling all of the $\psi_i$ except one. When you measure that specific particle's momentum again, you will again find the same momentum. It is now in a pure momentum state, for example $\psi_{123}$. 
Operating on a mixed-momentum wave function with the momentum operator without bra-keting (averaging over momenta) isn't something very useful, because it will give you a new wave function than the one you started with (not corresponding to a state that you are dealing with), multiplied by a quantity with the unit of a momentum.
A: Theoretical physics is about constructing a mathematical model which we hope describes the phenomena it's being modeled for and hence helps predicting stuff. In classical physics 
this mathematical model is based simply on the real numbers (at least locally) because of the nice behaviour of things. In quantum mechanics it's not the case. Experiments started giving discrete values as well as continuous values(such as energy of electrons,etc.). So there is a need for some class of mathematical objects that give equal importance to continuous and discrete cases.
linear operators have the property that they posses both discrete and continuous spectra and hence can act as the required class of mathematical objects. Hence we start identifying physical observables with appropriate linear operators.
Postulate 1. To each dynamic variable there exists a linear operator such that possible values are the eigenvalues of the operator.
We need some place where all the physics happens and where these operators act to give us the required results. So we construct a Hilbert space consisting of states of the system which we are observing.
Theorem (Wigner). Any transformation of the vector space onto itself that preserves the inner product must be unitary or anti unitary.
Hence to each symmetry transformation we have a unitary representation of it. The symmetries of nature are assumed to be the Galilean symmetries (Non relativistic). These are rotations, displacements, and transformations between uniformly moving frames of reference. Let $-J_\alpha, -P_\alpha,G_\alpha,H$ be the corresponding generators then the unitary transformation to the corresponding symmetry is,
$$U(s)=e^{isK}\approx I+isK$$
$K$ is one of the generators. If $Q$ represents the position operator,
$$Q|\Psi\rangle=x|\Psi\rangle$$
Write $\langle x|\Psi\rangle =\Psi(x)$. We want the momentum operator to be such that,
$$e^{ia\cdot P/\hbar}|x\rangle=|x+a\rangle$$
$$\implies \Psi(x+a)=\langle x+a|\Psi\rangle=\langle x|e^{ia\cdot P/\hbar}|\Psi\rangle\approx \left\langle x \bigg|1+\dfrac{i}{\hbar}a\cdot P\bigg|\Psi\right\rangle$$
$$=\Psi(x) +\left[\dfrac{i}{\hbar}a\cdot P\right]\Psi(x)$$
Don't worry about $\hbar$ it's just a conversion factor. We used $\langle x|\Psi\rangle=\Psi(x)$. Expanding $\Psi(x+a)$ in taylor series expansion ( we can do this because we know the wave function is a nice function from the continuity equation or the probability flux equation) we get,
$$\Psi(x+a)\approx\Psi(x)+a \frac{\partial \Psi(x)}{\partial x}$$
Comparing the two we get,
$$P_x=-i\hbar \frac{\partial }{\partial x}$$
