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

For some context, I am trying to assess the capacity that certain two qubit gates have to create entanglement. To do this I am using the idea of "entangling power", where one takes their favorite entanglement measure and takes some distribution over all the possible input states averaging over them to find the resulting concurrence.

I want to take a uniform distribution over all my possible inputs which would mean taking a uniform distribution over the tensor product of two Bloch Spheres.

I'm struggling to understand how to do this.

I feel the approach should be somewhat similar to what one would need to do to solve this Exercise from Preskill's notes but the jump of sampling from one to two Bloch Spheres does not seem obvious or simple to me.

In this paper, a sample over one Bloch Sphere is taken also.

Any ideas of how I should approach this?

• You can consider a random pure state of 2 qubits, distributed according to Haar measure. This is the natural way of formalizing a uniform distribution over all pure states of 2 qubits. Is this what you want to do? Oct 27 '20 at 12:43
• Just learning about the Haar measure today. Guess it's time to read Chapter 7 of your book :) Thanks @JohnWatrous. Oct 27 '20 at 13:01
• I would not want to discourage you from reading that chapter, but you may not need all of that stuff if you just want to perform some calculations, say. Feel free to clarify if I have not understood your question correctly. In particular, it is not clear to me if you want a random 2 qubit state or a random 2 qubit product (or unentangled) state. Oct 28 '20 at 15:03
• I want the state to be arbitrary so a random 2 qubit state. Yes, I do feel like I am probably reading up on too much stuff just to justify a simple calculation. Especially since most literature is abstract and I just want N=2 right now. But I keep running into some trouble understanding how to appropriately "wieght" the distribution within the integral. Oct 28 '20 at 15:49
• OK, great. Would numerical approximations be good enough? If so, you could just average many samples using a computer. To generate a Haar random state, just let the real and imaginary parts be iid standard normal random variables and then normalize. Alternatively, you could use QETLAB for MATLAB, which has a random state function built in. Oct 28 '20 at 16:03

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: $$|\psi\rangle= \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_0\ne 0$$): $$z_i = \frac{w_i}{w_0}$$ We obtain: $$|\psi\rangle= \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{d\text{Re}(w_k)d\text{Im}(w_k)}{2\pi}$$ 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{d\text{Re}(w_k)d\text{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}^3\frac{d\text{Re}(w_k)d\text{Im}(w_k)}{\pi} = \int |w_0|^6 (1+|z_1|^2+|z_2|^2+|z_3|^2)^{-1} \delta\left(|w_0|^2- \frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}\right) \frac{d\text{Re}(w_0d\text{Im}(w_0)}{\pi}\prod_{k=1}^3 \frac{d\text{Re}(z_k)d\text{Im}(z_k)}{\pi}$$ Performing the $$w_0$$ integration in polar coordinates: We have: $$\int |w_0|^6 \delta\left(|w_0|^2- \frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}\right)\frac{d\text{Re}(w_0)d\text{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}^3\frac{d\text{Re}(z_k)d\text{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: $$|\psi\rangle= \frac{[1 \, z]^t}{\sqrt{1+|z|^2}}$$
$$d_{\mu_{S^2}} = \frac{1}{(1+|z|^2)^2}\frac{d\text{Re}(z)d\text{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: $$|\psi\rangle= \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^2\times S^2}} = \frac{1}{(1+|z_1|^2)^2(1+|z_2|^2)^2}\prod_{k=1}^2\frac{d\text{Re}(z_k)d\text{Im}(z_k)}{\pi}$$
Practically, integrals of polynomial functions over the Haar measures $$d_{\mu_{\mathbb{C}P^3}}$$ $$d_{\mu_{S^2\times S^2}}$$ can be exactly evaluated e.g. by passing to polar coordinates.