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?