How to sample from the uniform distribution over the tensor product of two Bloch spheres?

The Hilbert space of a two-qubit system is $4$-dimensional complex vector space. An arbitrary normalized vector in this space can be written as: $$|psirangle= frac{[w_0 , w_1 , w_2 , w_3]^t}{sqrt{|w_0|^2+|w_1|^2+|w_2|^2+|w_3|^2}}$$ Since an overall phase doesn't change the state of the system, we may choose the first component $w_0$ as real, thus if we define (for $w_0ne 0$): $$z_i = frac{w_i}{w_0}$$ We obtain: $$|psirangle= frac{[1 , z_1 , z_2 , z_3]^t}{sqrt{1+|z_1|^2+|z_2|^2+|z_3|^2}}$$ The coordinate functions $z_1$, $z_2$, $z_3$ parametrize almost everywhere the complex projective space $mathbb{C}P^3$ which is the state space of two-qubit systems; but this fact will not be needed as the Haar probability measure on the space will be derived below from scratch. The Euclidean measure on the complex four-dimensional vector space $mathbb{C}^4$ is given by: $$d_{mu_{mathbb{C}^4}} = prod_{k=1}^4 frac{dtext{Re}(w_k)dtext{Im}(w_k)}{2pi}$$ Obviously, this measure is invariant under the $4$-dimensional unitary group. The normalization condition defines a seven-dimensional sphere $S^7$, a measure on $S^7$ can be constructed as: $$d_{mu_{S^7 }} = int delta(|w_0|^2+|w_1|^2+|w_2|^2+|w_3|^2-1) prod_{k=1}^4 prod_{k=1}^4 frac{dtext{Re}(w_k)dtext{Im}(w_k)}{pi}$$ ($delta$ is the Dirac delta function which restrict the measure to a unit spherical shell). Obviously, this measure is also invariant under the $4$-dimensional unitary group, thus it is a Haar-measure. Substituting the equations for $z_i$ and integrating over $w_0$, we obtain a Haar-measure over $mathbb{C}P^3$, i.e., the two-qubit state space: $$ d_{mu_{mathbb{C}P^3}} = int delta(|w_0|^2(1+|z_1|^2+|z_2|^2+|z_3|^2)-1) prod_{k=1}^3frac{dtext{Re}(w_k)dtext{Im}(w_k)}{pi} = int |w_0|^6 (1+|z_1|^2+|z_2|^2+|z_3|^2)^{-1} deltaleft(|w_0|^2- frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}right) frac{dtext{Re}(w_0dtext{Im}(w_0)}{pi}prod_{k=1}^3 frac{dtext{Re}(z_k)dtext{Im}(z_k)}{pi}$$ Performing the $w_0$ integration in polar coordinates: We have: $$int |w_0|^6 deltaleft(|w_0|^2- frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}right)frac{dtext{Re}(w_0)dtext{Im}(w_0)}{pi} = (1+|z_1|^2+|z_2|^2+|z_3|^2)^{-3}$$ Thus: $$d_{mu_{mathbb{C}P^3}} = frac{1}{(1+|z_1|^2+|z_2|^2+|z_3|^2)^4}prod_{k=1}^3frac{dtext{Re}(z_k)dtext{Im}(z_k)}{pi}$$ In summary, we have a representation of a random two-qubit state vector together with a Haar probability measure on the state space.

Of course, repeating this construction for a single qubit, we get the usual Bloch vector and the round measure on the Bloch sphere: $$|psirangle= frac{[1 , z]^t}{sqrt{1+|z|^2}}$$

$$d_{mu_{S^2}} = frac{1}{(1+|z|^2)^2}frac{dtext{Re}(z)dtext{Im}(z)}{pi}$$ The Cartesian product of two Bloch spheres is the state space of a random separable two-qubit state, the state vector is given by: $$|psirangle= frac{[1 , z_1]^t otimes [1 , z_2]^t }{sqrt{(1+|z_1|^2)( 1+|z_2|^2)}}$$ And the corresponding measure: $$d_{mu_{S^2times S^2}} = frac{1}{(1+|z_1|^2)^2(1+|z_2|^2)^2}prod_{k=1}^2frac{dtext{Re}(z_k)dtext{Im}(z_k)}{pi}$$

Practically, integrals of polynomial functions over the Haar measures $d_{mu_{mathbb{C}P^3}}$ $d_{mu_{S^2times S^2}}$ can be exactly evaluated e.g. by passing to polar coordinates.