How are LS and JJ coupling different if both involve just vector addition? In LS coupling the orbital angular momenta of particles $L_i$ couple together to form L. Similarly the spin angular momenta $S_i$ separately couple together to form S. Then S and L are coupled to get J.
In JJ coupling the $L_i$ and $S_i$ of each particle is coupled first and the resultant $J_i$s then combine to form J.
From the above, LS and JJ differ in the order of combining those vectors. If coupling is vector addition of the momenta (an associative operation) then how can it depend on the order of addition? 
 A: The addition of angular momenta is associative in the sense that the operators themselves, as operators, are the same, i.e. that the operators
$$
\mathbf J = (\mathbf L_1 + \mathbf L_2) + (\mathbf S_1+\mathbf S_2) = (\mathbf L_1+\mathbf S_1)+(\mathbf L_2+\mathbf S_2)
$$
are the same. However, when we refer to the addition of angular momenta in quantum mechanics, we mean quite a bit more than that; more specifically, we refer to a re-diagonalization process to find the eigenstates of the total angular momentum operator, and this is not an associative operation.
In the end, both the LS and JJ coupling schemes have the same final goal: to get an eigenbasis for $J^2$ and $J_z$ (and also of $L_1^2$, $L_2^2$, $S_1^2$ and $S_2^2$), but they achieve more than that:


*

*LS coupling produces a joint eigenbasis of $J^2$ and $J_z$ which is also an eigenbasis of $L^2$ and $S^2$.

*JJ coupling produces a joint eigenbasis of $J^2$ and $J_z$ which is also an eigenbasis of $J_1^2$ and $J_2^2$.


The reason this is possible is that  $J^2$ and $J_z$ commute with all of $L^2$, $S^2$, $J_1^2$ and $J_2^2$. However, they're not both simultaneously possible, because neither of $L^2$ or $S^2$ commutes with either of $J_1^2$ and $J_2^2$.
To see this explicitly, you first expand out the commutator:
\begin{align}
[L^2,J_1^2]
& =
\left[ (\mathbf L_1 + \mathbf L_2)^2 , (\mathbf L_1 + \mathbf S_1)^2 \right]
\\ & =
\left[ L_1^2 + L_2^2 + 2\mathbf L_1 \cdot\mathbf L_2, L_1^2 + S_1^2 + 2\mathbf L_1 \cdot\mathbf S_1 \right]
\\ & =
4 \left[ \mathbf L_1 \cdot\mathbf L_2, \mathbf L_1 \cdot\mathbf S_1 \right]
\\ & =
4 \sum_{k,n} \left[ L_{1,k} L_{2,k}, L_{1,n} S_{1,n} \right]
\\ & =
4 \sum_{k,n} \left[ L_{1,k}, L_{1,n} \right]L_{2,k} S_{1,n}
,
\end{align}
and here is the rub: each of those squares includes a bunch of terms linear in the components of $\mathbf L_1$ on either side of the commutator, and those stack up into a single commutator of the components of $\mathbf L_1$ with all the other components of $\mathbf L_1$, and this commutator is not zero.
(If you're really curious, the full commutator evaluates to $[L^2,J_1^2] = 4 (\mathbf L_2 \times \mathbf S_1) \cdot \mathbf L_1$.)
So, what's one to do in this situation? If you're adding up four angular momenta as
$$
\mathbf J = \mathbf L_1 + \mathbf L_2 + \mathbf S_1+\mathbf S_2,
$$
what you'd ideally want is a joint eigenbasis of


*

*the total angular momentum $J^2$ and one of its components, $J_z$,

*all of the (squares of the) individual angular momenta, $L_1^2$, $L_2^2$, $S_1^2$ and $S_2^2$, as well as

*all of the (squares of the) intermediate combinations, like $L^2$, $S^2$, $J_1^2$ and $J_2^2$, as well as the intermediate-coupling operators like $(\mathbf L+\mathbf S_1)^2$ and its ilk.


The first two goals are perfectly achievable but, because the operators in the third bullet point don't all commute with each other, we need to choose a finite subset of such operators to co-diagonalize with $J^2$ and $J_z$. This is the choice that's embedded in the dichotomy between the LS and JJ couplings.
A: The sequence in which you couple angular will simplify the evaluation of matrix elements of some Hamiltonian.  Thus, for instance, in constructing states of good total angular momentum $L$, it is natural to couple together the individual angular momentum $\ell_1,\ell_2\ldots$.  On the other hand, to compute a spin-orbit term $\vec \ell\cdot \vec s$, it is simplest to work with single particle states with good $j$.
Coupling angular momenta is not exactly just vector addition.  It involves algebraic Clebsch-Gordan coefficients and sums over intermediate states.  Thus, for instance, in coupling two $s=1/2$ spins, the states $s=1,m=0$ and $s=0,m=0$ are respectively
$$
\textstyle\frac{1}{\sqrt{2}}\vert \textstyle\frac{1}{2},\frac{1}{2}\rangle
\vert\frac{1}{2}-\frac{1}{2}\rangle \pm \frac{1}{\sqrt{2}}
\vert\frac{1}{2},-\frac{1}{2}\rangle\vert\frac{1}{2},\frac{1}{2}\rangle
$$
It is not enough to simply “add” the projections $m_1+m_2$ as the correct states involve sums over single particles states.  In more general cases, with say $\ell_1$ and $\ell_2$, the coefficients in front of the products are not equal and will usually have different signs.   
The situation becomes even more complex with three or more particles since there can be multiple values of a given total angular momentum or spin.  In the case of three spin-$1/2$ particle for instance, the possible totals spins are $3/2, 1/2$ and $1/2$ again. The rules for combining states and making sure they are orthogonal are intractable using the simple vector model and require the use of Clebsch-Gordan and Racah coefficients.
