Paradox about the Stern-Gerlach experiment with $B=0$, $\nabla B\ne 0$ In a modern interpretation of the historical Stern-Gerlach experiment, a beam of neutral silver atoms, each with spin 1/2, was sent in the z direction through a nonuniform magnetic field having both a nonzero value and a nonzero derivative, i.e., $B\ne0$ and $\nabla B\ne 0$. In reality the beam was sent through a wide rectangular collimator, but for simplicity let's consider a simpler version of the experiment in which the beam enters the spectrometer only along the z axis. The symmetry of the magnet poles is such that $B_x=0$, $\partial B_x/\partial y=0$ and $\partial B_y/\partial x=0$. However, the following are all nonzero: $B_y$, $\partial B_y/\partial y$, and $\partial B_x/\partial x$. Textbook treatments focus on $\partial B_y/\partial y$, but because the magnetic field is divergence-free, and the field is symmetric with respect to displacements along the z axis, we have $\partial B_x/\partial x=-\partial B_y/\partial y$. Let's use the notation $B$ for the value of $B_y$ experienced by the beam, and $k=-\partial B_x/\partial x=\partial B_y/\partial y$. These two parameters $B$ and $k$ completely describe the field.
Classically, if the dipole moment of one of the silver atoms is $\mu$, then the force acting on it is given by
\begin{align}
F_x&=-k\mu_x\\
F_y&=k\mu_y.
\end{align}
Now $B\ne0$, so the field along the y axis causes the dipoles to experience a torque. If the magnetic moment is proportional to the angular momentum, then the magnetic moment will precess about the y axis with some frequency $\Omega$. If $B$ is big enough, then $\Omega$ is big enough so that $<F_x>=0$, and therefore the classical version of the spectrometer only measures $\mu_y$. It gives no information about $\mu_x$.
The paradox seems to arise when $B=0$ but $k\ne0$, so that the field experienced by the beam is zero, but the gradient of the field is nonzero. This can be accomplished, for example, by superimposing a uniform field on top of the one used in the historical experiment, so as to null the component of the field along the y axis without changing its gradient. Classically, the spectrometer is now a device that measures both $\mu_x$ and $\mu_y$. But quantum-mechanically, we can't measure $\sigma_x$ and $\sigma_y$ at the same time, because there is no such thing as a simultaneous eigenstate of both these operators. Van Huele outlines the paradox in [Van Heule 2005a] (p. 6 of the pdf, section 5.1) and [Van Heule 2005b] (p. 6 of the pdf) but doesn't seem to resolve it. He references [Scully 1987], which is paywalled, but the abstract makes it sound like that paper does present an analysis that resolves the problem.
What would actually be observed in this version of the experiment, and how is the apparent paradox resolved?
References
[Scully 1987] Marlan O. Scully, Willis E. Lamb, Jr., and Asim Barut, "On the theory of the Stern-Gerlach apparatus;" paywalled, abstract at http://adsabs.harvard.edu/abs/1987FoPh...17..575S
[Van Huele 2005a] http://www.hrpub.org/download/20040201/UJPA-18490184.pdf
[Van Huele 2005b] Jean-Francois Van Huele and Jared Stenson, "Stern-
Gerlach experiments: past, present, and future," Journal of the
Utah Academy of Science (2005), https://www.physics.byu.edu/faculty/vanhuele/redundantheory/SGPaper.pdf
 A: Let's define the operator $G \equiv -p_x\sigma_x + p_y\sigma_y$. This $G$ is a Hermitian operator, so it constitutes a perfectly valid quantum measurement. I think this is essentially equivalent to the measurement you're describing. But when I write it this way it's less mysterious!
For example, here is $G$ in the z-spin-basis:
$$G=\left(\begin{matrix}0 & p_x-ip_y \\ p_x + ip_y & 0 \end{matrix}\right)$$
What does $G$ mean? What does it do? I don't have any concise answer, I think you just have to go through some examples. If the spin starts out such-and-such spatial and spin state, and I measure $G$, what are the possible measurement results and final states? I think you can get all kinds of interesting entanglements between x-components and y-components and spin.
Anyway, there's no paradox, and I can't follow exactly why you think there is, so I suggest just working through some examples of $G$ acting on different states yourself...
A: Steve's intuition for the measurement operator is in the right direction, but it's not quite the right interaction potential. Getting that right is simple, but it requires the complex task of taking the naive answer to its ultimate consequences and then understanding what it's really saying. The hamiltonian is the same as always: a direct coupling between the spin and the magnetic field:
$$
\hat H_\mathrm{int}=\mu \mathbf B(\hat{\mathbf r})\cdot \hat{\boldsymbol\sigma}.
$$
However, we know more about the magnetic field in this situation - namely, that it varies linearly with position - and we need to bring that to bear. Thus, we know that
\begin{align}
\mathbf B (\mathbf r) 
& = \mathbf B (\mathbf 0) +\mathbf r\cdot\nabla\mathbf B(\mathbf 0) +\mathcal O(r^2)
\\ & \approx -x\frac{\partial B_x}{\partial x}(\mathbf 0)\,\hat{\mathbf e}_x+y\frac{\partial B_y}{\partial y}(\mathbf 0)\,\hat{\mathbf e}_y
\\ & =k\left(-x\hat{\mathbf e}_x+y\hat{\mathbf e}_y\right),
\end{align}
to linear order, since you've specified everything else to be zero. If we put that back into our hamiltonian, then, we get
$$
\hat H_\mathrm{int}= k\mu \left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right).
$$
It's important to note here that we have not yet specified the spin of the system, so we do not know what representation of $\mathrm{SU}(2)$ the spin matrices are in - so they could be ten-by-ten matrices for all we know. Nevertheless, if we assume a spin-1/2 representation, one way to write down that hamiltonian is in the $\hat{\sigma}_z$ basis, where it reads
$$
\hat H_\mathrm{int} = k\mu \begin{pmatrix} 0 & -\hat{x}-i\hat{y} \\ -\hat{x} +i\hat{y}  & 0 \end{pmatrix}.
$$
However, you can also choose a representation along the $\hat{\sigma}_x$ basis, say, in which case it reads
$$
\hat H_\mathrm{int} = k\mu \begin{pmatrix} -\hat{x}& \phantom{+}\hat{y} \\ \phantom{+}\hat{y}& +\hat{x} \end{pmatrix},
$$
with a complete equivalence between the two.

To get some intuition for why this is right, it's helpful to step back a bit and pretend that we only had one of the two operators acting on our system, so the interaction hamiltonian reads
$$
\hat H_\mathrm{int}= k\mu\hat{x}\hat{\sigma}_x,
$$
as we like to use in the standard case. What does this operator do? Well, assuming that the measurement is impulsive, what happens is essentially that whatever our initial state $|\psi(0)\rangle$ was, after the measurement time $\tau$ it gets transferred to
$$
|\psi(\tau)\rangle 
= e^{-i\tau\hat{H}_\mathrm{int}/\hbar} |\psi(0)\rangle
= e^{-i\tau k\mu\hat{x}\hat{\sigma}_x/\hbar} |\psi(0)\rangle.
$$
That's a lot of operators inside an exponential, but it's actually rather simple to analyze if we assume that our initial state was in a $\hat{\sigma}_x$ eigenstate and it had definite momentum $p_x$, in which case the measurement acts as
$$
|\psi(\tau)\rangle
= e^{-i\tau k\mu\hat{x}\hat{\sigma}_x/\hbar} |\psi(0)\rangle
= e^{-i\tau k\mu\hat{x}\sigma_x/\hbar} |p_x,\sigma_x\rangle
= \left|p_x + k\mu\tau \, \sigma_x,\sigma_x\right\rangle,
$$
i.e. the $\hat{\sigma}_x$ in the exponent becomes a number, and $e^{-i\tau k\mu\hat{x}\sigma_x/\hbar}$ is just a displacement operator on momentum space. Thus, over the course of the measurement, the apparatus has exerted an impulse $\frac{k\mu \tau}{\hbar}\sigma_x$ on the particle, and off it shoots in that direction.

That's the minimal example, though, so now we need to look at our full hamiltonian with the tangle of $-\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y$ in it. This is a bit more complicated, because there is now no way to initialize a spin state that will reduce the spatial dynamics to an effectively-spinless interaction, as in the previous case. 
So, let's have a look at what happens in the general case, where the interaction with the magnet is now implemented by the unitary
$$
\hat{U}(\tau) = \exp\mathopen{}\left(
-i \frac{k\mu \tau}{\hbar}\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)
\right)\mathclose{}
.
$$
This is a bit of an intimidating beast but at least in the spin-1/2 case it's not that terrible to handle, because the big correlated matrix inside the exponential squares to a simple matrix:
\begin{align}
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)^2
&=
\hat{x}^2\hat{\sigma}_x^2 
-\hat{y}\hat{x}\hat{\sigma}_y\hat{\sigma}_x 
-\hat{x}\hat{y}\hat{\sigma}_x\hat{\sigma}_y
+ \hat{y}^2\hat{\sigma}_y^2
\\ &=
\hat{x}^2\hat{\sigma}_x^2 
+ \hat{y}^2\hat{\sigma}_y^2
\\ &=
\hat{x}^2 + \hat{y}^2,
\end{align}
where the first step is because spin operators anticommute (valid for arbitrary spin $j$), and the second step uses $\hat{\sigma}_x^2=\hat{\sigma}_y^2 = 1$, which is specific to spin 1/2. (For higher spins there are still useful things you can say, but they're technical and more complicated.)
This identity is useful because it lets us split any arbitrary power of the exponentiated matrix into two simple cases:
\begin{align}
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)^{2n}
&=
\left( \hat{x}^2 + \hat{y}^2\right)^n,
\quad \text{and}\\
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)^{2n+1}
&=
\left( \hat{x}^2 + \hat{y}^2\right)^n
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right),
\end{align}
which then lets us split the matrix exponential into two simple terms:
\begin{align}
\hat{U}(\tau) 
&= \exp\mathopen{}\left(
-i \frac{k\mu \tau}{\hbar} 
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)
\right)\mathclose{}
\\ & = 
\sum_{n=0}^\infty \frac{1}{n!} \left(-i\frac{k\mu \tau}{\hbar}\right)^n\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)^n
\\ & = 
\sum_{n=0}^\infty \frac{(-1)^n}{(2n)!} \left(\frac{k\mu \tau}{\hbar}\right)^{2n}\left( \hat{x}^2 + \hat{y}^2\right)^n
\\ & \qquad
-i
\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)!} \left(\frac{k\mu \tau}{\hbar}\right)^{2n+1}
\left( \hat{x}^2 + \hat{y}^2\right)^n
\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)
\\ & =
\cos\left(\frac{k\mu \tau}{\hbar} \sqrt{\hat{x}^2 + \hat{y}^2}\right)
-i\sin\left(\frac{k\mu \tau}{\hbar} \sqrt{\hat{x}^2 + \hat{y}^2}\right)
\frac{-\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y }{\sqrt{\hat{x}^2 + \hat{y}^2}}
.
\end{align}
For notational convenience, hereafter I set $\gamma = k\mu\tau/\hbar$ and $r=\sqrt{x^2+y^2}$.
Is this expression an improvement? For sure, it looks like it might not be, but it's really broken things down into really rather few important pieces, and those have relatively simple action: the operator $\hat{x}^2 + \hat{y}^2$ simply acts as the laplacian in momentum space, so it doesn't shift any directions, and there are now only two bits that act directly on the spin sector, and they both act directly instead of in some entangled exponential. But I digress.
So, how does this act on an actual spin-momentum state? Well, it's going to be complicated, that's for sure, and you're going to end up with some pretty complex entanglement between the momentum and the spin, which will then need to be traced out once the particle reaches the screen. That's going to make things complicated, but let's give it a shot. To do that, we describe the state of the spatial part of the system after the interaction, which is described by the density matrix
$$
\hat\rho = \mathrm{Tr}_S\left[\hat{U}(\tau)|\psi(0)\rangle\langle \psi(0)| \hat{U}(\tau)^\dagger\right],
$$
where we're allowed to trace out the spin since we won't be addressing it again. Here I don't see much of a substitute for diving straight in and expanding things out: assuming that the initial state $|\psi(0)\rangle = |\mathbf p\rangle|s\rangle$ had a well-defined initial momentum and it was in some initial spin state $|s\rangle$, the position representation of our reduced density matrix reads
\begin{align}
\langle \mathbf r|\hat\rho|\mathbf r'\rangle
&=
\mathrm{Tr}_S\left[\langle \mathbf r|\hat{U}(\tau)|\psi(0)\rangle\langle \psi(0)| \hat{U}(\tau)^\dagger|\mathbf r'\rangle\right]
\\&=
\mathrm{Tr}_S\left[
\left(\cos\left(\gamma r\right)
-\frac{i}{r}\sin\left(\gamma r\right)
\left(-{x}\hat{\sigma}_x + {y}\hat{\sigma}_y \right)\right)
\langle \mathbf r|\mathbf p\rangle
|s\rangle\langle s|
\right. \\ & \qquad \qquad \qquad \left.
\langle\mathbf p|\mathbf r'\rangle
\left(\cos\left(\gamma r'\right)
+\frac{i}{r'}\sin\left(\gamma r'\right)
\left(-{x}'\hat{\sigma}_x + {y}'\hat{\sigma}_y \right)\right)
\right]
\\&=
\mathrm{Tr}_S\left[
\left(\cos\left(\gamma r\right)
-\frac{i}{r}\sin\left(\gamma r\right)
\left(-{x}\hat{\sigma}_x + {y}\hat{\sigma}_y \right)\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
|s\rangle\langle s|
\right. \\ & \qquad \qquad \qquad \left.
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\left(\cos\left(\gamma r'\right)
+\frac{i}{r'}\sin\left(\gamma r'\right)
\left(-{x}'\hat{\sigma}_x + {y}'\hat{\sigma}_y \right)\right)
\right]
%\\&=
%\mathrm{Tr}_S\left[
%\left(\cos\left(\gamma r\right)
%-\frac{i}{r}\sin\left(\gamma r\right)
%\left(-{x}\hat{\sigma}_x + {y}\hat{\sigma}_y \right)\right)
%e^{i\mathbf p\cdot\mathbf r/\hbar}
%|s\rangle\langle s|
%e^{-i\mathbf p\cdot\mathbf r'/\hbar}
%\left(\cos\left(\gamma r'\right)
%+\frac{i}{r'}\sin\left(\gamma r'\right)
%\left(-{x}'\hat{\sigma}_x + {y}'\hat{\sigma}_y \right)\right)
%\right]
\\&=
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\mathrm{Tr}_S\left[
|s\rangle\langle s|
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}  \right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}' \right)
\mathrm{Tr}_S\left[
\hat{\sigma}_x|s\rangle\langle s|\hat{\sigma}_x
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\mathrm{Tr}_S\left[
\hat{\sigma}_y|s\rangle\langle s|\hat{\sigma}_y
\right]
\\ & \quad 
+
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_y
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\mathrm{Tr}_S\left[
\hat{\sigma}_x|s\rangle\langle s|
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}'\right)
\mathrm{Tr}_S\left[
\hat{\sigma}_y|s\rangle\langle s|\hat{\sigma}_x
\right]
\\ & \quad 
+
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}'\right)
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_x
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\mathrm{Tr}_S\left[
\hat{\sigma}_x|s\rangle\langle s|\hat{\sigma}_y
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\mathrm{Tr}_S\left[
\hat{\sigma}_y|s\rangle\langle s|
\right]
\end{align}
when you expand everything out. This is a mess for sure, but we can start to clean at least some parts up by calculating the trivial traces, giving
\begin{align}
\langle \mathbf r|\hat\rho|\mathbf r'\rangle
&=
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}  \right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}' \right)
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\\ & \quad 
+
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_y
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_x
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}'\right)
\mathrm{Tr}_S\left[
i|s\rangle\langle s|\hat{\sigma}_z
\right]
\\ & \quad 
+
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}'\right)
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_x
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\mathrm{Tr}_S\left[
-i|s\rangle\langle s|\hat{\sigma}_z
\right]
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\mathrm{Tr}_S\left[
|s\rangle\langle s|\hat{\sigma}_y
\right]
,
\end{align}
and then you can start putting in some definite spin states. For the sake of definiteness, let's specialize to $|s\rangle$ being an eigenstate of $\hat{\sigma}_x$, in which case $\mathrm{Tr}_S\left[|s\rangle\langle s|\hat{\sigma}_x\right]=s$, and the other two remaining traces are zero. (For a more general state, with spin $+1/2$ along some axis $\mathbf n$, the traces would return the cartesian components $\mathrm{Tr}_S\left[|s\rangle\langle s|\hat{\sigma}_j\right]=n_j$ of the axis.)
Under those assumptions, then, we have
\begin{align}
\langle \mathbf r|\hat\rho|\mathbf r'\rangle
&=
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-{x}  \right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-{x}' \right)
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\\ & \quad 
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-s{x}\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\cos\left(\gamma r'\right)
\\ & \quad 
+
\cos\left(\gamma r\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-s{x}'\right)
,
\end{align}
and we can start asking just how much we can simplify this expression. The ideal goal would be to have a pure momentum state $\langle \mathbf r|\hat\rho|\mathbf r'\rangle = e^{i \tilde{\mathbf p}\cdot\mathbf r/\hbar} e^{-i \tilde{\mathbf p}\cdot\mathbf r'/\hbar} $ with some shifted momentum $\mathbf p$, but we've lost a good deal of coherence to the entanglement with the spin, so we need to see what we can get.
To explore this, we can some rearrangement on our density matrix, and we can actually get pretty far:
\begin{align}
\langle \mathbf r|\hat\rho|\mathbf r'\rangle
&=
\left(
\cos\left(\gamma r\right)
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-s{x}  \right)
\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\left(
\cos\left(\gamma r'\right)
+
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-s{x}'\right)
\right)
\\ & \qquad 
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
.
\end{align}
If we didn't have that ugly last term, and we could ignore the fact that things are in terms of $r$ instead of $x$, then we'd be golden, and we'd have a definite shift in the momentum, as in the previous case. However, that's little less than wishful thinking, and instead we need to face up to the fact that over momentum space we have nothing close to a pure state.
It might not look like it, but we've actually come pretty far here, and we're not too far from being able to say useful things about what happens at the screen. What we really want, in fact, is the momentum representation of the density matrix in this state, whose diagonal values are what ultimately gets detected on the screen, and the only thing we need to get there is to take the Fourier transform of our current version of the result. Thus, again, let's dive right back into it:
\begin{align}
\langle \mathbf q |\hat\rho|\mathbf q'\rangle
& = 
\int
\mathrm d\mathbf r\,\mathrm d\mathbf r'
\langle \mathbf q |\mathbf r\rangle \langle \mathbf r |\hat\rho|\mathbf r'\rangle \langle \mathbf r' |\mathbf q'\rangle
\\ & = 
\int
\left[
\left(
\cos\left(\gamma r\right)
+
\frac{-i}{r}\sin\left(\gamma r\right)\left(-s{x}  \right)
\right)
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\left(
\cos\left(\gamma r'\right)
+
\frac{i}{r'}\sin\left(\gamma r'\right)\left(-s{x}'\right)
\right)
\right. \\ & \qquad \left.
+
\frac{-i}{r}\sin\left(\gamma r\right){y}
e^{i\mathbf p\cdot\mathbf r/\hbar}
e^{-i\mathbf p\cdot\mathbf r'/\hbar}
\frac{i}{r'}\sin\left(\gamma r'\right){y}'
\right]
e^{-i\mathbf q\cdot\mathbf r}
e^{i\mathbf q'\cdot\mathbf r'}
\mathrm d\mathbf r\,\mathrm d\mathbf r'
\\ & = 
\int
\left(
\cos\left(\gamma r\right)
+
is\frac{x}{r}\sin\left(\gamma r\right)
\right)
e^{i(\mathbf p-\mathbf q)\cdot\mathbf r}
\mathrm d\mathbf r
\\ & \qquad 
\times
\int
\left(
\cos\left(\gamma r'\right)
-
is\frac{x}{r}\sin\left(\gamma r\right)
\right)
e^{-i(\mathbf p-\mathbf q')\cdot\mathbf r'}
\mathrm d\mathbf r'
\\ & \qquad 
+ 
\int
-i\frac{y}{r}\sin\left(\gamma r\right)
e^{i(\mathbf p-\mathbf q)\cdot\mathbf r}
\mathrm d\mathbf r
\int
i\frac{y'}{r'}\sin\left(\gamma r'\right)
e^{-i(\mathbf p-\mathbf q')\cdot\mathbf r'}
\mathrm d\mathbf r'.
\end{align}
This is actually pretty remarkable: our huge, complicated mess of a density matrix has essentially come down to the incoherent mixture of just two states, one of which is independent of $s$,
$$
\langle \mathbf q|\psi_y\rangle 
=
-i
\int
\frac{y}{r}\sin\left(\gamma r\right)
e^{i(\mathbf p-\mathbf q)\cdot\mathbf r}
\mathrm d\mathbf r
,
$$
and which can rightfully be interpreted as a source of noise brought about by the fact that our measurement hamiltonian has $\hat{y}\hat{\sigma}_y$ terms that are at right angles to our spin state $|s\rangle$, and one which does depend on $s$, 
$$
\langle \mathbf q|\psi_x\rangle 
=
\int
\left(
\cos\left(\gamma r\right)
+
is\frac{x}{r}\sin\left(\gamma r\right)
\right)
e^{i(\mathbf p-\mathbf q)\cdot\mathbf r}
\mathrm d\mathbf r
,
$$
and which if you squint enough sort of begins to look like a momentum state displaced by $s\hbar \gamma$ along $q_x$. To put things in simpler language, for the states defined above we have
$$
\hat\rho = |\psi_x\rangle\langle \psi_x| +|\psi_y\rangle\langle \psi_y|,
$$
after the interaction with the magnet, and this can be used to find any spatial observable you care to name. If what you want is the pattern observed on the screen, you can just get it directly as $|\langle \mathbf q|\psi_x\rangle|^2 + |\langle \mathbf q|\psi_y\rangle|^2$, i.e. the incoherent mixture of the results from the two pure states above.
Now, those Fourier transforms do look like a pair of bruisers, but they're perfectly doable by going to polar coordinates: the angular integral will return a Bessel function of the form $J_\nu(|\mathbf p-\mathbf q|r)$, and the radial integral can be done in terms of hypergeometric $_2F_1$ functions at arguments that evaluate to rational functions with poles at $(\mathbf p-\mathbf q)^2 = \gamma^2$. (The convergence, on the other hand, will be pretty tricky.) Alternatively, those integrals are probably susceptible to a manageable integration by residues.
More physically, what that means is that the entanglement induced by the interaction with the magnet has spread what used to be a perfectly collimated plane wave $|\mathbf p\rangle|s\rangle$ into a blob that decays algebraically away from its center, but with a rather high likelihood this blob will be displaced along $q_x$ by the expected amount of $s\hbar \gamma$.
(If I have time I'll come back and chase those tails, but frankly it's all footwork now.)

OK, so that was a lot of work and a lot of words, but there's still some loose ends to tie up. In particular, you observe correctly that the vanishing divergence of the magnetic field forces us to have a nonzero gradient in the $y$ direction which we would very much rather ignore, as it would force us into the ugly shenanigans I've been labouring over for the past several pages. Nevertheless, most of the time it seems that you can just forget about that kind of thing and your Stern-Gerlach analysis will work just fine (or so we like to pretend), so - is that true? and if so, why? and more importantly, can that be justified in a solid quantum mechanical manner without handwaving references to "precession" or whatnot?
The answer to that is that, as you note, this kind of thing only really becomes relevant if we turn off the field at the position of the beam. In the more common situation that we have a rather strong magnetic field at the beam, we need to include that into the interaction hamiltonian, and the result of this is that the unitary effected by the magnet will have the form
$$
\hat{U}(\tau) = \exp\mathopen{}\left[
-i\left( b\hat{\sigma}_x + \gamma\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right)
\right)\right]\mathclose{}
,
$$
with a term in $b\tau \hat{\sigma}_x$ that will drive much of the dynamics. This term scuppers all of our work above, because now the square of the exponent does not simplify quite as well as it did above, so all of those manipulations need to be repeated with that in mind: you're still likely to be able to do e.g.
$$
\left( b\hat{\sigma}_x + \gamma\left( -\hat{x}\hat{\sigma}_x + \hat{y}\hat{\sigma}_y \right) \right)^2
=
(b-\gamma \hat{x})^2+ \gamma^2\hat{y}^2
,
$$
and to put that into a simplified expression for $\hat{U}(t)$, but at the end of the day in the limit of large $b/\gamma$ (which is a tricky limit as that quantity is a length) you're likely to get a good approximation to the limit by just ignoring the dynamics along $y$.
That's probably enough of a tome for now, though.
A: I came across an analysis of the evolution of the wavefunction in the Stern-Gerlach experiment in Bohm, Quantum Theory, sec. 22.6, p. 593. Although he does go into the details pretty well, he sort of punts on the issue raised in this question. He has a footnote saying, "Strictly speaking, since [div B=0], there is always an inhomogeneous component of the field, which should produce deflections in the y direction. Since these are of no interest to us here, we shall not include the y component of the field hereafter." Even so, it was helpful to have his treatment to refer to. His coordinate system is one in which the beam propagates in the x direction, and the field points in the z direction. Since I extended his treatment using his coordinate system, I'll keep using that, even though this is different from the coordinate system I used in the statement of the question.
Let $r=\sqrt{y^2+z^2}$ be the transverse radius and $\phi=\tan^{-1}(z/y)$. Let the unitless factor $A=B/kr$, and choose units such that $\mu=1$. Then the interaction Hamiltonian, including the term that Bohm left out, looks like this in the z spin basis:
$$
\frac{H}{k} =
\begin{pmatrix}
A+\sin\phi& i\cos\phi \\ -i\cos\phi& -A-\sin\phi
\end{pmatrix}r.
$$
I think this is consistent with the operator that Emilio Pisanty wrote down in his answer. It can be rewritten in a simpler form using the Pauli spin matrices,
$$
\frac{H}{k} = Ar\sigma_z+\textbf{r}'\cdot\boldsymbol{\sigma},
$$
where $\textbf{r}'=\langle 0,-y,z\rangle$ is the reflection of the transverse radius vector across the z axis. This expression depends implicitly on $\phi$, because $\textbf{r}'$ depends on $\phi$.
In the real historical Stern-Gerlach experiment, I assume that $A\gg 1$, so in this case the interaction Hamiltonian is nearly proportional to $\sigma_z$. In the impulse approximation, we basically just end up measuring a deflection that tells us the value of this operator. The semiclassical explanation for this is that $\sigma_y$ precesses rapidly about the field, so its expectation value is zero.
In a version of the experiment where the field is nulled, the situation is qualitatively different. Now $A=0$, so the Hamiltonian is simply proportional to $\textbf{r}'\cdot\boldsymbol{\sigma}$. So for example if we detect a particle emerging at the positive z axis, then we've measured it to have $\sigma_z=1$, but if we see it coming out on the positive y axis, we've measured it to have $\sigma_y=-1$. In general, if we see the particle at $\phi$, we've measured the component of its spin along the direction $\pi-\phi$.
If you form an unpolarized pencil beam with a circular collimator, then you should mainly see a uniform ring on the detector. However, portions of the beam lying at small values of $r$ will inevitably have nonnegligible values of $A$, no matter how hard you've tried to null the field, so you will still see some of the classical Stern-Gerlach two-spot pattern mixed in.
If you use an unpolarized beam with a wide, very narrow rectangular collimator, as in the historical experiment, you should see two separate spots, but the direction of the dispersion of the two spin states is defined by the long axis of the collimator. This is different from the historical experiment with $A\gg1$, where the dispersion is along the magnetic field axis.
The general case, where $A$ is nonzero but not much greater than unity, seems kind of messy and not very interesting. In the typical case where $\cos\phi\ne 0$, the eigenvectors are $\langle\cos\phi,i(A+\sin\phi\pm\delta)\rangle$, with eigenvalues $\mp r\delta$, where $\delta^2=A^2+2A\sin\phi+1$. There is a nontrivial dependence on $\phi$. There is also a special case where $\cos\phi=0$, in which case the eigenvectors and eigenvalues are not given by the expressions above but rather are just the ones for $\sigma_z$.
As for the resolution of the paradox posed by Van Huele for $A=0$, I think the idea is the following. The paradox arises basically because the semiclassical explanation involving $\langle\sigma_y\rangle=0$ fails. We don't necessarily need the semiclassical explanation, but without it, we need some other way to explain how a preferred axis for dispersion arises or doesn't arise in various versions of the experiment. When the field is nulled, there is no preferred axis defined by the spectrometer along which we can say we're measuring the spin. However, there is still a preferred axis defined by one beam particle's transverse separation from another. We can get rid of this preferred axis by making the beam sufficiently well collimated, but then we recover the classical Stern-Gerlach behavior, because the field can't be perfectly nulled, and for small enough values of $r$, $A=B/kr$ becomes large.
[EDIT]
Some time after writing this answer, I came across a 2005 master's thesis by Jared Rees Stenson that reproduces a lot of the same analysis. The PDF comes up in google searches, although I didn't find any web page that would serve as a permanent link. Stenson's equations 7.30-31 do the same change of variables I did. His eq. 7.2 gives an expression for the eigenvalues of the Hamiltonian that I've verified is the same as mine up to a constant factor. His figure 7.1 agrees with my prediction that a ring is observed. I didn't see any discussion of $\pi-\phi$ as opposed to $\phi$ for the direction of polarization. This initially made me wonder whether I was wrong on that point, or whether I was right but he just doesn't highlight this very clearly in his presentation. I then came across a blog post by someone named Marty Green who fiddled with Stenson's work and came up with a diagram showing exactly the $\pi-\phi$ geometry.
[EDIT]
Since Emilio Pisanty didn't seem convinced by this answer, I decided to try to confirm this analysis by a completely different approach. I wrote a program that numerically simulates the Schrodinger equation on a y-z grid for the Hamiltonian being discussed here. The code is here on github. The Hamiltonian is the same as in Bohm, except that where Bohm ignores both the x and y coordinates, I only ignore x, and I include the gradient of the magnetic field in the y direction, which has to exist by Gauss's law. Some results are shown below for a case where $B=0$ but $k\ne0$, which is the main case of interest in this question. The distance units are such that each square is 2 units by 2 units. The other units are such that $\hbar=1$, the mass of the electron is 1, and the magnetic dipole moment of the electron is 1 (not -1). In these units, the output below has $k=\partial B_z/\partial z=300$. The wave packet is initialized as a Gaussian profile with an rms of 0.15 distance units. In this particular simulation, I used a beam polarized along the x axis. The grid for the numerical computation is 41x41, shown as 41x41 pixels below. The boxes superimposed for reference are 5 pixels square, and are offset by 0.5 pixels to the right and 0.5 pixels downward from the center, so that the lines lie on the boundary between pixels.

Here's how to interpret this diagram. Each box shows one spin component of the wavefunction, in the z basis. The left side is |z+>, the right |z->. The axes are y and z, as labeled. You have to imagine the two squares superimposed on each other. The amplitude of the wavefunction, with an arbitrary normalization, is indicated by brightness, and the phase by a cyclic colormap defined in the legend on the diagram.
It looks to me like this mostly confirms my original pen-and-paper analysis and my interpretation of that analysis. As predicted, we get a ring on the detector. As predicted, particles observed on the positive z axis are polarized as |z+>, those on the negative z axis as |z->, those on the positive y axis as |y->, and those on negative y as |y+>. It's easy to verify the ones on the z axis: there is appreciable amplitude only for one component or the other in the z basis. For the y, you have to look carefully at the phase relationship between the two graphs. To help in doing this, I've selected two points and drawn white boxes around them. The one on the right is yellow for |z+> and red for |z->, which means that, ignoring normalization, it equals approximately i|z+>+|z->. Ignoring an arbitrary phase, this is the same as |y->. A similar analysis for the point selected on the left shows that it is is i|z+>-|z->, which is |y+>.
Here are a couple more plots for comparison. The initial wavepacket is polarized along the x axis, and you can see that because the red color for both components indicates that its wavefunction is of the form |z+>+|z->:

Here is the final state of a simulation under conditions more like those of the historical Stern-Gerlach experiment, with $B=300$ and $k=300$:

This is pretty much what we expect to see. The beam splits into two wave packets, one heading up and one heading down. The rainbow stripes show the character of a complex exponential $e^{\pm iz}$. The wavepackets' rms radii have grown somewhat, which makes sense due to dispersion.
It would be fun to play with this some more, and anyone who wants to run the code and get more insight is welcome to do download the code from github and give it a try. The code is unfortunately not very optimized at all, so it takes a long time to run. You need to do very small time steps, otherwise you get unphysical results. (The code implements a rough check for whether you're doing enough time steps, and if it thinks you're not, it refuses to run.)
