Neutral quantum particle in inhomegeneous magnetic field

Question

I'm trying to understand Stern-Gerlach experiment on a computational level. Suppose we have a neutral particle with magnetic moment (e.g. a neutron), and apply an inhomogeneous magnetic field to it (let it change linearly with coordinate). As I understand, its Hamiltonian would look like:

$$\hat H=-\frac{\hbar^2}{2m}\nabla^2+\left(\frac e{mc}\right)\hat{\vec s}\vec B$$

Now the spin operator is $$\hat s_i=\frac{\hbar}2\sigma_i,$$

where $\sigma_i$ is $i$th Pauli matrix.

So, for magnetic field $\vec B=\vec e_x B_0 x$ we'd have Schrödinger 1D (Y and Z directions can be separated due to translation symmetry) equation:

$$-\frac{\hbar^2}{2m}\frac{\partial^2\psi}{\partial x^2}+\left(\frac {\hbar e}{2mc}\right)\sigma_x B_0 x\psi=i\hbar \frac{\partial\psi}{\partial t}.$$

I now try to solve this equation numerically, taking initial wave function in the following form:

$$\psi(x,t=0)=\begin{pmatrix}\psi_0(x)\\ \psi_0(x)\end{pmatrix},$$

where $\psi_0(x)$ is a gaussian wave packet with zero average momentum.

The problems start when I select $\sigma_x$ as is usually given:

$$\sigma_x=\begin{pmatrix}0&1\\1&0\end{pmatrix}.$$

The solution appears to look like showed below. I.e. both wave function components accelerate left!

I thought, what if I choose another axis as $x$, so I tried doing the same with $\sigma_y$:

$$\sigma_y=\begin{pmatrix}0&-i\\i&0\end{pmatrix}.$$

The result in the animation below. Now it's a bit better: the wavefunction at least splits into two parts, one going left, another right. But still, both parts are composed of a mix of spin-up and spin-down states, so not really what one would expect from Stern-Gerlach experiment.

Finally, I tried the last option — using $\sigma_z$:

$$\sigma_z=\begin{pmatrix}1&0\\0&-1\end{pmatrix}.$$

The result is again showed below. Finally, I get the splitting into "independent" spin parts, i.e. one spin part goes left, another one goes right.

Now, the question: how to interpret these results? Why does choice of active axis result in such drastic differences in results? How should I have done instead to get meaningful results? Shouldn't permutation of Pauli matrices not affect results?

Answer 1

I'm going to make some cosmetic changes to your equation:

1. Your magnetic field should be $\vec B = \vec e_x B_0 x/x_0$, so that $\vec B$ and $B_0$ have the same dimensions. This field doesn't obey $\vec\nabla\cdot\vec B=0$, but we'll leave that alone for now.

2. I'll write $\mu \vec\sigma\cdot\vec B$ for the magnetic energy. Notice that the neutron's magnetic moment is pretty tiny, about 50 neV/T. So-called "ultra-cold" neutrons, with energies below about 100 neV, can be confined by a magnetic field minimum, but for cold or thermal ~meV neutrons, Stern-Gerlach steering is generally negligible. (If it weren't negligible, you could use Stern-Gerlach separation to turn one neutron beam into two polarized beams going in different directions; alas, real neutron polarizers all absorb one spin state.)

3. Your equation of motion is also separable in time, so I'll consider only the spatial part of it; we can tack on a factor of $e^{-i\omega t}$ later.

4. I'll use $\phi$ and $\chi$ to give separate names to the two components of your spinor.

In matrix notation, your time-independent Schrödinger equation is $$ \left(\begin{array}{cc} \frac{-\hbar^2}{2m}\partial_x^2 & \frac{\mu B_0 x}{x_0} \\ \frac{\mu B_0 x}{x_0} & \frac{-\hbar^2}{2m}\partial_x^2 \end{array}\right) {\phi\choose\chi} = E {\phi\choose\chi}. $$

This makes it a little more obvious what's happening. When you defined your initial ensemble to have $\phi = \chi = \psi_0(x)$, you set your system into an eigenstate of $\sigma_x$! So of course the entire ensemble moves together: it's polarized along the field direction! If your initial condition is $\phi = -\chi$, you should get a packet moving in the other direction.

When you replace $\sigma_x$ with $\sigma_{y,z}$, you're effectively changing the direction of the field. Remember that the energy is $\mu\vec\sigma\cdot\vec B$; if the only nonzero term in this product is $\sigma_y B_0 x/x_0$, you're telling your model that the field varies in strength with $x$ but points along the $y$ direction. (Since those fields do obey $\vec\nabla\cdot\vec B=0$, you should prefer them anyway.) So in your second and third figure you're putting an $x$-polarized sample into a field along $y$ or along $z$, and it separates. The separation is along the direction of the field strength gradient, rather than along the direction of the field; the UCN trappers talk about "strong field seekers" and "weak field seekers".

In the field pointing along $z$, your sample separates into your spinor basis, which gives you the nice colors.

In the field pointing along $y$, your sample still separates. But instead of separating into the $\phi,\chi$ basis that you're using for colors, the diagonal states are $\phi\pm i\chi$. This gives you interference when the two wavepackets overlap.