Neutral quantum particle in inhomegeneous magnetic field 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?
 A: I'm going to make some cosmetic changes to your equation:


*

*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.

*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.)

*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.

*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.
