Why must the electron's electric dipole moment (EDM) always be aligned with the spin?

As explained in the comments, this is due to the Wigner-Eckart Theorem. This is a bit of a hard one to really grok, but what it really says is that if you have a quantum-mechanical system with well-defined directional characteristics (in the sense that it is in a state with a well-defined angular momentum), and you're studying the properties of an observable that has some kind of directionality (like, say, a vector-valued operator like the electric dipole moment), then there are some tight constraints on how the orientation of the observable and the orientation of the state can interact.

That's about the best I can sum it up without immediately getting more technical. So, having done that bit of hand-waving, and since the only thing I can do is get technical, I guess I'll do that.

More specifically, to use the Wigner-Eckart theorem you need to have:

• a system in an angular momentum state $|ell mrangle$, and
• an observable that "transforms like a spherical tensor", i.e. a set of $2k+1$ observables $T_{-k}^{(k)}$, $T_{-k+1}^{(k)}$, $ldots$, $T_{k-1}^{(k)}$, $T_{k}^{(k)}$, the set of whose linear transformations is closed under spatial rotations, and which follow the same rules under rotation as the spherical harmonics $Y_{kq}$.
  • An explicit example helps: vector operators fit that bill, with $k=1$, by setting $T_{0}^{(1)}=v_z$ and $T_{pm1}^{(1)}=frac{1}{sqrt{2}}(v_xpm i v_y)$ (up to a sign).

Once you have that, then the theorem dictates that the expected value of your operator in that state, $$ langle ell m'|T_{q}^{(k)}|ell mrangle, $$ (possibly including a transition to some other orientation $m'$), will split into

• a "meaningful" part, denoted by $langle ell ||T^{(k)} ||ell rangle$, which depends on which representation the state and the observable live in, i.e. on $ell$ and $k$, but not on the specific "orientation", i.e. on the specific component $q$ of the observable that you're asking for, or the orientation $m$ of the state, and which carries all of the true dynamical information about the expectation value; and
• a factor that encodes all of the dependence on the orientation $m$ and $m'$ and on the choice of the observable's component $q$, known as a Clebsch-Gordan coefficient and denoted $langle ell m' kq|ell mrangle$, but with no knowledge of what $T$ actually is.

If you put that all together for your operator $T_{q}^{(k)}$, it reads as the equation $$ langle ell m'|T_{q}^{(k)}|ell mrangle = langle ell m' kq|ell mrangle : langle ell ||T^{(k)} ||ell rangle . $$ So, let's specialize this to some vector operator $v$, like an electric dipole moment, for a spin-1/2 particle, giving $$ langle tfrac12 m'|v_q|tfrac12 mrangle = langle tfrac12 m' 1q|tfrac12 mrangle : langle tfrac12 ||v ||tfrac12 rangle , $$ and let's compare that with how the angular momentum behaves in this context: $$ langle tfrac12 m'|S_q|tfrac12 mrangle = langle tfrac12 m' 1q|tfrac12 mrangle : langle tfrac12 ||S||tfrac12 rangle , $$ where $langle tfrac12 ||S||tfrac12 rangle$ is some numerical constant.

With this, we now have enough tools to tackle the claim as you posed it:

the electron electric dipole moment (EDM) has to be collinear with the spin.

What this really means is that, as regards orientation, our vector operator is pretty much indistinguishable from the spin, i.e. $$ langle tfrac12 m'|v_q|tfrac12 mrangle = frac{ langle tfrac12 ||v ||tfrac12 rangle }{ langle tfrac12 ||S||tfrac12 rangle } langle tfrac12 m'|S_q|tfrac12 mrangle . $$ Or, multiplying by the basis vectors $hat{mathbf e}_q$ and summing over $q$, we can recover the vector character of our equation: $$ langle tfrac12 m'|mathbf v|tfrac12 mrangle = frac{ langle tfrac12 ||v ||tfrac12 rangle }{ langle tfrac12 ||S||tfrac12 rangle } langle tfrac12 m'|mathbf S|tfrac12 mrangle , $$ which simplifies to $$ mathbf v = frac{ langle tfrac12 || v ||tfrac12 rangle }{ langle tfrac12 || S ||tfrac12 rangle } mathbf S $$ as an operator equality, since the matrix elements it considers span a basis for the space. And that puts in some context into what is meant by the claim: formally speaking, they're not "parallel" as such, but for all the measurable matrix elements that matter, and for all possible components (or linear combinations of components), the two operators yield the same result modulo a multiplicative constant. Since that is about as strong a statement of parellelism as you can make in quantum mechanics about two vector operators (which, generically, won't even commute), then we just take that as-is and keep the claim in its simplified form, which is easier to remember.

---

Having said all that, though, there is more that you can say without getting all that technical, at least in the common case of spin-$1/2$ systems, and indeed without invoking the Wigner-Eckart theorem at all. More specifically, consider the following observation observation:

For a spin-$1/2$ system in an arbitrary pure state $|psirangle$, there is always a direction $hat{mathbf n} propto langlepsi| mathbf S |psirangle $ such that the state $|psirangle$ is an eigenstate of the spin component $S_{hat{mathbf n}}={mathbf S}cdothat{mathbf n}$ along that direction, with eigenvalue $+1/2$.

This is relatively easy to show via a variety of routes, but the most important part is that it is false for any higher spin. (As an example the $m=0$ of a spin-$1$ system will never be an $m=+1$ eigenstate of any other axis orientation, and any state $a|m=1rangle+b|m=-1rangle$ with unequal nonzero weights $|a|neq |b|neq 0$ is precluded from being an eigenstate of any component of the system's spin.)

Moreover, that observation has a few direct consequences:

• The state $|psirangle$ is therefore rotationally invariant about the axis $hat{mathbf n}$.
• That means that the two components of ${mathbf S}$ orthogonal to $hat{mathbf n}$ must have vanishing expectation values, or they would break the rotational invariance.
• The same is true for any vector operator $mathbf v$.

In other words, that's enough to conclude that $$ langlepsi| mathbf v |psirangle propto langlepsi| mathbf S |psirangle $$ for all states $|psirangle$, and indeed we can go further and conclude that the proportionality constant $K$ in that relationship must be independent of $|psirangle$, because all states (in spin $1/2$) are unitarily equivalent through a rotation of the coordinate axes. Putting in some convenient notation for that proportionality constant, we get that $$ langlepsi| mathbf v |psirangle = frac{ langle tfrac12 || v ||tfrac12 rangle }{ langle tfrac12 || S ||tfrac12 rangle } langlepsi| mathbf S |psirangle $$ for all states $|psirangle$.

Now, that's not quite enough to conclude that $mathbf vpropto mathbf S$ as operators, as we concluded in the rigorous Wigner-Eckart section above, but the full operator identification is not that far off: to get it, you simply have to do as you do with the polarization identities and consider the multiple equations you get when you replace $|psirangle$ with some other arbitrary state $|phirangle$ as well as with the various superpositions $|psirangle pm |phirangle$ and $|psirangle pm i|phirangle$, and you will get enough equations to conclude that $$ langlephi| mathbf v |psirangle = frac{ langle tfrac12 || v ||tfrac12 rangle }{ langle tfrac12 || S ||tfrac12 rangle } langlephi| mathbf S |psirangle $$ for all states $|psirangle$ and $|phirangle$, and therefore that $$ mathbf v = frac{ langle tfrac12 || v ||tfrac12 rangle }{ langle tfrac12 || S ||tfrac12 rangle } mathbf S $$ as operators on that spin-$1/2$ space, completing this second version of the proof.

So: is this proof better? It is certainly as rigorous as the Wigner-Eckart one (or it can be made to be), but it doesn't really inscribe itself into a larger framework, and it suggests that the result is restricted to spin $1/2$ when the Wigner-Eckart argument is far more general. So, there's some play on both sides, and both arguments are worth understanding and exploring.