Can 2-qubit systems be represented in the Bloch sphere representation? I am studying a system with 2 qubits, so I need, for a given state, a Bloch representation for each qubit.
I am having difficulties because I get results that do not have sense at all. For example if we are studying a singlet (and we assume separability),
$$ |\psi\rangle=\frac{1}{\sqrt{2}}(|{\uparrow}\rangle |{\downarrow}\rangle-|{\downarrow}\rangle |{\uparrow}\rangle)  $$
For calculating the Bloch representation the logical step is to calculate the expected value of the Pauli operators,
$$\vec{n}=(n_x,n_y,n_z) \quad \text{where} \quad n_i=\left\langle \sigma_i\otimes 1 \right\rangle \quad \text{or} \quad n_i=\left\langle 1 \otimes\sigma_i \right\rangle$$
But if we do this calculations with the singlet, we get that the two qubits point to $n=(0,0,0)$, but I am a bit confused, because I do not know if the process is well done, or if in some way this result makes sense.
If someone knows if this result is logical or how to do it, will be very helpful.
 A: This is an intricate problem and there are multiple directions people have taken to solve it. The main issue is that, while a single qubit can be described by three real degrees of freedom, there are more degrees of freedom for any multi-qubit system. In general, a single Bloch sphere is insufficient to describe the state of two qubits.
What can we do? For two qubits there are four basis states. There is the singlet state $\propto |{\uparrow\downarrow}\rangle-|{\downarrow\uparrow}\rangle$ and the three triplet states $|{\uparrow\uparrow}\rangle$, $|{\downarrow\downarrow}\rangle$, and $\propto|{\uparrow\downarrow}\rangle+|{\downarrow\uparrow}\rangle$. This means that a general superposition has four complex degrees of freedom, minus one for normalization and minus another one for the global phase, leaving six remaining real degrees of freedom.  This also means that a general mixed state has 15 real degrees of freedom! We definitely can't represent everything as a point inside of the Bloch sphere.

A first thought for generalizing the Bloch sphere is to use hyperspheres. Maybe we can just look at spheres in higher dimensions, since our states are higher dimensional?  Alas, this does not work. Pure states should be on the surface of this hypersphere, so we'd expect at least $S^6$ (where the Bloch sphere is $S^2$), but then we'd have no way of accounting for the other nine free parameters in a mixed state. In the single-qubit case, two pure-state parameters told us the orientation and the third parameter, present only in mixed states, told us how close/far we were from the surface. If we have a hypersphere $S^6$, we get six parameters for the orientation and a seventh for the radius, but that's it - there's not enough freedom to describe all of the parameters of a state.
It turns out that you can represent a pair of qubits by a non-spheriodal six-dimensional surface that is much more complicated than the Bloch sphere. This is related to representations of the Lie algebra $\mathfrak{su}(4)$ (because you can get from any pure two-qubit state to any other pure two-qubit state via an SU$(4)$ operation) but is not so easy to visualize - there have been a number of attempts but not is nearly as useful/intuitive as the Bloch sphere.

When the overall state is pure, we can do some slightly more intuitive things. Ideally, one would think that the six free parameters correspond to two sets of three parameters, so that we could represent a two-qubit pure state using two Bloch spheres. We could look at the reduced density matrix for each qubit, calculate the coordinates of each qubit on the Bloch sphere, and be done. But this is not good enough, because both reduced density matrices will have the same purity, thus have the same radius on their Bloch sphere, meaning that the Bloch vectors for the pair of reduced density matrices only encode five of the six parameters of the state. So a first candidate for representing two-qubit pure states would be a pair of Bloch vectors for the reduced density matrices plus another parameter, where that one parameter has something to do with the entanglement properties of the state but is not so easy to directly interpret.
Another method for tackling two-qubit pure states is to use a "Hopf fibration" by rewriting the four complex amplitudes as quaternions parametrizing the surface of $S^7$ and then mapping this onto $S^4$ in a rather abstract fashion. This leads to an overall parametrization of the state in terms of three spheres $S^2$; two of them correspond to the Bloch spheres of each qubit and the third parametrizes the entanglement between the two (see, e.g. here).

When the overall state is pure and only exists in the triplet subspace, we can connect a bit more directly to the Bloch sphere. Then, the Majorana representation for spin systems lets us uniquely parametrize the state as a pair of points on the surface of the Bloch sphere. This method generalizes to symmetric $N$-qubit pure states, which can be represented by $N$ points on the surface of the Bloch sphere.
The way this representation works is by writing the states in a Fock basis that counts the number of qubits in each symmetrized state:
$$|\psi\rangle= \sum_{m=0}^N \psi_m |m\rangle_{\uparrow}\otimes |N-m\rangle_{\downarrow}.$$ We can make a polynomial out of the coefficients $\psi_m$, find its $N$ complex roots, plot those $N$ roots in the complex plane, then do an inverse steriographic projection to find $N$ points on the surface of a sphere.
Now, how is this in any way intuitive? If we use the second-quantized language in which creation operators create symmetric states from the vacuum, we can write each symmetric basis state as
$$|m\rangle_{\uparrow}\otimes |N-m\rangle_{\downarrow}=\frac{\left(\hat{a}^\dagger_{\uparrow}\right)^m}{\sqrt{m!}}\otimes \frac{\left(\hat{a}^\dagger_{\downarrow}\right)^{N-m}}{\sqrt{(N-m)!}}|\mathrm{vacuum}\rangle.$$ Then, our state $|\psi\rangle$ can be uniquely written as
$$\sum_{m=0}^N \psi_m |m\rangle_{\uparrow}\otimes |N-m\rangle_{\downarrow}=\mathcal{N}\prod_{k=1}^N\left(\cos\frac{\theta_k}{2}\hat{a}^\dagger_{\uparrow}+e^{i\phi_k}\sin\frac{\theta_k}{2}\hat{a}^\dagger_{\downarrow}\right)|\mathrm{vacuum}\rangle.$$ The $N$ coordinates $(\theta_k,\phi_k)$ uniquely parametrize the state and correspond exactly to the $N$ angular coodinates of the inverse steriographic projection. These $N$ Majorana "stars" on the surface of the Bloch sphere are the generalization of a single qubit to an $N$-qubit symmetric pure state. In some sense, we can interpret this as though there are $N$ qubits, each which coordinates $(\theta_k,\phi_k)$, but they behave in some symmetric way such that you can never say which qubit has which coordinates.

In all: the Bloch sphere is unique to single qubits. Any generalization to multi-qubit states tends to be less intuitive, less mathematically convenient, or only applicable to special cases (eg pure states). The next question one could ask is how to visualize superpositions of states with different numbers of qubits...
