I wanted to try to understand density operators better so I thought to try to actually solve a simple problem with them. I started with the Neumann equation


From which I got the time-independent one


I defined the potential

$$V(x)=0 \text{ for}, x \geq 0, x \leq a$$ $$V(x)=\infty \text{ for}, x<0, x>a$$

I defined the density operator like this because I want to get the probabilities of seeing the particle in a space x. I am not really sure how to define $\rho$ for a continuous basis. I never saw someone using an integral for $\rho$ so I keep it as a sum. $$\rho_x=\sum_{x=0}^a p_x|x\rangle\langle x|$$

After I did some calculation I got stuck here

$$\frac{d^2\rho_x}{dx^2}+\alpha^2\rho_x=0$$ $$\alpha^2=\frac{2mE}{\hbar^2}$$

The above equation I rewriten to solve it and I got it to this form, and my question is how exactly do I continue from here? $$\lambda^2+\alpha^2=0$$

If I was to try to solve this normally I will get $$\rho_x=Ae^{i\lambda x}+Be^{-i\lambda x}$$

But I don't think that I should get the exact form as if I was using the wavefunction. How to I get a result in a matrix form? Since the density operator is a matrix.

  • 1
    $\begingroup$ You can solve for wave functions and assemble the density matrix from the wave functions. $\endgroup$
    – Roger V.
    Nov 29, 2022 at 10:08

1 Answer 1



I think it's easier to work in the eigenbasis of the Hamiltonian and do the solution there, which I will show below. If you want the position basis, you need to expand the density operator as

$$ \hat{\rho}(t) = \int_{-\infty}^{\infty} \mathrm{d}x \int_{-\infty}^{\infty} \mathrm{d} x' \rho(x,x',t) | x \rangle \langle x' |, $$

substitute it into the von Neumann equation and then you end up with

$$ i \hbar \frac{\partial}{\partial t} \rho(x,x',t) = \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right] \rho(x,x',t) - \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x'^2} + V(x') \right] \rho(x,x',t).$$

The form of this equation is suggestive of a separation of variables $\rho(x,x',t) = c(t)\psi(x)\phi(x')$. It turns out that it indeed works and you will just get back the stationary Schrödinger equation for both $\psi(x)$ and $\phi(x')$:

$$ \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right] \psi(x) = E \psi(x).$$ (Actually, you get the complex conjugate equation for $\phi(x)$ but the eigenfunctions of this problem are real.) You can now find the spectrum $E_n,~n=1,2,\ldots$ and the eigenfunctions $\psi_n(x)$.

Putting it all together, write the density matrix elements in position basis in series form

$$ \rho(x,x',t) = \sum_{m,n=1}^{\infty} c_{mn}(t) \psi_n(x) \psi_m(x') $$

If you plug this back in to the von Neumann equation you get that

$$ \frac{\mathrm{d}}{\mathrm{d}t}c_{mn}(t) = \frac{1}{i \hbar} (E_m-E_n) c_{mn}(t) $$

which you can solve:

$$ c_{mn}(t) = c_{mn}(0) e^{-i(E_m-E_n)t/\hbar} $$

where the $c_{mn}(0)$ are set by the initial density matrix.

1) Solution in the energy basis

This is pretty simple actually, just expand the density operator $\hat{\rho}(t)$ as

$$ \hat{\rho}(t) = \sum_{m,n=1}^{\infty} c_{mn}(t) | E_m \rangle\langle E_n | $$

where the eigenvalues and vectors of the Hamiltonian $\hat{H}$ satisfy $\hat{H}| E_n \rangle = E_n | E_n \rangle,~n=1,2,\ldots$ and I assumed that the spectrum of the Hamiltonian is discrete, as is the case for the potential well. Then you substitute this into the von Neumann equation of motion

$$ \frac{\mathrm{d}}{\mathrm{d}t} \hat{\rho}(t) = \frac{1}{i \hbar} [\hat{H},\hat{\rho}]$$

to get

$$ \sum_{m,n=1}^{\infty} \frac{\mathrm{d}c_{mn}(t)}{\mathrm{d}t} | E_m \rangle \langle E_n | = \frac{1}{i \hbar} \sum_{m,n=1}^{\infty} c_{mn}(t) (\hat{H}| E_m \rangle \langle E_n | - | E_m \rangle \langle E_n |\hat{H}).$$

Using $\hat{H}| E_m \rangle = E_m | E_m \rangle$ and $\langle E_n |\hat{H} = \langle E_n | E_n$ you get

$$ \sum_{m,n=1}^{\infty} \frac{\mathrm{d}c_{mn}(t)}{\mathrm{d}t} | E_m \rangle \langle E_n | = \frac{1}{i \hbar} \sum_{m,n=1}^{\infty} c_{mn}(t) (E_m - E_n) | E_m \rangle \langle E_n |.$$

Now you multiply both sides of the equation from the left by $\langle E_{m'} |$ and from the right by $| E_{n'} \rangle$ and supposing the eigenbasis is orthonormal, i.e. $\langle E_n | E_m \rangle = \delta_{nm}$, you end up with

$$\frac{\mathrm{d}c_{m'n'}(t)}{\mathrm{d}t} = \frac{1}{i \hbar} c_{m'n'}(t) (E_{m'} - E_{n'}).$$

From here on the solution is easy, and it is

$$ c_{m'n'}(t) = c_{m'n'}(0) e^{-i(E_{m'}-E_{n'})t/\hbar} $$

where the $c_{m'n'}(0)$ are set by your initial conditions.

2) Reduction to a partial differential equation in the position basis

Now you write down the density operator in position representation as

$$ \hat{\rho}(t) = \int_{-\infty}^{\infty} \mathrm{d}x \int_{-\infty}^{\infty} \mathrm{d} x' \rho(x,x',t) | x \rangle \langle x' | $$ and again substitute it into the von Neumann equation of motion to get

$$ \frac{\mathrm{d}}{\mathrm{d}t}\hat{\rho}(t) = \int_{-\infty}^{\infty} \mathrm{d}x \int_{-\infty}^{\infty} \mathrm{d} x' | x \rangle \langle x' | \frac{\partial}{\partial t}\rho(x,x',t) = \frac{1}{i \hbar} \int_{-\infty}^{\infty} \mathrm{d}x \int_{-\infty}^{\infty} \mathrm{d} x' \rho(x,x',t) (\hat{H}| x \rangle \langle x' | - | x \rangle \langle x' |\hat{H}).$$

Again you use the same trick as before, multiply both sides on the left by $\langle x'' |$ and on the right by $| x''' \rangle$. Hopefully your position basis vectors are normalised as $\langle x | x' \rangle = \delta(x-x')$ and remembering that $\int_{-\infty}^{\infty} f(x) \delta(x-x') \mathrm{d}x = f(x')$, you get

$$ \frac{\partial}{\partial t} \rho(x'',x''',t) = \frac{1}{i \hbar} \int_{-\infty}^{\infty} \mathrm{d}x \rho(x,x''',t) \langle x'' | \hat{H} | x \rangle - \frac{1}{i \hbar} \int_{-\infty}^{\infty} \mathrm{d}x' \rho(x'',x',t) \langle x' | \hat{H} | x''' \rangle.$$

Before proceeding, I relabel variables for slightly easier handling, so we have

$$ \frac{\partial}{\partial t} \rho(x,x',t) = \frac{1}{i \hbar} \int_{-\infty}^{\infty} \mathrm{d}x'' \rho(x'',x',t) \langle x | \hat{H} | x'' \rangle - \frac{1}{i \hbar} \int_{-\infty}^{\infty} \mathrm{d}x'' \rho(x,x'',t) \langle x'' | \hat{H} | x' \rangle.$$

The slightly tricky part now is to figure out what the "matrix elements" $\langle x | \hat{H} | x' \rangle$ are.

To do this, go back to the Schrödinger equation in the form

$$ i \hbar \frac{\mathrm{d}}{\mathrm{d}t} | \Psi(t) \rangle = \hat{H} | \Psi(t) \rangle $$ where the ket $| \Psi(t) \rangle$ is not a wavefunction but a vector in the Hilbert space of the problem. However, it can be expanded in the position basis as $| \Psi(t) \rangle = \int_{-\infty}^{\infty} \mathrm{d}x | x \rangle \langle x | \Psi(t) \rangle$ where the coefficients $\Psi(x,t) = \langle x | \Psi(t) \rangle$ are traditionally known as the wavefunction of the system. Playing the same game as with the von Neumann equation, you can rewrite the Schrödinger equation as

$$ i \hbar \frac{\partial}{\partial t} \Psi(x,t) = \int_{-\infty}^{\infty} \mathrm{d}x' \langle x |\hat{H}| x' \rangle \Psi(x',t).$$ Standard wave mechanics tells you however that that the right-hand side (RHS) must be

$$ \int_{-\infty}^{\infty} \mathrm{d}x' \langle x |\hat{H}| x' \rangle \Psi(x',t) = \left[-\frac{\hbar}{2m} \frac{\partial^2}{\partial x^2} + V(x)\right] \Psi(x,t)$$

which then tells you what happens when you perform integrals of the form $ \int_{-\infty}^{\infty} \mathrm{d}x' \langle x | \hat{H} | x' \rangle f(x') $

Now that we have what we need, go back to the von Neumann equation

$$ i \hbar \frac{\partial}{\partial t} \rho(x,x',t) = \int_{-\infty}^{\infty} \mathrm{d}x'' \rho(x'',x',t) \langle x | \hat{H} | x'' \rangle - \int_{-\infty}^{\infty} \mathrm{d}x'' \rho(x,x'',t) \langle x'' | \hat{H} | x' \rangle $$

where we can now evaluate the first term on the RHS:

$$ i \hbar \frac{\partial}{\partial t} \rho(x,x',t) = \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right] \rho(x,x',t) - \int_{-\infty}^{\infty} \mathrm{d}x'' \rho(x,x'',t) \langle x'' | \hat{H} | x' \rangle.$$

To evaluate the second term on the RHS, note that $\langle x'' | \hat{H} | x' \rangle = \langle x' | \hat{H}^\dagger | x'' \rangle^*$ where the $*$ denotes complex conjugation and the $\dagger$ symbol means taking the adjoint. The Hamiltonian is self-adjoint (Hermitian), so we have $\langle x'' | \hat{H} | x' \rangle = \langle x' | \hat{H} | x'' \rangle^*$. But the matrix elements are real (a consequence of the time-reversal symmetry of the problem) so nothing really changes, and you can write

$$ i \hbar \frac{\partial}{\partial t} \rho(x,x',t) = \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right] \rho(x,x',t) - \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x'^2} + V(x') \right] \rho(x,x',t).$$

  • $\begingroup$ Thank you so much for the answer. A see that a big problem that I have is that I don't understand how to write the density operator correctly. Could you please explain why you chose that expansion of the density operator in the position basis? $\endgroup$ Nov 29, 2022 at 14:24
  • $\begingroup$ Depends on the level of explanation you need. If you take e.g. real square matrices $\mathbf{A} \in \mathbb{R}^{N \times N}$, they live in an $N^2$ dimensional vector space. A basis $\mathbf{B}_{ij}$ of this space can be formed from any basis $\mathbf{b}_1,\mathbf{b}_2,\ldots,\mathbf{b}_N$ of $\mathbb{R}^N$ by $\mathbf{B}_{ij}=\mathbf{b}_i\mathbf{b}_j^T$ (the outer product of every two basis vectors). So in the end you can write down any matrix as $\mathbf{A}=\sum_{i,j=1}^{N} c_{ij}\mathbf{B}_{ij}$. $\endgroup$ Nov 29, 2022 at 15:41
  • $\begingroup$ If you now consider instead operators $\hat{A}$ acting on an infinite dimensional Hilbert space with basis $|E_k\rangle$, you can write them down as $\hat{A} = \sum_{k,k'=1}^{\infty} a_{kk'} |E_k\rangle \langle E_{k'}|$. If the basis is indexed by a continuous label such as the position basis vectors $| x \rangle$, instead of the sums you write integrals and have $\hat{A} = \int \mathrm{d}x \int \mathrm{d}x' A(x,x') | x \rangle \langle x' |$. This is really heuristic, and to be fair, I never looked up the proper maths behind it, but it does at least work. $\endgroup$ Nov 29, 2022 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.