How to compute the root of the "bit-flip" linear map $\epsilon(\rho) = U_1\rho\,U_1^* + U_2\rho\,U_2^*$

First I'll give an answer and then I'll explain how I got that answer.

For the generators $V_1,V_2$ of the square root of the bit-flip (with probability $(1-p)\in[0,1]$) we have to distinguish two cases:

• If $p\in(\frac12,1]$ then the square root is generated by the Kraus operators $$ V_1= \sqrt{ \frac{1+\sqrt{2p-1}}{2} }\begin{pmatrix}1&0\\0&1\end{pmatrix} \qquad V_2=\sqrt{ \frac{1-\sqrt{2p-1}}{2} }\begin{pmatrix}0&1\\1&0\end{pmatrix}\,.\tag{1} $$
• If $p\in[0,\frac12]$ then the square root is generated by the Kraus operators $$ V_1= \frac{\sqrt{1+\sqrt{1-2p}}}{2}\begin{pmatrix}1&-i\\-i&1\end{pmatrix} \qquad V_2=\frac{\sqrt{1-\sqrt{1-2p}}}{2}\begin{pmatrix}1&i\\i&1\end{pmatrix}\,.\tag{2} $$

Given $p\in[0,1]$ our goal is to find a completely positive, trace-preserving map $S:\mathbb C^{2\times 2}\to\mathbb C^{2\times 2}$ and corresponding Kraus operators such that $S^2=\epsilon_p$, where \begin{align} \epsilon_p:\mathbb C^{2\times 2}&\to \mathbb C^{2\times 2}\\ \rho&\mapsto U_1\rho U_1^*+\sigma_x\rho \sigma_x^*=p\rho+(1-p)\begin{pmatrix}0&1\\1&0\end{pmatrix}\rho\begin{pmatrix}0&1\\1&0\end{pmatrix} \end{align} is the bit-flip with probability $1-p$.

Before we start let us make the following observation: if $p\in(\frac12,1]$ then $\epsilon_p$ is time-independent Markovian, meaning it is the solution of a (time-independent) master equation $ \frac{d}{dt}T(t)=L\circ T(t)$, $T(0)=\operatorname{id} $ for some $t\in[0,\infty)$; here $L:\mathbb C^{2\times 2}\to\mathbb C^{2\times 2}$ is the linear map $$A\mapsto \begin{pmatrix}0&1\\1&0\end{pmatrix}A\begin{pmatrix}0&1\\1&0\end{pmatrix}-A\,.$$ In other words $L$ is the Lindbladian generated by the Pauli-x matrix which by means of a direct computation leads to $e^{tL}\equiv\epsilon_{(1+e^{-2t})/2}$ for all $t\geq 0$. This has two immediate consequences:

1. The sqaure root of $e^{tL}$ is obviously given by $e^{tL/2}$, so for the Markovian regime the bit-flip will be its own inverse -- albeit with a different probability. More precisely for all $p\in(\frac12,1]$ we find $\epsilon_p=e^{-\ln(\sqrt{2p-1})L}$; therefore given $p\in(\frac12,1]$, one solution of $$ e^{-\ln(2q-1)L}=\epsilon_q^2=\epsilon_p=e^{-\ln(\sqrt{2p-1})L} $$ is given by $q=\frac{1+\sqrt{2p-1}}2\in(\frac12,1]$. This yields the original Kraus operators but with a different weight, i.e. equation (1).
2. The composition law $T(t)\circ T(s)=e^{tL}\circ e^{sL}=e^{(t+s)L}=T(t+s)$ of the exponential suggests that the bit-flip satisfies a more general composition rule. Indeed, for all $p,q\in[0,1]$ one finds $ \epsilon_p\circ\epsilon_q=\epsilon_{2pq-p-q+1} $ which implies $$ \epsilon_p\circ\epsilon_0=\epsilon_{1-p}=\epsilon_0\circ\epsilon_p \tag{3} $$ for all $p\in[0,1]$.

Hence the case $p\in[0,\frac12]$ could potentially be reduced to the known Markovian case by means of the pure bit-flip $\epsilon_0$ using equation (3). Usually there is no easy way to compute the square root of $AB$ symbolically (as far as I know) -- but in our case $\epsilon_0$ and $\epsilon_p$ commute by (3) so our idea will be to write $$ \sqrt{\epsilon_p}=\sqrt{\epsilon_0\circ\epsilon_{1-p}}\overset{[\epsilon_0,\epsilon_{1-p}]=0}=\sqrt{\epsilon_0}\circ\sqrt{\epsilon_{1-p}}\overset{\text{(1)}}=\sqrt{\epsilon_0}\circ\epsilon_{(1+\sqrt{1-2p})/2}\tag{4} $$

The next step is to find the root of $\epsilon_0$, that is, the root of the unitary channel $\rho\mapsto\sigma_x\rho\sigma_x$. This is particularly easy as it suffices to find a (unitary) sqaure root $\sqrt{\sigma_x}$ of $\sigma_x$ because $$ (\sqrt{\sigma_x}(\cdot)\sqrt{\sigma_x}^*)^2=\sqrt{\sigma_x}\sqrt{\sigma_x}(\cdot)\big(\sqrt{\sigma_x}\sqrt{\sigma_x}\big)^*=\sigma_x(\cdot)\sigma_x^*\,. $$ The (complex) square roots of any normal matrix $A=\sum_j a_j|\psi_j\rangle\langle\psi_j|$ are obtained via functional calculus, that is, one applies the (continuous) function in question to the eigenvalues of $A$ by means of a unitary diagonalization. We compute \begin{align} \sqrt{\begin{pmatrix}0&1\\1&0\end{pmatrix}}&=\sqrt{\frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}1&0\\0&-1\end{pmatrix}\frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix}}\\ &=\frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix} \begin{pmatrix}\sqrt{1}&0\\0&\sqrt{-1}\end{pmatrix} \frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\,. \end{align} While there are four possible solutions to this -- two for each square root -- we will apply the principal branch of the square root to obtain $$ \sqrt{\sigma_x}=\frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix} \begin{pmatrix}1&0\\0&i\end{pmatrix} \frac1{\sqrt{2}}\begin{pmatrix}1&1\\1&-1\end{pmatrix}=\frac12\begin{pmatrix}1+i&1-i\\1-i&1+i\end{pmatrix}\,. $$ The corresponding unitary channel $\sqrt{\epsilon_0}$ then acts like $$ \sqrt{\sigma_x}\begin{pmatrix}a&b\\c&d\end{pmatrix}\sqrt{\sigma_x}^*=\frac12\begin{pmatrix} (a+d)+i(b-c)&(b+c)+i(a-d)\\(b+c)-i(a-d)&(a+d)-i(b-c) \end{pmatrix}\,. $$ All that is left is to insert this into equation (4) to obtain an explicit form of $\sqrt{\epsilon_p}$ for $p\in[0,\frac12]$. This leads to \begin{align} \sqrt{\epsilon_p}:\mathbb C^{2\times 2}&\to\mathbb C^{2\times 2}\\ \begin{pmatrix}a&b\\c&d\end{pmatrix}&\mapsto \frac12\begin{pmatrix} (a+d)+i\sqrt{1-2p}(b-c)&(b+c)+i\sqrt{1-2p}(a-d)\\(b+c)-i\sqrt{1-2p}(a-d)&(a+d)-i\sqrt{1-2p}(b-c) \end{pmatrix}\,. \end{align} If anyone does not trust our arguments until here, one can of course verify $\sqrt{\epsilon_p}^2=\epsilon_p$ by hand.

Finally we have to find Kraus operators which generate $\sqrt{\epsilon_p}$ for $p\in[0,\frac12]$. For this we simply write down the Choi matrix of $\sqrt{\epsilon_p}$ $$ C(\sqrt{\epsilon_p})=\frac12\begin{pmatrix} 1&i\sqrt{1-2p}&i\sqrt{1-2p}&1\\-i\sqrt{1-2p}&1&1&-i\sqrt{1-2p}\\-i\sqrt{1-2p}&1&1&-i\sqrt{1-2p}\\1&i\sqrt{1-2p}&i\sqrt{1-2p}&1 \end{pmatrix}\,. $$ Now the Kraus operators can be computed from the Choi matrix by diagonalizing the latter and by rearranging the resulting eigenvectors (by means of vectorization). One computes $C(\sqrt{\epsilon_p})=(1+\sqrt{1-2p})|\psi_+\rangle\langle\psi_+|+(1-\sqrt{1-2p})|\psi_-\rangle\langle\psi_-\rangle$ where $\{\psi_+,\psi_-\}$ is the orthonormal system $$ |\psi_+\rangle=\frac12\begin{pmatrix}1\\-i\\-i\\1\end{pmatrix}\qquad |\psi_-\rangle=\frac12\begin{pmatrix}1\\i\\i\\1\end{pmatrix}\,. $$ The Kraus operators of $\sqrt{\epsilon_p}$ then read $$ K_+=\frac{\sqrt{1+\sqrt{1-2p}}}2\begin{pmatrix}1&-i\\-i&1\end{pmatrix}\qquad K_-=\frac{\sqrt{1-\sqrt{1-2p}}}2\begin{pmatrix}1&i\\i&1\end{pmatrix}\,. $$ (This is the standard procedure to "extract" Kraus operators from any completely positive map, but if you don't trust that $K_+,K_-$ do the job feel free to check it by hand!).