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.
I will write
$$
\vert\alpha,\ell,s;j,m_j\rangle=\sum_{m_\ell m_s} C^{jm_j}_{\ell m_\ell; s m_s}
\vert \alpha \ell m_\ell\rangle \vert \alpha s m_s\rangle\, , \tag{1}
$$
and will assume that the tensor $T^{(k)}_q$ acts only on the orbital part, i.e. on the kets $\vert \alpha \ell m_\ell\rangle $. An expansion similar to (1) can be done for the
bra $\langle \alpha',\ell',s';j',m'_j\vert $ and the Wigner Eckart theorem will produce a
sum of the type
\begin{align}
\langle \alpha';\ell',s';j'm'_j\vert T^{(k)}_q \vert\alpha,\ell,s;j,m_j\rangle
&=\sum_{m_\ell m_s m_{\ell'}}
C^{jm_j}_{\ell m_\ell; s m_s} C^{j'm'_j}_{\ell' m'_\ell; s m_s}\delta_{ss'}\delta_{m_sm'_s}
\frac{\langle \alpha' \ell'\Vert T^{(k)}\Vert \alpha \ell\rangle}{\sqrt{2\ell'+1}}
C^{\ell' m_{\ell'}}_{kq;\ell m_\ell}\, \\
&=\frac{\langle \alpha j'\Vert T^{(k)}\Vert \alpha j\rangle}{\sqrt{2j'+1}}C^{j'm'}_{kq;jm}
\end{align}
The triple sum of Clebsch's on the right is actually proportional to the product of $ C^{\ell' m_{\ell'}}_{kq;\ell m_\ell}$ and a $6j$ symbol. Alternatively, one can multiply both sides by $C^{j''m''}_{kq;jm}$ and sum over $j'',m''$ to obtain a quadruple product of Clebsch's, proportional to a single $6j$ symbol. These summations can be found in
Varshalovich, Dmitriĭ Aleksandrovich, Anatolï Nikolaevitch Moskalev, and Valerii Kel'manovich Khersonskii. Quantum theory of angular momentum. 1988.
Either way, once these operations are done (they will involve permuting indices in the CGs ) one finally obtains
\begin{align}
\langle \alpha ' j' \ell' s'\Vert T^{(k)}\Vert \alpha j \ell s\rangle &=\delta_{ss'}
\sqrt{\frac{2j'+1}{2\ell'+1}}U(s \ell j' k; j \ell') \langle \alpha' \ell'\Vert T^{(k)}\Vert \alpha \ell\rangle
\tag{2}\\
&= \delta_{ss'}(-1)^{s+\ell+j'+k}
\sqrt{(2j '+1)(2j+1)}
\left\{\begin{array}{ccc}s&\ell&j\\ \ell &j&k \end{array}\right\}\langle \alpha' \ell'\Vert T^{(k)}\Vert \alpha \ell\rangle\, ,
\end{align}
where $\left\{\begin{array}{ccc}s&\ell&j\\ \ell &j&k \end{array}\right\}$ is a $6j$ symbol. There are several
version of Eq.(2), based on symmetries of the $6j$ symbols. The version given here is from
Rowe, David J., and John L. Wood. Fundamentals of nuclear models: foundational models. World Scientific Publishing Company, 2010.
Various authors use various symbols such as the $W$ or $U$ symbols of Wigner and Racah respectively.
Best Answer
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.