A magnetic dipole transition can be modelled as a time-dependent perturbation $V_{\text{md}}(t) = {e\over 2 m}(\vec{L} + 2\vec{S})\cdot \vec{B}e^{-i \omega t}$. Fermi's Golden Rule tells us that the transition rate for $b-X,1$ is proportional to the matrix element of the perturbation between the initial and final states,
$$W \propto \langle \psi_b|{e\over 2 m}(\vec{L} + 2\vec{S})\cdot \vec{B}|\psi_{X}\rangle,$$
where $|\psi_b\rangle$ is the excited state and $|\psi_{X}\rangle$ is the ground state (with three possible $M_S$ values.)
The effect of $\vec L$ and $\vec S$ will be to turn the final state into some combination of the triplet states, but it won't change $J$. Therefore we might expect the transition to be 'spin-forbidden':
$$W \propto \langle b^1\Sigma_g^+ |X^3\Sigma_{g,M_S=0,\pm1}^-\rangle = \langle J=0 | J=1 \rangle = 0.$$
This is where the spin-orbit coupling comes into play. Spin-orbit coupling is the reason why the singlet $(b)$ state has a higher energy than the triplet $(X)$ states. It is a perturbation of the form $V_\text{SO}={\mu\over\hbar}\vec{L}\cdot\vec{S}$, which can be rewritten as $V_\text{SO}={\mu \over 2\hbar}(J^2-L^2-S^2)$. In a spherically symmetric system like the helium atom, this perturbation commutes with the Hamiltonian, so all you get is a shift in the energy of the triplet (L=1) and singlet (L=0) states. However, in a linear molecule like $O_2$ you lose the spherical symmetry, so $[L^2,H]\neq0$ and in addition to an energy shift, you also get some mixing of the unperturbed eigenstates, so that the excited state is not exactly $|b^1\Sigma_g^+\rangle$, but rather $|\psi_b\rangle = c_1|b^1\Sigma_g^+\rangle + c_2|X^3\Sigma_{g,M_S=0}^-\rangle$. This mixing of J=0 and J=1 states is what allows $W$ to have a nonzero value. Since we can write $S_x = S_+ + S_-$, there will be a term in the transition rate like
$$W\propto c_2^*\langle X^3\Sigma_{g,M_S=0}^-|S_{\pm}|X^3\Sigma_{g,M_S=\mp1}^-\rangle+\cdots \neq 0.$$
Does this help? I know this is a bit hand-wavy so let me know if I can clarify anything.
There are five relevant quadrupole moment operators, and when labeled by the change $\Delta m$ in angular momentum projection they read
$$
\begin{array}{c|ccccc}
\Delta m & -2 & -1 & 0 & 1 & 2\\
\hat Q_{2m}& (x-iy)^2 & (x-iy)z & x^2+y^2-2z^2 & (x+iy)z & (x+iy)^2
\end{array}
$$
These operators arise as (a basis for) all the homogeneous second-degree polynomials that are 'traceless' in the sense that they integrate to zero over the unit sphere.
As such, quadrupole transitions with $\Delta m=0$ arise from the transition quadrupole moment operator $\hat Q_{20}=x^2+y^2-2z^2$, for which it is easy to check that
$$⟨lm|\hat Q_{20}|l'm'⟩=0 \quad\text{whenever}\ m≠m',$$
since then the longitudinal integral over $\phi$ vanishes.
A bit more physically, perhaps, transitions caused by $x^2+y^2-2z^2$ represent absorption of photons that are propagating in the $z$ direction, driven not by the electric field of the wave (as is the case for dipole transitions) but by its (generally small, but nonzero) spatial change across the atom in the direction of propagation.
Similarly, if you try to do a quadrupole transition between $l=0$ and $l'=0$, then both $m$ and $m'$ are forced to be zero, so transition amplitudes are proportional to $⟨00|\hat Q_{20}|00⟩$, and when integrated explicitly over the unit sphere this gives
\begin{align}
⟨00|\hat Q_{20}|00⟩
& =
\frac{1}{4\pi}\int
x^2+y^2-2z^2
\:\mathrm d\Omega
\quad\text{(since the states have a flat wavefunction)}
\\ & =
\frac{1}{4\pi}\int_0^{2\pi}\int_0^\pi
(\sin^2(\theta)-2\cos^2(\theta))\sin(\theta)\mathrm d\theta\mathrm d\phi
\\ & =
\frac{1}{2}\int_0^\pi
(1-3\cos^2(\theta))\sin(\theta)\mathrm d\theta
\\ & =
\frac{1}{2}\int_{-1}^1
(1-3u^2)\mathrm du
=\frac12 \left[u-u^3\right]_{-1}^1
\\ & = 0.
\end{align}
Similarly, the transition $l=0\leftrightarrow l'=0$ is always forbidden from dipole radiation upwards, because it leads essentially to the inner product between the relevant transition operator $Y_{lm}$ and the flat $Y_{00}$ that arises from the wavefunctions, and those spherical harmonics are always orthogonal.
Best Answer
Selection rules in atomic transitions (including the virtual transitions which characterize Raman scattering) occur because the incident photon is changing the internal state of the atom, and must do so in ways which conserve energy, angular momentum, parity, and perhaps some other quantum numbers.
Neutron scattering is generally done with milli-eV neutrons, but the lowest excitations in nuclei typically have mega-eV energies above the ground state. So a scattered neutron typically can’t drive any excitations in the target nucleus, and therefore doesn’t care about any of the low-lying nuclear quantum numbers. To the extent that there are quantized states, they are in the virtual compound nucleus $(\text{target}+\text n)$ whose decay towards its ground state is the mechanism for neutron capture. However, most nuclei have an extremely dense set of very broad, short-lived states that overlap with their neutron separation energy, so you don’t end up with selection rules there, either. If you need some quantum numbers in a compound nucleus, there’s almost certainly an available state which has them.
If your neutrons are energetic enough to scatter from nuclear resonances, their wavelengths are probably too short to be interesting for materials science. And there are selection rules for some capture reactions, such as the $\rm^3He(n,p)^3H$ reaction whose spin dependence is important for neutron polarimetry. But that’s nucleon scattering, which is much more exciting than neutron scattering.