How does one describe a state with a density matrix after measuring position? My question is about position measurement in non relativistic quantum mechanics. I've been taught that when you measure the value of an observable for some state of a system described by $|\psi\rangle$ then the state of the sistem "collapses" to the eigenvector associated with the eigenvalue measured. 
Following this logic, when you measure the position of a particle described by $|\psi\rangle$ and you get a value $x_0$ then the state should collapse to an eigenvector $|x_0\rangle$. However, as is obvious, this eigenvector is the dirac delta function in position basis, $\langle x|x_0\rangle = \delta (x-x_0)$. This seems to me unacceptable since this state is non normalizable and furthermore simple things, such as $\langle \hat{x}\rangle$ are undefined for this state.
The first thing I tried when I came up with this is to try to describe the state with a density matrix $\hat{\rho} = \sum_{i} p_i \vert \psi \rangle \langle \psi \vert$ but I couldn't manage to find a way in which to do this satisfactorily. On the one hand, if I try to do this in the obvious way I get something like $\mathrm{Tr}(\hat{\rho} \hat{x})=\sum_i p_i \langle \hat{x} \rangle_{x_i}$ which is not helpful because I've got no way to individually evaluate the $\langle \hat{x} \rangle_{x_i}$ terms. On the other hand I think I would need some way to assign a probability distribution to the different eigenstates, it seems to me that in a real world scenario you would assign, for example, a gaussian centered in some value $x_0$ to the value you measure (I apologize if this seems vague but I hope the point comes across) and not probabilties to the individual eigenstates in a discrete way, as is usual when one defines the density matrix operator.
In short, how is this usually done in practice? I believe there must be a straightforeward way of going about this that I'm not seeing, or a generally agreed upon manner to handle it. I'm sure experimentalists make measurements of position all the time, and I guess its only natural to describe the state after measurement in some way. How is this done? It bugs me that this seems such a simple, basic question, yet I haven't been able to find the answer anywhere.
Note: I've been taught the Copenhagen interpretation (or wharever you want to call the fact that the state "collapses" after measurement) and understand quantum mechanics based on this. I would really appreciate answers to stick with this idea unless for some reason its unavoidable to talk about other interpretations in this context. My guess is that, if its true that all interpretations (or most) are actually "equal" in the sense that they have the same underlying mathematics and give the same answer then this question must have some satisfactory answer in the context of this interpretation. If this is not the case for some reason please explain why, and please take into account that I'm not really familiar with other interpretations of QM. Thanks in advance!!
 A: 
How do you describe a state with a density matrix after measuring position?

If you already have a density matrix, you use the Heisenberg picture.  In the Heisenberg picture the state never changes. Instead, the operators corresponding to observables change in time. So if you measure position at time $t_1$, then density matrix and state does not change, but now if later you want to measure anything else, such as position at time $t_2$ then you need to consider operators such as $\hat x(t_2)-\hat x(t_1)$ depending on exactly what you are measuring.
However, maybe your position operator was part of your state preparation. So you measured your system at time $t_0$ and got a result $x_0$ and all of that was state preparation. In that case your density matrix is the pure state $|x_0\rangle\langle x_0 |$. If measuring was part of the preparation, but the result $x_0$ was not, then if $\Psi$ is the raw material state, then you might consider the mixed state $\displaystyle \int\mathrm dx\,\overline{\Psi(x)}\Psi(x)|x\rangle\langle x |.$
OK, those are generally the choices for using a density matrix, leave the state alone and adjust the operators corresponding to observables over time (the usual approach) or to make a pure or mixed state depending exactly on how you prepared your original state. But the body of your question seemed to be concerned with the eigenstates of the position operator, so let's deal with that.
One approach is to use Rigged Hilbert Spaces or something similar so that eigenfunctions of position and momentum are mathematically legitimate things.  This is an entirely additional layer of abstraction, the bra and ket become functions themselves so the density matrix is an even higher order function.  So maybe making sure you know what a density matrix is backwards and forwards is helpful to do first to avoid getting lost.
Another approach is to look at the actual experimental setup where the measurement is not perfect, instead you have a detector that reacts over a finite region, so the eigenstate is the post interaction state, which has support over a finite region.  Make sure to realize that a physically realizable measurement might not be complete in the sense that there might be other operators that commute with it, so don't make the mathematical error of assuming it is complete.
You can also do lattice approaches where continuous space is replaced with a discrete lattice of points.  Or replace infinite space with a giant box of length $L$ and take the limit as $L$ goes to infinity later.  Or any other approach that replaces the mathematical fundamentals but is designed to give the same answer as what we see in reality when you take the final limit.
The philosophy behind all these approaches is that getting a perfect eigenfunction of momentum or position isn't something that you really do, it is an approximation of what you really do, an idealization.  So you can either desrcribe what you really do, or you can adjust the system in ways that don't affect how the idealization differs from what you really do. You have some freedom to adjust the math and the model to accommodate things that you don't really do, to make them convenient for you.  Convenience for you is really the only legitimate reason to talk about things you don't (and can't) actually do.
A: You can use the vector, just expressing it at a "lower level" is difficult
When you say $\langle x| x_0\rangle = \delta(x - x_0)$ you are already succeeding, even if you're a little confused about what the "square root of a $\delta$-function" looks like in terms of a real, honest-to-goodness wavefunction. So for example once you have the above result you can easily write $$\hat x = \int_{\mathbb R} \mathrm dx'~ x' ~ |x'\rangle\langle x'|,$$ whence you discover that $$\hat x |x_0\rangle = x_0 |x_0\rangle,$$ hence if $\rho_0 = |x_0\rangle\langle x_0|$ then $\langle \hat x \rangle_{\rho_0} = x_0$, and if $\rho = |\Psi\rangle\langle\Psi|$ for $|\Psi\rangle = \displaystyle \int_{\mathbb R} \mathrm dx~\Psi(x)|x\rangle$ then $$\langle \hat x \rangle_\rho = \iiint_{\mathbb R^3} ~\mathrm dx~\mathrm dx'~\mathrm dx''~\langle x'' | \Psi^*(x'') ~ x' |x'\rangle\langle x'| ~ \Psi(x) |x\rangle = \int_{\mathbb R} \mathrm dx~\Psi^*(x) ~ x ~ \Psi(x).$$ Doing the analogue for $\hat p$ also seems feasible but you'd write it by writing out the Fourier transform to and from the basis rather than acting on the parameter $x$.
So there is nothing wrong with this; you are merely, I think, uncomfortable with writing $\Psi_{x_0}(x) = \sqrt{\delta(x - x_0)}$ or so, because the Dirac $\delta$-function is not really a function. You should feel uncomfortable about that. Let's talk about how to rectify it. (Disclaimer: I'm stealing half an answer from someone I already replied to.)
How to do continuous systems
For continuous systems, we want a family of solutions based on some parameters which I'll collectively identify as $\alpha \in A$; the solutions are then labeled simply as $|\alpha\rangle$. They do not need to be eigenfunctions of a Hamiltonian or orthogonal or any such thing to use them as a "basis", they just need to satisfy one major requirement: there must be a normalization kernel, a function $\mathcal N(\alpha)$, such that they, with the kernel, "resolve the identity":$$\hat 1 = \int_{A} d\alpha ~\mathcal N(\alpha)~ |\alpha\rangle\langle\alpha|,$$ where $\hat 1$ is the identity operator.
This generalization of the usual probability rule $|\psi(\alpha)|^2 = \operatorname{Pr}(a = \alpha|\text{state} = \psi)$is necessary in certain contexts, for example when you want to talk about coherent states as a basis for your wavefunctions and they overlap so that $\langle \alpha | \beta \rangle \ne 0$ when $\alpha \ne \beta$. Or, as you're seeing, it's necessary if you want to use Gaussians centered at various positions as a position basis.
Once you have the above property, you can calculate any expectation value with those coordinates. Because trace is cyclically permutative, your expectation values $$\langle A \rangle_\rho = \operatorname{Tr} \left(\hat A~\rho\right)$$ become $$\langle A \rangle_\rho = \int_A \mathrm d\alpha~\mathcal N(\alpha)~\langle \alpha | ~\hat A ~ \rho~ |\alpha\rangle.$$ One choice which we can choose is "inspired" by the coherent states: let $\alpha = (q, p)$, $A = \mathbb R^2$, and $\lambda$ be an arbitrary constant; we can now define $$|q, p\rangle = \int_{-\infty}^\infty \mathrm dx~ \left(2 \pi \lambda^2\right)^{-1/4}~ \exp\left({(x - q)^2 \over 4 \lambda^2} + i ~ \frac{p ~ x}{\hbar}\right)~|x\rangle$$ where indeed we just treat $\langle x | x_0\rangle = \delta(x - x_0)$. Then it's not too hard to see that $\langle q, p | \hat x | q, p\rangle = q$ and $\langle q, p | \hat p | q, p\rangle = p.$ If I've done everything correct then $\mathcal N(q, p) = 1 / (2 \pi \hbar)$ no matter what $\lambda$ actually is, so you can set it to the error tolerance of your position measurement if you like.
After the measurement
We now know that the state $|\phi\rangle$ is the same as the state:$$|\phi\rangle = \hat 1 |\phi\rangle = \iint \frac{\mathrm dp~\mathrm dq}{h} |q, p\rangle \langle q, p | \phi\rangle = \iint \frac{\mathrm dp~\mathrm dq}{h} ~\phi(q, p)~ |q, p\rangle $$We project this onto the measurement outcome by accepting any $p$ but rejecting any $q \ne x_0$: $$\alpha |\phi'\rangle = \int \mathrm dp' ~|x_0, p'\rangle~ \iint \frac{\mathrm dp~\mathrm dq}{h} \langle x_0, p'| q, p\rangle ~ \phi(q, p).$$We choose $\alpha$ to normalize the resulting state and $\lambda$ to reflect our measurement accuracy and this resulting state $|\phi'\rangle$ is the state of the system after position $x_0$ is measured. Some of the momentum should then be "squeezed" by $\lambda$ in various ways as a result of this transfer-function $\langle x_0, p'|q, p\rangle$: but for very large $\lambda$ the $p, p'$ part of this transfer looks like a Fourier transform of a very wide Gaussian, and it is a very narrow Gaussian with respect to $p - p'$, thus behaving like $\delta(p - p').$
A: 
I've been taught that when you measure the value of an observable for some state of a system described by $|\psi\rangle$ then the state of the sistem "collapses" to the eigenvector associated with the eigenvalue measured. 
  ...
Following this logic, when you measure the position of a particle described by $|\psi\rangle$ and you get a value $x_0$ then the state should collapse to an eigenvector $|x_0\rangle$. However, as is obvious, this eigenvector is the dirac delta function in position basis, $\langle x|x_0\rangle = \delta (x-x_0)$.
  ...
This seems to me unacceptable since this state is non normalizable and furthermore simple things, such as $\langle \hat{x}\rangle$ are undefined for this state.

Exactly. The resolution of this puzzle is that the original assumption - the first quoted sentence -  is not universally valid. It is quite usable for discrete spectrum, for example for spins (although practice is more complicate than the simple projection gymnastics), but for continuous spectra it is wrong when taken literally. (some interpretations do not even consider projection postulate necessary in quantum theory - e. g. Ballentine's statistical interpretation).
Mathematically, this means there is no Born-normalizable wave function that would describe certainty in position or momentum.
Physically, one way out is to realize results of measurements of continuous variables have always non-vanishing uncertainty and use peaks with corresponding width to describe the system after the measurement. It is not clear to me how much useful this has been so far, though.
A: I remember asking myself the same question about ten years ago. How is it possible to have a definite particle position at all times?
The answer (and I'm not sure how it fits the density matrix scheme, as it's been a while) is related to the dispersion relation of the system).
Indeed $\delta(x-x_o)$ is not an eigenstate of any system (but for a dirac well if I'm not mistaken) and the wavepacket will disperse rapidly once the measurement was taken. 
In other words, time in and of itself will cause the initial superimposition to flatten.
Basically, you have to invoke the time-dependent Schrödinger equation to solve the problem and eventually a time dependent density matrix.
Does that make sense?
