$ J^2=(L+S)^2 $ versus $ (P - \frac{q}{c} A)^2 $ The total angular momentum J can be written as $$ J^2 = L^2 + S^2 + 2  L \cdot S .$$ Now, I assume this is a simplification of a more general tensor rule that $$ (M + N )^2 = M^2 + N^2 + M \cdot N + N \cdot M ,$$ where $ M = A \otimes B $, for instance, and $ N = C \otimes D .$
I'm questioning my understanding of this rule now. I thought it came from something as simple as $$( M + N ) (M + N ) = A^2 \otimes B^2 + A C \otimes B D + C A \otimes D B + C^2 \otimes D^2 = M^2 + N^2 + M \cdot N + N \cdot M .$$
This "proof" does not extend to 3 dimensions. For example, for the exercise to compute $ (\bar P - \frac{q}{c} A)^2 $ , one may write this as $ P^2 + A'^2 + P \cdot A' + A' \cdot P $ according to Sakurai, where I've allowed the $ A' $ to swallow the constants to its left.
So, what gives? How would you go about showing this identity? If I blindly multiply as I did in the 2D case, you would receive several cross terms. Even just computing $ (\bar P)^2 $ would result in a 9-term sum since there are 3 separate terms in $ \bar P = P_x \otimes I \otimes I + I \otimes P_y ... $ etc.
I'm guessing I'm either getting the "tensor multiplication" rules wrong, or I'm getting the definition of $ ^2 $ wrong for tensor multiplication.
EDIT: after clarification for a different question, I see that $ P_x \ne P_x \otimes I \otimes I$  Idk why but I associated "different directions" to "different tensor spaces". I think this was what Zero was getting at (thanks). But the two answers below still hold so that you don't get ugly cross terms like $ P_x P_y $ or something when doing the square.
 A: Let's go slowly and explicitly spell everything out. First, $\mathbf{J}$ is a vector operator, which means it's actually a set of three operators $$\mathbf{J} = (J_x, J_y, J_z).$$
Now, as I explain here, the definition of the square of a vector operator is
$$J^2 = J_x^2 + J_y^2 + J_z^2.$$
The components $J_i$ are defined to be the sum of the angular momentum and spin components. However, these act on the position and spin degrees of freedom, which are completely independent; particles with both position and spin have states which live in the tensor product of these two spaces. To make this tensor product structure explicit, we have
$$J_i = L_i \otimes I_s + I_p \otimes S_i$$
where $I_s$ is the identity operator in spin space, and $I_p$ is the identity operator in position space. Then
$$J_i^2 = (L_i \otimes I_s + I_p \otimes S_i)^2 = L_i^2 \otimes I_s + 2 L_i \otimes S_i + I_p \otimes S_i^2$$
where we used the basic tensor product identity $(A \otimes B)(C \otimes D) = AC \otimes BD$. Finally, by definition,
$$L^2 \equiv \sum_i L_i^2 \otimes I_s, \quad S^2 \equiv \sum_i I_p \otimes S_i^2, \quad L \cdot S \equiv \sum_i L_i \otimes S_i$$
which gives the final result. The reason it looks different from $\pi^2 = (p - qA/c)^2$ is because, while both $\pi$ and $\mathbf{J}$ are vector operators, the individual operators $J_i$ are themselves built from tensor product, while the $\pi_i$ aren't.
It's understandable to be confused, because the equation you asked about combines two subtleties in notation: inner products of vector operators, and operators built from tensor products. In both cases, textbooks tend to use compact notation that suppresses the implementation details, which is perfectly reasonable, but definitely rough the first time you run into it.
A: I've accepted Dr. Zhou's answer, since I think he was getting at the following, and this is what I think would fill in the dots for my remaining challenges if anyone else is curious. I'm referring to Dr. Z's answer in the linked post as well.
When someone is saying $ whatever^2 $, they really mean you have a vector of equations $\{ whatever_x , whatever_y, whatever_z\} $ that are indexed by $ x, y, z$ even though it may not actually have 3 separate Cartesian directions (for example, $ J = L \otimes I  + I \otimes S $ does not have three separate "sub-Hilbert-spaces", it has two. $ L_x, L_y, L_z $ all live in the same space).
Now, you can list the vector equations for $ \bar J $ as follows: $ J_x = L_x \otimes I + I \otimes S_x $ and so on. You'd compute $ J_x^2 $ the way I explain in the question and Dr. Z also explains in his answer, and so forth.
You could in principle (and somewhat more naturally, since it comes equipped with directions) do it for the canonical momentum.
So, for $ \pi_x^2 $ you'd get $ ((P_x - \frac{q}{c} A_x) \otimes I \otimes I )^2 =(P_x - \frac{q}{c} A_x)^2 \otimes I \otimes I  $ and then you'd expand it as necessary.
And as a last note, my rambling about Sakurai is either deeply misinterpreted on my end or a questionable use of notation from him when he writes $ A \cdot P - P \cdot A$ (I'm using the latest edition, equation 2.341).
This is because writing $ A \cdot P $ should imply $ A_x P_x \otimes A_y P_y \otimes A_z P_z$, which is not what we should get. We should get $ \{ A_x , P_x \} \otimes I \otimes I + I \otimes \{ A_y ,P_y \} ... etc $.
One can see this if you try to work out the answer for $ A = Y \hat x $. If you think Sakurai's got the equation right, then you should get something along the lines of $ \bar P^2 + Y^2 \otimes I \otimes I + \{Y , P_x \} \otimes P_y \otimes P_z  $ which actually should look something like $ P_y^2 + P_z^2 + (P_x - \frac{q}{c} Y \hat x )^2 \otimes I \otimes I $ where I've dropped some of the $ \otimes $ for convenience.
