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 m\rangle$, 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_x\pm 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 m\rangle,
$$
(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 m\rangle$, 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 m\rangle
=
\langle \ell m' kq|\ell m\rangle \:
\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 m\rangle
=
\langle \tfrac12 m' 1q|\tfrac12 m\rangle \:
\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 m\rangle
=
\langle \tfrac12 m' 1q|\tfrac12 m\rangle \:
\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 m\rangle
=
\frac{
\langle \tfrac12 ||v ||\tfrac12 \rangle
}{
\langle \tfrac12 ||S||\tfrac12 \rangle
}
\langle \tfrac12 m'|S_q|\tfrac12 m\rangle
.
$$
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 m\rangle
=
\frac{
\langle \tfrac12 ||v ||\tfrac12 \rangle
}{
\langle \tfrac12 ||S||\tfrac12 \rangle
}
\langle \tfrac12 m'|\mathbf S|\tfrac12 m\rangle
,
$$
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 $|\psi\rangle$, there is always a direction $\hat{\mathbf n} \propto \langle\psi| \mathbf S |\psi\rangle $ such that the state $|\psi\rangle$ is an eigenstate of the spin component $S_{\hat{\mathbf n}}={\mathbf S}\cdot\hat{\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=1\rangle+b|m=-1\rangle$ 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 $|\psi\rangle$ 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
$$
\langle\psi| \mathbf v |\psi\rangle
\propto
\langle\psi| \mathbf S |\psi\rangle
$$
for all states $|\psi\rangle$, and indeed we can go further and conclude that the proportionality constant $K$ in that relationship must be independent of $|\psi\rangle$, 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
$$
\langle\psi| \mathbf v |\psi\rangle
=
\frac{
\langle \tfrac12 || v ||\tfrac12 \rangle
}{
\langle \tfrac12 || S ||\tfrac12 \rangle
}
\langle\psi| \mathbf S |\psi\rangle
$$
for all states $|\psi\rangle$.
Now, that's not quite enough to conclude that $\mathbf v\propto \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 $|\psi\rangle$ with some other arbitrary state $|\phi\rangle$ as well as with the various superpositions $|\psi\rangle \pm |\phi\rangle$ and $|\psi\rangle \pm i|\phi\rangle$, and you will get enough equations to conclude that
$$
\langle\phi| \mathbf v |\psi\rangle
=
\frac{
\langle \tfrac12 || v ||\tfrac12 \rangle
}{
\langle \tfrac12 || S ||\tfrac12 \rangle
}
\langle\phi| \mathbf S |\psi\rangle
$$
for all states $|\psi\rangle$ and $|\phi\rangle$, 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.
There are two possible sources of strong CP-violation, one from the quark sector and the other in the gluon sector. The QCD Lagrangian is
$$
\mathcal L_\mathrm{QCD}=-\frac14 F_{\mu\nu}^a F^{a\mu\nu}+\frac{g^2\Theta_\mathrm{YM}}{16\pi^2 N_c}F_{\mu\nu}^a\tilde F^{a\mu\nu}+\bar qi\not Dq-\bar qM e^{i\Theta_Q\gamma^5/N_f}q
$$
with terms corresponding to gluon field strength, Yang-Mills vacuum parameter, quark kinetic term and quark mass term (coming from the Yukawa terms after Higgs SSB and a field redefinition) respectively. $\Theta_\mathrm{YM}$ and $\Theta_Q$ are in general complex numbers, and their respective terms individually violate CP-symmetry for non-zero values.
However, in the quantum path integral, we can perform a field redefinition on the quarks, eliminating the theta term in the quark mass matrix. This brings with it a chiral anomaly due to the non-invariance of the path integral measure under such a shift, and a change in $\Theta_Q$ induces a corresponding change in $\Theta_\mathrm{YM}$. Thus only their difference is a meaningful contributor to strong-CP violation, and this parameter is denoted by $\bar\Theta\equiv\Theta_\mathrm{YM}-\Theta_\mathrm{Q}$. For unknown reasons, this value is very finely-tuned to zero, causing little to no CP-violation in the strong sector - this is the strong-CP problem. Since $\bar\Theta$ appears in the factor $\exp(in\bar\Theta)$ in the path integral, it is periodic with period $2\pi$
You can also perform this analysis on U(1) theory to show that there is no CP-violation in quantum electrodynamics even without fine-tuning: there are no instantons to prevent the Yang-Mills theta term from decoupling.
tl;dr - there appear to be two CP-violating terms in the strong Lagrangian, only the difference $\bar\Theta$ of their respective parameters has a measurable effect.
Best Answer
The implication would be that there's yet unknown (and not accounted for in the Standard model) source of the $\mathcal{P}$ and $\mathcal{CP}$-violation likely caused by new particles with masses near $\sim 100\,\mathrm{TeV}$.
$\mathcal{P}$ violation means that the natural laws differ when you reflect spatial coordinates. $\mathcal{C}$ violation means that the natural laws are not symmetric under exchange of all particles to their corresponding antiparticles. $\mathcal{CP}$ violation means that if you combine the spatial reflection with exchange of the particles to the antiparticles, this is still not a symmetry. Because of the $\mathcal{CPT}$ theorem, it is equivalent to the natural laws changing with respect to the time reversal.
All these symmetries are violated already in the Standard model but only through a well-known mechanism. There is also no way to add the extra sources of the $\mathcal{CP}$ violation with the particles present in the Standard model. So the discovery of the much larger electron EDM would mean that there is a new physics beyond the Standard model and already within a reach of the next generation of the particle accelerators.
This may partly explain why there are much more baryons (protons and neutrons) than antibaryons in the universe. Partly, because some unknown phase transition in the early universe would also be required to provide the conditions for the disbalance between the baryon-producing and baryon-destroying processes.
While in case of such discovery we will know that there is some new physics beyond the Standard model, we sadly will not know what exactly happens at those energy scales (maybe supersymmetry, maybe something else!) For this we will need new collider, much more powerful than LHC. Nevertheless this will be a very precious piece of information and will likely stimulate the construction of this new collider.
P.S. Don't treat too seriously all this popular talk about shape of the electron.