PDF change of variables
The goal of this module is to numerically evaluate the probability density at a
point on the Bloch sphere. The probability density is known analytically on
\(\mathbb{R}^2\) for pairs of real variables \((q_1,q_2)\) as:
\[\begin{split}\begin{align}
G(q_1,q_2)&=\frac{\epsilon}{4\pi}\exp\left\{-\frac{\epsilon(q_1^2+q_2^2+1)}
{2}\right\}(\cosh(\epsilon(q_1+q_2))+\cosh(\epsilon(q_1-q_2)))
\end{align}\end{split}\]
We also have an analytic expression for mapping points in \(\mathbb{R}^2\)
to points on the Bloch sphere:
\[\begin{split}\begin{align}
\cos\theta&=\frac{2}{\cosh(\epsilon(q_1+q_2))+\cosh(\epsilon(q_1-q_2))} \\
\phi&=\operatorname{atan2}\Big(2\sinh(\epsilon q_2),\,
\sinh(\epsilon(q_1+q_2))+\sinh(\epsilon(q_1-q_2))\Big)
\end{align}\end{split}\]
where \((\theta,\phi)\) are spherical coördinates about the y-axis:
\[\begin{split}\begin{align}
z&=\sin\theta\cos\phi \\
x&=\sin\theta\sin\phi \\
y&=\cos\theta
\end{align}\end{split}\]
Let \(\vec{\omega}:=(\cos\theta,\phi)\) and \(\vec{q}:=(q_1,q_2)\), and
define \(\Omega:\vec{q}\mapsto\vec{\omega}\). The differentials in these two
sets of variables transform as:
\[\begin{split}\begin{align}
d(\cos\theta)d\phi&=\vert\operatorname{det}(\mathrm{D}\Omega)\vert dq_1dq_2
\end{align}\end{split}\]
where \(\mathrm{D}\Omega\) is the Jacobian:
\[\begin{split}\begin{align}
\begin{bmatrix}
\frac{d(\cos\theta)}{dq_1} & \frac{d(\cos\theta)}{dq_2} \\
\frac{d\phi}{dq_1} & \frac{d\phi}{dq_2}
\end{bmatrix}
\end{align}\end{split}\]
The probability density function we want is \(\tilde{G}\) such that:
\[\begin{split}\begin{align}
\tilde{G}(\cos\theta,\phi)d(\cos\theta)d\phi&=
G(\Omega^{-1}(\cos\theta,\phi))dq_1dq_2
\end{align}\end{split}\]
Using the change of variables formula we can write:
\[\begin{split}\begin{align}
\tilde{G}(\cos\theta,\phi)&=\frac{G(\Omega^{-1}(\cos\theta,\phi))}
{\vert\operatorname{det}(\mathrm{D}\Omega)\vert}
\end{align}\end{split}\]
Jacobian
We calculate the matrix elements of the Jacobian to be:
\[\begin{split}\begin{align}
\frac{d(\cos\theta)}{dq_1}&=-2\epsilon\frac{\sinh(\epsilon q_+)+
\sinh(\epsilon q_-)}{\big(\cosh(\epsilon q_+)+\cosh(\epsilon q_-)\big)^2} \\
\frac{d(\cos\theta)}{dq_2}&=-2\epsilon\frac{\sinh(\epsilon q_+)-
\sinh(\epsilon q_-)}{\big(\cosh(\epsilon q_+)+\cosh(\epsilon q_-)\big)^2} \\
\frac{d\phi}{dq_1}&=-2\epsilon\frac{\sinh(\epsilon q_2)\big(
\cosh(\epsilon q_+)+\cosh(\epsilon q_-)\big)}{\big(\sinh(\epsilon q_+)+
\sinh(\epsilon q_-)\big)^2+4\sinh^2(\epsilon q_2)} \\
\frac{d\phi}{dq_2}&=-2\epsilon\frac{\sinh(\epsilon q_2)\big(
\cosh(\epsilon q_+)-\cosh(\epsilon q_-)\big)-\cosh(\epsilon q_2)\big(
\sinh(\epsilon q_+)+\sinh(\epsilon q_-)\big)}{\big(\sinh(\epsilon q_+)+
\sinh(\epsilon q_-)\big)^2+4\sinh^2(\epsilon q_2)}
\end{align}\end{split}\]
using the shorthand \(q_\pm:=q_1\pm q_2\).
Limiting case
For some of our plots, we want to calculate \(\tilde{G}(0,\phi)\). We can
find \(\Omega^{-1}(0,\phi)\), since:
\[\begin{split}\begin{align}
0&=\theta&\Leftrightarrow \\
0&=\big(\cosh(\epsilon(q_1+q_2))-1\big)+\big(\cosh(\epsilon(q_1-q_2))-1\big)
&\Leftrightarrow \\
q_1&=q_2=0
\end{align}\end{split}\]
We easily calculate \(G(0,0)=\epsilon e^{-\epsilon/2}/2\pi\). However,
\(\vert\operatorname{det}(D\Omega)\vert\) evaluates to the indeterminate
form \(0/0\) at \(q_1=q_2=0\). The full expression is given below:
\[\vert\operatorname{det}(D\Omega)\vert=4\epsilon^2\left\vert
\frac{2s_2(s_-c_+-s_+c_-)-(s_++s_-)^2c_2}{(c_++c_-)^2((s_++s_-)^2+4s_2^2)}
\right\vert\]
where the notation uses further shorthands \(s_a:=\sinh(\epsilon q_a)\) and
\(c_a:=\cosh(\epsilon q_a)\). We can rewrite in a different (perhaps more
suggestive) form:
\[\vert\operatorname{det}(D\Omega)\vert=4\epsilon^2\frac{c_2}{(c_++c_-)^2}
\left\vert 1-\frac{4s_2^2+2\frac{s_-c_+-s_+c_-}{c_2}s_2}
{4s_2^2+(s_++s_-)^2}\right\vert\]
From this and numerical calculations, I believe
\(\lim_{(q_1,q_2)\to(0,0)}=\epsilon^2\), but I haven’t managed to prove it
yet.