text:
# The Real Density Matrix

## I Prologue

The density matrix plays a central role in the modern theory of quantum mechanics, and an equally important role in its applications to optics, spectroscopy, and condensed matter physics. Viewed abstractly, it is a self-adjoint operator $`\rho `$ on system’s Hilbert space, the expectation values $`0\psi |\rho |\psi 1`$ of which give the probability of observing the system in the state $`|\psi `$. As a matrix, however, it is generally represented versus the operator (or “Liouville”) basis $`|ij|`$ induced by a choice of a complete orthonormal basis $`\{|ii=0,1,\mathrm{}\}`$ in the underlying Hilbert space. These complex-valued matrices $`𝝆[i|\rho |j]_{i,j}`$ are necessarily Hermitian and positive semi-definite. Their diagonal entries are the probabilities of these mutually exclusive basis states, whereas their off-diagonal entries prescribe the amounts by which the probabilities of their coherent superpositions deviate from the corresponding classical mixtures due to interference.

Another option is to use a operator basis the elements of which have rank exceeding one, so that it is not induced by any Hilbert space basis. The most common example here is the representation of operators on a two-dimensional Hilbert space by real linear combinations of Pauli matrices $`\{𝝈_0(𝐈_{2𝖣}),𝝈_1,𝝈_2,𝝈_3\}`$. In this case the basis elements themselves are self-adjoint and so can be given a physical interpretation, e.g. as the components of the Bloch or Stokes vector Bloch (1946); Feynman et al. (1957). Although arbitrary bases of self-adjoint operators could be used, for multi-particle systems it is desirable that the overall basis be induced by identical bases on each particle’s Hilbert subspace. With the Pauli matrices this leads to the so-called *product operator* representation Ernst et al. (1987), in which density operators are represented by linear combinations of all possible tensor (*ergo* Kronecker) products of the Pauli matrices, e.g. $`𝝈_k^1𝝈_{\mathrm{}}^2(𝝈_k𝝈_0)(𝝈_0𝝈_{\mathrm{}})`$. In contrast to the Hermitian case, it has not been widely recognized that this tensor product structure is reflected by the real-valued coefficients in the expansion of any density operator in terms of product operators (also known as the *coherence vector* Mahler and Weberruss (1998)). Thus if one properly arranges these coefficients in a matrix one obtains a *real* but *nonsymmetric* analog of the Hermitian density matrix with the *same* tensor product structure. This fact holds for any number of Hilbert spaces $`_k`$ ($`k=1,2,\mathrm{}`$) of arbitrary (even infinite) dimension $`L>0`$ and self-adjoint bases $`\{B_{\mathrm{}}^k\}_{\mathrm{}=1}^L`$ for the space of bounded linear operators on each.

The purpose of this paper is show how, in the case of multi-qubit (*ergo* two-state quantum) systems, one can perform all the usual Hermitian density matrix calculations entirely with the real density matrix. While the formulae are more complicated in most cases than they are with the Hermitian density matrix, we argue that they are in many respects closer to the underlying physics than the Hermitian formulae, simply because the entries of the real density matrix correspond to (expectation values of) observables. Indeed, it is well-known that the single-qubit Pauli algebra is nothing but a complex matrix representation of the geometric (or Clifford) algebra of a three-dimensional Euclidean vector space Havel and Doran (2002a). This *real* algebra in turn has been demonstrated to be a concise but versatile formalism within which to analyze and teach much of modern physics Hestenes (2003); Doran and Lasenby (2003). It makes a certain amount of sense to use a representation wherein the “reverse-even” entities in the algebra, i.e. scalars and vectors, are real while the less-familiar “reverse-odd” entities, i.e. “pseudo-scalars” and “pseudo-vectors”, are purely imaginary (cf. Baylis (1999)). The drawback of our representation is that the matrices no longer form a representation of the underlying geometric algebra, i.e. the geometric product no longer corresponds to matrix multiplication. We will leave it to the community to decide if or when the advantages outweigh the disadvantages, and offer our results simply as the outcome of an intellectual exercise.

## II Metamorphosis

We begin by introducing a bit of notation which will considerably simplify the remainder of our presentation. First, instead of the above bra-ket notation, let us write the $`2\times 2`$ elementary matrices as

$$𝐄_{00}|\mathrm{\hspace{0.17em}0}\mathrm{\hspace{0.17em}0}|,𝐄_{10}|\mathrm{\hspace{0.17em}1}\mathrm{\hspace{0.17em}0}|,𝐄_{01}|\mathrm{\hspace{0.17em}0}\mathrm{\hspace{0.17em}1}|,𝐄_{11}|\mathrm{\hspace{0.17em}1}\mathrm{\hspace{0.17em}1}|,$$
(1)
where $`|\mathrm{\hspace{0.17em}0}𝐞_0`$, $`|\mathrm{\hspace{0.17em}1}𝐞_1`$ denote an orthonormal basis for a two-dimensional Hilbert space. Then it is easily seen that, for any nonnegative integers $`i,jM2^N1`$, the $`(M+1)\times (M+1)`$ elementary matrix $`𝐄_{ij}`$ is the Kronecker product of $`2\times 2`$ elementary matrices the indices of which are the bits $`i_n,j_n\{0,\mathrm{\hspace{0.17em}1}\}`$ in the binary expansions of $`i,j`$, respectively, i.e.

$$𝐄_{ij}\left[\delta _{ik}\delta _j\mathrm{}\right]_{k,\mathrm{}=0}^{M,M}=𝐄_{i_1j_1}\mathrm{}𝐄_{i_Nj_N},$$
(2)
where the $`\delta `$’s are Kronecker deltas. In an analogous fashion, we will denote the usual $`2\times 2`$ Pauli matrices by

$$𝐏_{00}𝝈_{\mathrm{\hspace{0.17em}0}},𝐏_{10}𝝈_{\mathrm{\hspace{0.17em}1}},𝐏_{01}𝝈_{\mathrm{\hspace{0.17em}2}},𝐏_{11}𝝈_{\mathrm{\hspace{0.17em}3}}.$$
(3)
In this notation it may readily be verified that the multiplication table among the Pauli matrices may be expressed succinctly as

$$𝐏_{ij}𝐏_k\mathrm{}=ı^{(i\mathrm{}jk)(12ij)(12k\mathrm{})}𝐏_{(i+k2ik),(j+\mathrm{}2j\mathrm{})}\left(i,j,k,\mathrm{}\{0,\mathrm{\hspace{0.17em}1}\}\right),$$
(4)
where $`ı^2=1`$. This indexing scheme may be extended to all Kronecker products of these matrices in the same way as for the $`2\times 2`$ elementary matrices, i.e.

$$𝐏_{ij}=𝐏_{i_1j_1}\mathrm{}𝐏_{i_Nj_N}(0i,jM).$$
(5)
For example, if $`M=3`$ we have $`𝐏_{01}=𝝈_0𝝈_2`$, $`𝐏_{02}=𝝈_2𝝈_0`$, and $`𝐏_{03}=𝝈_2𝝈_2`$.

The Hermitian density matrix, of course, can be expanded relative to either the elementary matrix basis or the Pauli matrix basis, e.g.

$$𝝆=\underset{i,j=0}{\overset{\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1}}{}}\rho _{ij}𝐄_{ij}=\frac{1}{2}\underset{i,j=0}{\overset{1,\mathrm{\hspace{0.17em}1}}{}}\sigma _{ij}𝐏_{ij},$$
(6)
for a single qubit with $`\rho _{ij}`$ but $`\sigma _{ij}`$. Both of these bases are orthogonal relative to the Hilbert-Schmidt inner product, which is given by $`𝐗|𝐘\mathrm{tr}(\mathrm{𝐗𝐘}^{})2^N\mathrm{𝐗𝐘}^{}`$ for any number of qubits $`N>0`$. Thus there is a unique unitary superoperator $`2^{1/2}𝒰`$ which carries the Pauli to the elementary matrix basis, where the factor of $`\sqrt{2}`$ comes from $`𝐏_{ij}^2𝐏_{ij}|𝐏_{ij}=2𝐄_{ij}^2`$. This superoperator, moreover, is nearly self-adjoint since $`𝐄_{ij}|𝐏_k\mathrm{}=𝐄_k\mathrm{}|𝐏_{ij}`$ for all $`0i,j,k,\mathrm{}1`$ with the sole exception of $`𝐄_{10}|𝐏_{01}=ı𝐄_{01}|𝐏_{10}`$ and its complex conjugate. On applying $`𝒰`$ to both sides of Eq. (6), therefore, we obtain

$$𝒰(𝝆)=\frac{1}{2}\underset{i,j=0}{\overset{\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1}}{}}\rho _{ij}\sqrt{ı^{ji}/2^{|ji|}}\left(𝐏_{ji}+(1)^{(1i)j}𝐏_{ij}\right)=\underset{i,j=0}{\overset{1,\mathrm{\hspace{0.17em}1}}{}}\sigma _{ij}𝐄_{ij}.$$
(7)
We will take this as our definition of the *real density matrix* for a single qubit. Henceforth, we shall denote this by

$$\left[\begin{array}{cc}\sigma _{00}& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]𝝈_{}𝒰(𝝆)=2\left[\begin{array}{cc}𝝆𝝈_0& 𝝆𝝈_2\\ 𝝆𝝈_1& 𝝆𝝈_3\end{array}\right].$$
(8)
Note that our choice of normalization gives $`\sigma _{00}=2𝝆=1`$, so that although $`𝒰`$ is otherwise unitary the Hilbert-Schmidt norm is scaled by a factor of $`\sqrt{2}`$; specifically $`2𝝆^2=𝝈_{}^22𝝈_{}^{}𝝈_{}=1+\sigma _{10}^2+\sigma _{01}^2+\sigma _{11}^2`$.

Let us now look at some explicit representations of the superoperator $`𝒰`$. To begin with, the mapping from the Pauli to the basis elementary is clearly

$$\begin{array}{cc}\hfill 𝐄_{00}=& \frac{1}{2}\left(𝐏_{00}+𝐏_{11}\right)\hfill \\ \hfill 𝐄_{10}=& \frac{1}{2}\left(𝐏_{10}ı𝐏_{01}\right)\hfill \end{array}\begin{array}{cc}\hfill 𝐄_{01}=& \frac{1}{2}\left(𝐏_{10}+ı𝐏_{01}\right)\hfill \\ \hfill 𝐄_{11}=& \frac{1}{2}\left(𝐏_{00}𝐏_{11}\right),\hfill \end{array}$$
(9)
and hence (since coordinates are contravariant)

$$|𝝈_{}\left[\begin{array}{c}1\\ \sigma _{10}\\ \sigma _{01}\\ \sigma _{11}\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& ı& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 1& 0& 0& 1\end{array}\right]\left[\begin{array}{c}\rho _{00}\\ \rho _{10}\\ \rho _{01}\\ \rho _{11}\end{array}\right]𝓥𝓦|𝝆𝓤|𝝆,$$
(10)
where we have factored the overall superoperator’s matrix $`𝓤`$ into the product of a diagonal matrix $`𝓥`$ and a purely real one $`𝓦`$. An operator sum representation for the superoperator $`𝒲`$ may be derived from the singular value decomposition of its *Choi matrix*, i.e.

$$\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 1& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& 1\\ 0& 1& 0& 1\\ 1& 0& 1& 0\end{array}\right]\left[\begin{array}{cccc}1& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\end{array}\right]\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& 1\\ 0& 1& 0& 1\\ 1& 0& 1& 0\end{array}\right]^{},$$
(11)
where the Choi matrix (left) is obtained simply by swapping certain pairs of the entries in $`𝓦\left[w_{ij}\right]_{i,j=0}^{\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}3}}`$ Havel (2003); specifically $`w_{20}w_{01}`$, $`w_{30}w_{11}`$, $`w_{22}w_{03}`$ and $`w_{32}w_{13}`$. Observe that the left-singular vectors associated with the nonzero singular values of $`𝓦`$ can be written as “columnized” Pauli matrices, specifically $`|𝐏_{00}`$ and $`|𝐏_{10}`$ in the notation of Eq. (10), while the corresponding right-singular vectors are $`|𝐏_{11}`$ and $`|𝐏_{10}`$. It follows that the operator sum form of $`𝒲`$ is (Havel, 2003, Proposition 3)

$$𝒲(𝝆)=𝐏_{00}𝝆𝐏_{11}+𝐏_{10}𝝆𝐏_{10}=\left[\begin{array}{cc}\rho _{00}+\rho _{11}& \rho _{10}\rho _{01}\\ \rho _{10}+\rho _{01}& \rho _{00}\rho _{11}\end{array}\right].$$
(12)
Although an operator sum form for $`𝒱`$ could be obtained by this same approach, since it is diagonal a more compact representation of its action may be obtained by packing its nonzero entries into a single $`2\times 2`$ matrix which acts via the “entrywise” or *Hadamard product*$``$Havel et al. (2001), as follows:

$$\begin{array}{cc}\hfill 𝝈_{}=𝐐\left(𝝆𝐏_{11}+𝐏_{10}𝝆𝐏_{10}\right)& \left[\begin{array}{cc}1& ı\\ 1& 1\end{array}\right]\left[\begin{array}{cc}\rho _{00}+\rho _{11}& \rho _{10}\rho _{01}\\ \rho _{10}+\rho _{01}& \rho _{00}\rho _{11}\end{array}\right]\hfill \\ \hfill & \left[\begin{array}{cc}\rho _{00}+\rho _{11}& ı(\rho _{01}\rho _{10})\\ \rho _{01}+\rho _{10}& \rho _{00}\rho _{11}\end{array}\right].\hfill \end{array}$$
(13)
Since $`𝒲`$ is self-adjoint and the overall superoperator $`𝒰`$ is unitary (up to a factor of $`\sqrt{2}`$), it is easily seen that the inverse $`𝒰^1`$ can be written as

$$𝝆=\frac{1}{2}\left(\left(\overline{𝐐}𝝈_{}\right)𝐏_{11}+𝐏_{10}\left(\overline{𝐐}𝝈_{}\right)𝐏_{10}\right),$$
(14)
where the overbar indicates the complex conjugate of all the matrix entries.

The beauty of this operator sum form for the superoperator $`𝒰`$ is that the Hadamard product obeys the mixed product formula with the Kronecker product,

$$(𝐀𝐁)(𝐂𝐃)=(𝐀𝐂)(𝐁𝐃),$$
(15)
just like the usual matrix product does. Thus if we extend $`𝒰`$ to factorable multi-qubit density matrices in the obvious way,

$$𝒰(𝝆^1\mathrm{}𝝆^N)𝒰(𝝆^1)\mathrm{}𝒰(𝝆^N),$$
(16)
and thence to arbitrary multi-qubit density matrices by linearity, we immediately obtain a general expression. Explicitly, in the case of two qubits, we get

$$\begin{array}{cc}\hfill 𝝈_{}& 𝝈_{}^1𝝈_{}^2𝒰(𝝆^1)𝒰(𝝆^2)\hfill \\ \hfill =& \left(𝐐(𝝆^1𝐏_{11}+𝐏_{10}𝝆^1𝐏_{10})\right)\left(𝐐(𝝆^2𝐏_{11}+𝐏_{10}𝝆^2𝐏_{10})\right)\hfill \\ \hfill =& \left(𝐐𝐐\right)(\left(𝝆^1𝐏_{11}+𝐏_{10}𝝆^1𝐏_{10}\right)(𝝆^2𝐏_{11}+𝐏_{10}𝝆^2𝐏_{10}))\hfill \\ \hfill =& (𝐐𝐐)(\begin{array}{cc}& (𝝆^1𝐏_{11})(𝝆^2𝐏_{11})+(𝝆^1𝐏_{11})(𝐏_{10}𝝆^2𝐏_{10})\mathrm{}\hfill \\ & +(𝐏_{10}𝝆^1𝐏_{10})(𝝆^2𝐏_{11})+(𝐏_{10}𝝆^1𝐏_{10})(𝐏_{10}𝝆^2𝐏_{10}))\hfill \end{array}\hfill \\ \hfill =& (\begin{array}{cc}& 𝐐𝐐)((𝝆^1𝝆^2)(𝐏_{11}𝐏_{11})+(𝐏_{00}𝐏_{10})(𝝆^1𝝆^2)(𝐏_{11}𝐏_{10})\mathrm{}\hfill \\ & +(𝐏_{10}𝐏_{00})(𝝆^1𝝆^2)(𝐏_{10}𝐏_{11})+(𝐏_{10}𝐏_{10})(𝝆^1𝝆^2)(𝐏_{10}𝐏_{10}))\hfill \end{array}\hfill \\ \hfill & 𝐐^2\left(𝝆𝐏_{33}+𝐏_{10}𝝆𝐏_{32}+𝐏_{20}𝝆𝐏_{31}+𝐏_{30}𝝆𝐏_{30}\right),\hfill \end{array}$$
(17)
where $`𝝆𝝆^1𝝆^2`$ and $`𝐐^N`$ ($`N>0`$) denotes the $`N`$-fold Kronecker power of $`𝐐`$. It is readily verified that the general formula for $`N`$ qubits is

$$𝒰(𝝆)=𝐐^N\underset{m=\mathrm{\hspace{0.17em}0}}{\overset{M}{}}𝐏_{m,0}𝝆𝐏_{M,(Mm)},$$
(18)
where $`𝝆`$ is any (not necessarily factorable) density matrix. Similarly, the inverse is given by

$$𝒰^1(𝝈_{})=2^N\underset{m=\mathrm{\hspace{0.17em}0}}{\overset{M}{}}𝐏_{m,0}\left(\overline{𝐐}^N𝝈_{}\right)𝐏_{M,(Mm)}.$$
(19)
Evidently $`𝐐^N𝒰^1(𝝈_{})=\mathrm{\hspace{0.17em}2}^N𝒰(\overline{𝐐}^N𝝈_{})`$.

## III Reality Check

The following properties of the real density matrix $`𝝈_{}`$ are worth noting explicitly:

* In addition to being real, it is nonsymmetric with one fixed element $`\sigma _{00}=1`$;
* It contains all the same information that the Hermitian density matrix does (since they are related by the bijection $`𝒰`$).
* It is diagonal if and only if the Hermitian density matrix is diagonal (which is why we defined $`\sigma _{11}`$ to be the coefficient of the diagonal Pauli matrix $`𝝈_3`$).
* It has the same tensor product structure as the Hermitian density matrix, since (as shown by Eq. (17)) $`𝒰`$ maps Kronecker products to Kronecker products.

In addition to these nice analytic features, the real density matrix can also be quite useful for displaying the results of *quantum state tomography*: the determination of density matrices from experimental data. In most cases to date, the real or imaginary parts of the Hermitian density matrix have been displayed using two-dimensional bar graphs (see e.g (Nielsen and Chuang, 2000, §7.7.4)). Although useful, such a plot must both omit information and exhibit redundant information. Real density matrices are definitely superior in this respect, and sometimes may also exhibit the underlying symmetry of a state more clearly. Bar graphs of the real density matrix are shown below for both the diagonal Hilbert space basis as well as the Bell basis, illustrating how easily these states may be distinguished. Further examples with experimental NMR data may be found in Havel et al. (2002).

The fact that the mapping between Pauli matrix coefficients and the entries of the Hermitian density operator preserves the tensor product structure has been noted earlier by Pittenger and Rubin (2000) (and without doubt by many other researchers as well). In our present notation, their observation was based upon the following simple relation:

$$\left[\begin{array}{cc}\sigma _{00}& \sigma _{10}\\ \sigma _{11}& ı\sigma _{01}\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]\left[\begin{array}{cc}\rho _{00}& \rho _{10}\\ \rho _{11}& \rho _{01}\end{array}\right].$$
(20)
While this is certainly a simpler relation than our operator sum, the “density matrix” on the right-hand side is not the usual Hermitian one, and the mapping between the two can be written explicitly only by using operator sums, supermatrices, or the like. Our goal here is to translate the usual operations and relations on Hermitian density matrices into the real domain, and the reordering of the entries of the real density matrix as above offers no advantage for this purpose.

As our first example of such a translation, let us show how the usual criterion for the purity of the Hermitian density matrix can be carried over to the real domain:

$`1=2^N𝝆^2=`$ $`𝝆{\displaystyle \underset{m=\mathrm{\hspace{0.17em}0}}{\overset{m=M}{}}}𝐏_{m,0}\left(\overline{𝐐}^N𝝈_{}\right)𝐏_{M,Mm}`$
$`=`$ $`\left(\overline{𝐐}^N𝝈_{}\right){\displaystyle \underset{m=0}{\overset{M}{}}}𝐏_{M,Mm}𝝆𝐏_{m,0}`$
$`=`$ $`{\displaystyle \underset{m=\mathrm{\hspace{0.17em}0}}{\overset{m=M}{}}}𝐏_{m,0}𝝆𝐏_{M,Mm}\left(\overline{𝐐}^N𝝈_{}\right)^{}`$ (21)
$`=`$ $`\left(\overline{𝐐}^N𝝈_{}\right)\left(\overline{𝐐}^N𝝈_{}\right)^{}`$
$`=`$ $`\left(𝐐^N\overline{𝐐}^N𝝈_{}\right)^{}𝝈_{}=𝝈_{}^{}𝝈_{}.`$
In going to the last line, we have used the general relation $`(𝐀𝐁)𝐂^{}=(𝐀𝐂)^{}𝐁`$ for arbitrary conformant matrices $`𝐀,𝐁,𝐂`$ Lütkepohl (1996). This derivation easily generalizes to a formula for the ensemble-average expectation values of any observable with Hermitian matrix $`𝝁`$ and corresponding real matrix $`𝝂=𝒰(𝝁)`$, showing that

$$𝝁|𝝆2^N𝝁𝝆=𝝂^{}𝝈_{}=2^N𝝂|𝝈_{}.$$
(22)
We close this section by noting that the *partial trace* operation corresponds simply to extracting a principal submatrix of the real density matrix Somaroo et al. (1998).

## IV Life in the Real World

While the expectation values of observables carry over to the real domain without significant complication, things become distinctly more challenging when it comes to integrating the equations of motion. In the case of a single qubit, it is readily verified that the commutator with an arbitrary Hamiltonian $`𝝁=𝒰^1(𝝂)`$ becomes

$$\begin{array}{cc}\hfill \left[[𝝈_{},𝝂]\right]𝒰\left([𝝆,𝝁]\right)& =ı𝐏_{01}\left(𝝂𝐄_{11}𝝈_{}𝝈_{}𝐄_{11}𝝂\right)𝐏_{01}\hfill \\ & =ı\left[\begin{array}{cc}0& \sigma _{11}\nu _{10}\sigma _{10}\nu _{11}\\ \sigma _{01}\nu _{11}\sigma _{11}\nu _{01}& \sigma _{10}\nu _{01}\sigma _{01}\nu _{10}\end{array}\right],\hfill \end{array}$$
(23)
wherein the matrix entries are the components of the usual vector cross product.

This equation of motion is most simply integrated by considering the matrix representation of the commutation superoperator defined by $`𝝂`$, which we henceforth assume without loss of generality has $`\nu _{00}=0`$. Letting $`|𝐗`$ denote the column vector of height $`(M+1)^2`$ obtained by stacking the columns of the $`(M+1)\times (M+1)`$ matrix $`𝐗`$ on top of one another in left-to-right order, and applying the well-known identity

$$|\mathrm{𝐀𝐗𝐁}=\left(𝐁^{}𝐀\right)|𝐗$$
(24)
(see e.g. Lütkepohl (1996)), we find that111The factor of $`1/2`$ here does not mean that the rotation is spinorial, but rather that the rate of rotation is $`1/\sqrt{2}`$ times the Hilbert-Schmidt norm of the Pauli matrix that generates it, while we pick up another factor of $`1/\sqrt{2}`$ on transforming to the real domain.

$$\begin{array}{cc}\hfill \dot{𝝈_{}}=|\frac{ı}{2}\left[[𝝈_{},𝝂]\right]=& \frac{1}{2}\left(𝐏_{01}𝐏_{01}\right)\left(𝐏_{00}𝝂𝐄_{11}𝝂^{}𝐄_{11}𝐏_{00}\right)|𝝈_{}\hfill \\ \hfill =& \frac{1}{2}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& \nu _{11}& \nu _{01}\\ 0& \nu _{11}& 0& \nu _{10}\\ 0& \nu _{01}& \nu _{10}& 0\end{array}\right]\left[\begin{array}{c}\sigma _{00}\\ \sigma _{10}\\ \sigma _{01}\\ \sigma _{11}\end{array}\right]\frac{1}{2}𝓡_𝝂|𝝈_{}.\hfill \end{array}$$
(25)
Since $`𝓡_𝝂^{\mathrm{\hspace{0.17em}3}}=𝝂^2𝓡_𝝂`$ ($`𝝂^22𝝂^{}𝝂𝝂|𝝂`$), this one-sided matrix differential equation is easily integrated by using the Cayley-Hamilton theorem to exponentiate the (lower-right $`3\times 3`$ block of the) coefficient matrix Najfeld and Havel (1995), obtaining

$$\begin{array}{cc}& |𝝈_{}(t)=\left(𝐏_{00}𝐏_{00}+\frac{\mathrm{sin}(𝝂t/2)}{𝝂}𝓡_𝝂+\frac{1\mathrm{cos}(𝝂t/2)}{𝝂\text{}^2}𝓡_𝝂^{\mathrm{\hspace{0.17em}2}}\right)|𝝈_{}(0)\hfill \\ \hfill =& |𝝈_{}(0)+ı\mathrm{sin}(𝝂t/2)\left[[𝝈_{}(0),\widehat{𝝂}]\right](1\mathrm{cos}(𝝂t/2))\left[[\left[[𝝈_{}(0),\widehat{𝝂}]\right],\widehat{𝝂}]\right],\hfill \end{array}$$
(26)
wherein $`\widehat{𝝂}𝝂/𝝂`$ and it is readily shown that

$$𝓡_𝝂^{\mathrm{\hspace{0.17em}2}}=|𝝂𝝂|𝝂^2\left(𝐏_{00}𝐏_{00}𝐄_{00}𝐄_{00}\right)$$
(27)
is $`𝝂^2`$ times the projection onto the plane orthogonal to the unit vector $`|\widehat{𝝂}=\widehat{𝝂}|^{}`$. The whole formula can thus be expressed more geometrically as

$$|\stackrel{ˇ}{𝝈}_{}(t)=\widehat{𝝂}|\stackrel{ˇ}{𝝈}_{}|\widehat{𝝂}+\mathrm{sin}(𝝂t/2)|\widehat{𝝂}\times |\stackrel{ˇ}{𝝈}_{}+\mathrm{cos}(𝝂t/2)|\widehat{𝝂}\times \left(|\widehat{𝝂}\times |\stackrel{ˇ}{𝝈}_{}\right),$$
(28)
where $`\stackrel{ˇ}{𝝈}_{}𝝈_{}(0)𝐄_{00}`$ and “$`\times `$” is the cross product of the (last three components of the) vectors it connects. This is a standard expression for rotation of the three-dimensional vector $`|\stackrel{ˇ}{𝝈}_{}`$ about the axis $`|\widehat{𝝂}`$ by an angle $`𝝂t/2`$.

The extension of these formulae to general multi-particle commutators is not straightforward, since the tensor product of commutation superoperators is not simply related to the tensor products of their underlying operators. Nevertheless, we can give a reasonably simple formula for the two-particle commutator with a factorable Hamiltonian, which is the most important case in practice. This is based upon a geometric algebra expression for the commutator of a tensor product of two three-dimensional vectors (or equivalently in the present context, traceless $`2\times 2`$ Hermitian matrices), which is derived in Havel and Doran (2002b):

$$2[𝐀𝐂,𝐁𝐃]=𝐀|𝐁\left(𝐏_{00}[𝐂,𝐃]\right)+\left([𝐀,𝐁]𝐏_{00}\right)𝐂|𝐃.$$
(29)
Letting $`𝐚𝒰(𝐀)`$, etc. be the corresponding real matrices, this translates to:

$$2\left[[𝐚𝐜,𝐛𝐝]\right]=𝐚|𝐛\left(𝐄_{00}[[𝐜,𝐝]]\right)+\left([[𝐚,𝐛]]𝐄_{00}\right)𝐜|𝐝$$
(30)
This formula is easily extended to the case in which $`a_{00},\mathrm{},d_{00}0`$ by multilinearity; in the following, however, we will need only the case in which $`a_{00}=c_{00}=1`$, which introduces two additional terms:

$$\begin{array}{cc}\hfill \left[[𝐄_{00}𝐜,𝐛𝐝]\right]\frac{1}{2}[𝐏_{00}𝐂,𝐁𝐃]=\frac{1}{2}𝐁[𝐂,𝐃]& \frac{1}{2}𝐛[[𝐜,𝐝]];\hfill \\ \hfill \left[[𝐚𝐄_{00},𝐛𝐝]\right]\frac{1}{2}[𝐀𝐏_{00},𝐁𝐃]=\frac{1}{2}[𝐀,𝐁]𝐃& \frac{1}{2}[[𝐚,𝐛]]𝐝.\hfill \end{array}$$
(31)
Finally, we shall need the general result Lütkepohl (1996); Havel (2003):

$$|𝐗𝐘=\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)\left(|𝐗|𝐘\right),$$
(32)
for $`2\times 2`$ matrices $`𝐗`$ & $`𝐘`$, where the two-particle *commutation matrix* $`𝓚_{22}`$ is given by

$$𝓚_{22}=\underset{i,j=0}{\overset{\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}3}}{}}𝐄_{ij}𝐄_{ji}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right].$$
(33)
We are interested in the case that $`𝐚=𝝈_{}^1`$, $`𝐛=𝝂^1`$, $`𝐜=𝝈_{}^2`$ and $`𝐝=𝝂^2`$, i.e. we have a factorizable two-particle state $`𝝈_{}^1𝝈_{}^2`$ evolving under a bi-axial interaction $`𝝂^1𝝂^2`$. To express this more compactly, we define the matrix $`𝓢_{𝝂^1}`$ via

$$𝝂^1|𝝈_{}^1|𝐄_{00}+|𝝂^1=\left[\begin{array}{cccc}0& \nu _{10}^1& \nu _{01}^1& \nu _{11}^1\\ \nu _{10}^1& 0& 0& 0\\ \nu _{01}^1& 0& 0& 0\\ \nu _{11}^1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}1\\ \sigma _{10}^1\\ \sigma _{01}^1\\ \sigma _{11}^1\end{array}\right]𝓢_{𝝂^1}|𝝈_{}^1$$
(34)
with an analogous definition in $`𝝂^2|𝝈_{}^2|𝐄_{00}+|𝝂^2𝓢_{𝝂^2}|𝝈_{}^2`$. Then Eqs. (25), (30), (31) & (32) give us

$`_t|𝝈_{}^1𝝈_{}^2=`$ $`|\frac{ı}{4}\left[[𝝈_{}^1𝝈_{}^2,𝝂^1𝝂^2]\right]`$
$`=`$ $`\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)\frac{ı}{4}\begin{array}{cc}\hfill (& 𝝈_{}^1|𝝂^1|𝐄_{00}|\left[[𝝈_{}^2,𝝂^2]\right]+\hfill \\ \hfill \mathrm{}& |\left[[𝝈_{}^1,𝝂^1]\right]|𝐄_{00}𝝈_{}^2|𝝂^2+\hfill \\ \hfill \mathrm{}& |𝝂^1|\left[[𝝈_{}^2,𝝂^2]\right]+\hfill \\ \hfill \mathrm{}& |\left[[𝝈_{}^1,𝝂^1]\right]|𝝂^2)\hfill \end{array}`$ (35)
$`=`$ $`(𝐏_{00}𝓚_{22}𝐏_{00})\frac{1}{4}(𝓢_{𝝂^1}𝓡_{𝝂^2}\begin{array}{c}\hfill +𝓡_{𝝂^1}𝓢_{𝝂^2})\left(\right|𝝈_{}^1|𝝈_{}^2)\end{array}`$
$`=`$ $`\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)\frac{1}{4}\begin{array}{cc}\hfill \left(𝓢_{𝝂^1}𝓡_{𝝂^2}+𝓡_{𝝂^1}𝓢_{𝝂^2}\right)\mathrm{}& \\ \hfill \mathrm{}\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)|𝝈_{}^1𝝈_{}^2& .\hfill \end{array}`$
Since left or right multiplication of a $`4`$D column or row vector by $`𝓡_𝐱`$ gives the cross product of the last three components of that vector with $`[x_{10},x_{01},x_{11}]`$, it may be seen that $`𝓢_{𝝂^1}`$, $`𝓡_{𝝂^1}`$ are mutually annihilating (i.e. $`𝓢_{𝝂^1}𝓡_{𝝂^1}=𝓡_{𝝂^1}𝓢_{𝝂^1}=\mathrm{𝟎}`$), and similarly for $`𝓢_{𝝂^2}`$, $`𝓡_{𝝂^2}`$. As a result, the two terms on the last line of Eq. (IV) commute and their exponential factorizes. It is moreover easily shown that $`𝓢_𝐱^3=𝐱^2𝓢_\text{}𝐱`$ and $`𝓡_𝐱^3=𝐱^2𝓡_\text{}𝐱`$, so the overall integral is

$$\begin{array}{cc}\hfill |𝝈_{}(t)=\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)\mathrm{𝐄𝐱𝐩}\left(𝓢_{𝝂^1}𝓡_{𝝂^2}t/4\right)\mathrm{𝐄𝐱𝐩}\left(𝓡_{𝝂^1}𝓢_{𝝂^2}t/4\right)\mathrm{}& \\ \hfill \mathrm{}\left(𝐏_{00}𝓚_{22}𝐏_{00}\right)|𝝈_{}(0),& \end{array}$$
(36)
where $`𝝈_{}(0)=𝝈_{}^1𝝈_{}^2`$ and (letting $`𝐏_{00}^4`$ be the $`2\times 2`$ identity tensored with itself $`4`$ times)

$$\begin{array}{cc}\hfill \mathrm{𝐄𝐱𝐩}\left(𝓢_{𝝂^1}𝓡_{𝝂^2}t/4\right)=𝐏_{00}^4+\frac{\mathrm{sin}\left(𝝂^1𝝂^2t/4\right)}{𝝂^1𝝂^2}\left(𝓢_{𝝂^1}𝓡_{𝝂^2}\right)+& \\ \hfill \mathrm{}\frac{1\mathrm{cos}\left(𝝂^1𝝂^2t/4\right)}{𝝂^1^2𝝂^2^2}\left(𝓢_{𝝂^1}𝓡_{𝝂^2}\right)^2& \end{array}$$
(37)
with an almost identical expression for $`\mathrm{𝐄𝐱𝐩}\left(𝓡_{𝝂^1}𝓢_{𝝂^2}t/4\right)`$. The square in the last term may be evaluated by combining Eq. (27) with

$$𝓢_𝐱^2=|𝐱𝐱|+𝐱^2𝐄_{00}𝐄_{00}.$$
(38)
Because $`𝓢_{𝝂^1}𝓡_{𝝂^2}`$ and $`𝓡_{𝝂^1}𝓢_{𝝂^2}`$ are mutually annihilating, the product of their exponentials expands to only five terms, two pairs of which have identical trigonometric coefficients. On pulling the right-hand column operator “$`|`$” back to the left in Eq. (36), essentially reversing what we did to derive the differential version in Eq. (IV), we obtain the integrated equation of motion we have been seeking:

$$\begin{array}{cc}\hfill 𝝈_{}(t)=𝝈_{}^1𝝈_{}^2+\mathrm{sin}\left(𝝂^1𝝂^2t/4\right)\left[[𝝈_{}^1𝝈_{}^2,\widehat{𝝂}^1\widehat{𝝂}^2]\right]& \\ \hfill \mathrm{}\left(1\mathrm{cos}\left(𝝂^1𝝂^2t/4\right)\right)\left[[\left[[𝝈_{}^1𝝈_{}^2,\widehat{𝝂}^1\widehat{𝝂}^2]\right],\widehat{𝝂}^1\widehat{𝝂}^2]\right].\end{array}$$
(39)
Equations (3031) tell us (more or less) what the geometric interpretation of the two particle commutator is, so we turn our attention to the double commutator. Geometric algebra shows that the double commutator of tensor products of the corresponding traceless $`2\times 2`$ Hermitian matrices reduces to

$$\begin{array}{c}2[[𝐀𝐂,𝐁𝐃],𝐁𝐃]\hfill \\ \hfill \begin{array}{cc}\hfill =& [𝐀|𝐁\left(𝐏_{00}[𝐂,𝐃]\right)+\left([𝐀,𝐁]𝐏_{00}\right)𝐂|𝐃,𝐁𝐃]\hfill \\ \hfill =& 𝐀|𝐁𝐁[[𝐂,𝐃],𝐃]+[[𝐀,𝐁],𝐁]𝐃𝐂|𝐃,\hfill \end{array}\end{array}$$
(40)
so we will be done once we figure out what the double commutator of $`2\times 2`$ matrices is. On expanding the commutator and using the fact that for such matrices the anticommutator satisfies $`\mathrm{𝐀𝐁}+\mathrm{𝐁𝐀}=𝐀|𝐁𝐏_{00}`$, we get (including the real analogs):

$$\begin{array}{cc}\hfill [[[[𝐚,𝐛]],𝐛]]& [[𝐀,𝐁],𝐁]=\mathrm{𝐀𝐁}^2\mathrm{\hspace{0.17em}2}\mathrm{𝐁𝐀𝐁}+𝐁^2𝐀\hfill \\ \hfill =& 𝐀𝐁^2\mathrm{\hspace{0.17em}2}(\mathrm{𝐁𝐀}+\mathrm{𝐀𝐁}\mathrm{𝐀𝐁})𝐁\hfill \\ \hfill =& 2\left(𝐀𝐁^2𝐀|𝐁𝐁\right)𝐚𝐛^2𝐚|𝐛𝐛.\hfill \end{array}$$
(41)
Back-substitution of this and the corresponding expression for $`[[𝐂,𝐃],𝐃]`$ into the preceding equation now yields:

$$\begin{array}{c}[[𝐀𝐂,𝐁𝐃],𝐁𝐃]\hfill \\ \hfill \begin{array}{cc}& =𝐀|𝐁𝐁\left(𝐂𝐃^2𝐂|𝐃𝐃\right)+\left(𝐀𝐁^2𝐀|𝐁𝐁\right)𝐃𝐂|𝐃\hfill \\ & =𝐀|𝐁𝐃^2(𝐁𝐂)+𝐂|𝐃𝐁^2(𝐀𝐃)2𝐀|𝐁𝐂|𝐃(𝐁𝐃).\hfill \end{array}\end{array}$$
(42)
Since no matrix products occur in this expression, we may transliterate directly to the corresponding real expression including the additional terms arising from setting $`a_{00}=c_{00}=1`$ (cf. Eq. (31)):

$$\begin{array}{c}4\left[[\left[[(𝐚+𝐄_{00})(𝐜+𝐄_{00}),𝐛𝐝]\right],𝐛𝐝]\right]=\hfill \\ \hfill \begin{array}{cc}\hfill 𝐚|𝐛𝐝^2(𝐛𝐜)+𝐜|𝐝𝐛^2(𝐚𝐝)\mathrm{}& \\ \hfill 2𝐚|𝐛𝐜|𝐝(𝐛𝐝)\mathrm{}& \\ \hfill +2\left[[𝐛[[𝐜,𝐝]]+[[𝐚,𝐛]]𝐝,𝐛𝐝]\right]& .\hfill \end{array}\end{array}$$
(43)
The first of these additional terms (given on the last line of the above equation) may be evaluated via Eqs. (30), (31) & (41) as indicated below,

$$\begin{array}{cc}\hfill 2\left[[𝐛[[𝐜,𝐝]],𝐛𝐝]\right]=& 𝐛^2\left(𝐄_{00}[[[[𝐜,𝐝]],𝐝]]\right)+\left([[𝐛,𝐛]]𝐄_{00}\right)[[𝐜,𝐝]]|𝐝\hfill \\ \hfill =& 𝐛^2\left(𝐄_{00}(𝐝^2𝐜𝐜|𝐝𝐝)\right),\hfill \end{array}$$
(44)
with an analogous expression for the remaining term. We now put in our previous values $`𝐚=\stackrel{ˇ}{𝝈}_{}^1𝝈_{}^1𝐄_{00}`$, $`𝐛=𝝂^1`$, $`𝐜=\stackrel{ˇ}{𝝈}_{}^2𝝈_{}^2𝐄_{00}`$ and $`𝐝=𝝂^2`$ and expand the result fully to get:

$$\begin{array}{c}4\left[[\left[[𝝈_{}^1𝝈_{}^2,𝝂^1𝝂^2]\right],𝝂^1𝝂^2]\right]=\hfill \\ \hfill \begin{array}{cc}\hfill 𝝈_{}^1|𝝂^1𝝂^2^2\left(𝝂^1\stackrel{ˇ}{𝝈}_{}^2\right)+𝝂^1^2𝝈_{}^2|𝝂^2\left(\stackrel{ˇ}{𝝈}_{}^1𝝂^2\right)\mathrm{}& \\ \hfill 2𝝈_{}^1|𝝂^1𝝈_{}^2|𝝂^2\left(𝝂^1𝝂^2\right)\mathrm{}& \end{array}\\ \hfill \begin{array}{cc}\hfill +𝝂^1^2𝝂^2^2\left(𝐄_{00}\stackrel{ˇ}{𝝈}_{}^2+\stackrel{ˇ}{𝝈}_{}^1𝐄_{00}\right)\mathrm{}& \\ \hfill \left(𝝂^1^2𝝈_{}^2|𝝂^2𝐄_{00}𝝂^2+𝝈_{}^1|𝝂^1𝝂^2^2𝝂^1𝐄_{00}\right).& \end{array}\end{array}$$
(45)
Finally, we divide through by $`𝝂^1^2𝝂^2^2`$ to get the normalized “vectors” $`\widehat{𝝂}^1`$ and $`\widehat{𝝂}^2`$, replace the $`𝝈_{}`$’s by $`\stackrel{ˇ}{𝝈}_{}`$’s inside the traces (which doesn’t change their values) and recombine terms to obtain:

$$\begin{array}{c}4\left[[\left[[𝝈_{}^1𝝈_{}^2,\widehat{𝝂}^1\widehat{𝝂}^2]\right],\widehat{𝝂}^1\widehat{𝝂}^2]\right]=\hfill \\ \hfill \begin{array}{cc}& \left(\stackrel{ˇ}{𝝈}_{}^1\stackrel{ˇ}{𝝈}_{}^1|\widehat{𝝂}^1\widehat{𝝂}^1\right)\left(\stackrel{ˇ}{𝝈}_{}^2|\widehat{𝝂}^2\widehat{𝝂}^2+𝐄_{00}\right)\mathrm{}\hfill \\ \hfill +& \left(\stackrel{ˇ}{𝝈}_{}^1|\widehat{𝝂}^1\widehat{𝝂}^1+𝐄_{00}\right)\left(\stackrel{ˇ}{𝝈}_{}^2\stackrel{ˇ}{𝝈}_{}^2|\widehat{𝝂^2}\widehat{𝝂}^2\right).\hfill \end{array}\end{array}$$
(46)
This has a fairly simple interpretation: The rejection of $`\stackrel{ˇ}{𝝈}_{}^1`$ from $`\widehat{𝝂}^1`$ is tensored with the projection of $`\stackrel{ˇ}{𝝈}_{}^2`$ onto $`\widehat{𝝂}^2`$ (plus the usual scalar part of $`1`$) and added to the projection of $`\stackrel{ˇ}{𝝈}_{}^1`$ onto $`\widehat{𝝂}^1`$ (plus the scalar part) tensored with the rejection of $`\stackrel{ˇ}{𝝈}_{}^2`$ from $`\widehat{𝝂}^2`$. This should be contrasted with the single commutator (Eq. (30)), wherein the inner and outer products of each pair are tensored together both ways and added.

## V Meta-Metamorphosis

In the previous section we have shown how the Hamiltonians usually assumed for quantum computing with qubits can be integrated entirely within the real domain. For a single qubit, the results could be interpreted as a simple Bloch vector rotation. With a bi-axial interaction between two qubits, we also found that that the integrated expression had a reasonably nice geometric interpretation. Algebraically, however, it is usually easier to integrate in the Hermitian domain, simply because Hermitian matrices are easier to diagonalize. In this section, therefore, we shall derive formulae by which matrix representations of general superoperators can be translated from the Hermitian into the real domain, along with some specific examples of their utility.

Using the identity given in Eq. (24) together with our operator sum expressions for $`𝒰`$ and its inverse (Eqs. (18) & (19)), it is straightforward to show that an arbitrary superoperator $`𝒮`$ with matrix representation $`𝓢`$ acting on $`𝝆`$ via $`𝝆𝓢|𝝆`$ transforms into the real domain according to

$$𝓢\stackrel{𝒰}{}2^N𝓠\left(\underset{m^{}=0}{\overset{M}{}}𝐏_{M,(Mm^{})}𝐏_{m^{},0}\right)𝓢\left(\underset{m=0}{\overset{M}{}}𝐏_{M,(Mm)}𝐏_{m,0}\right)\overline{𝓠},$$
(47)
where $`𝓠\mathrm{𝐃𝐢𝐚𝐠}(|𝐐^N)`$. The superoperator $`𝒮`$ may be written in operator sum form versus the basis of elementary matrices Havel (2003) as

$$𝓢|𝝆=|\underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}s_k\mathrm{}^{ij}𝐄_{ki}𝝆𝐄_j\mathrm{}=\underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}s_k\mathrm{}^{ij}\left(𝐄_\mathrm{}j𝐄_{ki}\right)|𝝆.$$
(48)
On substituting the second of these equations (sans $`|𝝆`$) into the first and rearranging things a bit, the transformed superoperator becomes

$$2^N𝓠\left(\underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}s_k\mathrm{}^{ij}\underset{m,m^{}=0}{\overset{M,M}{}}\left(𝐏_{M,(Mm^{})}𝐄_\text{}\mathrm{}j𝐏_{M,(Mm)}\right)(𝐏_{m^{},0}𝐄_{ki}𝐏_{m,0})\right)\overline{𝓠}.$$
(49)
Unfortunately, because each factor in the above Kronecker product depends on both indices $`m`$ and $`m^{}`$, this cannot be regarded as a transformation of the elementary matrix basis into the real domain.

Better insight can be obtained by looking at how the Choi matrix $`𝑪𝒉𝒐𝒊(𝓢)`$ of the superoperator transforms. This may be obtained from the propagating matrix $`𝓢`$ simply by replacing Kronecker products of the elementary matrices by dyadic products of the corresponding columnized basis, but in the opposite order Havel (2003). This task is facilitated by expressing the left- and right-multiplication by diagonal matrices as a Hadamard product, using the well-known formula Havel et al. (2001)

$$\mathrm{𝐃𝐢𝐚𝐠}(𝐚)𝐗\mathrm{𝐃𝐢𝐚𝐠}^{}(𝐛)=(\mathrm{𝐚𝐛}^{})𝐗,$$
(50)
which is essentially a special case of Eq. (32). This allows the transformed superoperator to be rewritten as

$$\begin{array}{c}2^N\left(|𝐐^N𝐐^N|\right)\underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}s_k\mathrm{}^{ij}\mathrm{}\hfill \\ \hfill \mathrm{}\underset{m,m^{}=0}{\overset{M,M}{}}(\left(𝐏_{M,(Mm^{})}𝐄_\text{}\mathrm{}j𝐏_{M,(Mm)}\right)(𝐏_{m^{},0}𝐄_{ki}𝐏_{m,0})).\end{array}$$
(51)
The advantage of this form is that the Hadamard product commutes with the $`𝑪𝒉𝒐𝒊`$ operator (since it rearranges the entries of the product’s operands identically), giving us

$`𝑪𝒉𝒐𝒊(𝓢)\stackrel{𝒰}{}`$ $`2^N𝑪𝒉𝒐𝒊\left(|𝐐^N𝐐^N|\right){\displaystyle \underset{i,j=0}{\overset{M,M}{}}}{\displaystyle \underset{k,\mathrm{}=0}{\overset{M,M}{}}}s_k\mathrm{}^{ij}\mathrm{}`$
$`\mathrm{}{\displaystyle \underset{m,m^{}=0}{\overset{M,M}{}}}𝑪𝒉𝒐𝒊(\left(𝐏_{M,(Mm^{})}𝐄_\text{}\mathrm{}j𝐏_{M,(Mm)}\right)(𝐏_{m^{},0}𝐄_{ki}𝐏_{m,0}))`$
$`=`$ $`2^N\left(\overline{𝐐}^N𝐐^N\right){\displaystyle \underset{i,j=0}{\overset{M,M}{}}}{\displaystyle \underset{k,\mathrm{}=0}{\overset{M,M}{}}}s_k\mathrm{}^{ij}\mathrm{}`$
$`\mathrm{}{\displaystyle \underset{m,m^{}=0}{\overset{M,M}{}}}\left(|𝐏_{m^{},0}𝐄_{ki}𝐏_{m,0}𝐏_{M,(Mm^{})}𝐄_\text{}\mathrm{}j𝐏_{M,(Mm)}|\right)`$
$`=`$ $`2^N\left(\overline{𝐐}^N𝐐^N\right){\displaystyle \underset{i,j=0}{\overset{M,M}{}}}{\displaystyle \underset{k,\mathrm{}=0}{\overset{M,M}{}}}s_k\mathrm{}^{ij}{\displaystyle \underset{m,m^{}=0}{\overset{M,M}{}}}\left(𝐏_{m,0}𝐏_{m^{},0}\right)|𝐄_{ki}\mathrm{}`$
$`\mathrm{}𝐄_\text{}\mathrm{}j|\left(𝐏_{M,(Mm)}𝐏_{M,(Mm^{})}\right)`$ (52)
$`=`$ $`2^N\left(\overline{𝐐}^N𝐐^N\right){\displaystyle \underset{m,m^{}=0}{\overset{M,M}{}}}\left(𝐏_{m,0}𝐏_{m^{},0}\right){\displaystyle \underset{i,j=0}{\overset{M,M}{}}}{\displaystyle \underset{k,\mathrm{}=0}{\overset{M,M}{}}}s_k\mathrm{}^{ij}\mathrm{}`$
$`\mathrm{}\left(𝐄_{ij}𝐄_k\mathrm{}\right)\left(𝐏_{M,(Mm)}𝐏_{M,(Mm^{})}\right)`$
$``$ $`2^N\left(\overline{𝐐}^N𝐐^N\right){\displaystyle \underset{m=0}{\overset{M^{}}{}}}𝐏_{m,0}𝑪𝒉𝒐𝒊(𝓢)𝐏_{M^{},(M^{}m)},`$
where $`M^{}(M+1)^21=2^{2N}1`$ and we have used the relation $`|𝐄_{ki}𝐄_\mathrm{}j|=(𝐞_i𝐞_k)(𝐞_j𝐞_{\mathrm{}})^{}=𝐄_{ij}𝐄_k\mathrm{}`$.

In other words, the Choi matrix of a superoperator maps into the real domain much like a density matrix on twice as many qubits. In fact we can write the *real* transformation matrix $`𝓣`$, which acts on the real density matrix as $`𝓣|𝝈_{}`$, in the following compact form:

$$𝓣=2^N𝑪𝒉𝒐𝒊\left(𝒰\left(𝑪𝒉𝒐𝒊(𝓢)\right)\right)\overline{𝓠}^{\mathrm{\hspace{0.17em}2}}.$$
(53)
Turning this around, we also find that we can express the Choi matrix of $`𝒮`$ in the Pauli basis as

$$\begin{array}{cc}& 2^N𝑪𝒉𝒐𝒊(𝓢)=𝒰^1\left(𝑪𝒉𝒐𝒊(𝓣𝓠^2)\right)\hfill \\ \hfill & \underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}𝒰^1\left(t_k\mathrm{}^{ij}𝑪𝒉𝒐𝒊\left((𝐄_\mathrm{}j𝐄_{ki})𝓠^2\right)\right)\hfill \\ \hfill =& \underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}t_k\mathrm{}^{ij}𝒰^1\left(𝑪𝒉𝒐𝒊(𝐄_\mathrm{}j𝐄_{ki})𝑪𝒉𝒐𝒊(|\mathbf{\hspace{0.17em}1\hspace{0.17em}1}^{}𝐐^N𝐐^N|)\right)\hfill \\ \hfill =& \underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}t_k\mathrm{}^{ij}𝒰^1(𝐄_{ij}𝐄_k\mathrm{})((𝐐^N𝐐^N)\left(\mathrm{𝟏}\mathrm{𝟏}^{}\right)))\hfill \\ \hfill =& 2^{2N}\left(\left(𝐐^N𝐐^N\right)\left(\mathrm{𝟏}\mathrm{𝟏}^{}\right)\right)\underset{i,j=0}{\overset{M,M}{}}\underset{k,\mathrm{}=0}{\overset{M,M}{}}t_k\mathrm{}^{ij}\left(𝐏_{ij}𝐏_k\mathrm{}\right),\hfill \end{array}$$
(54)
where $`\mathrm{𝟏}`$ denotes a column vector of ones of the appropriate size.

It is time for our examples! We shall begin with operator sums for single qubit rotations, and go on to show how rotations about the $`𝝈_3`$ axis as well as the $`𝝈_3^1𝝈_3^2`$ interaction between two qubits can be compactly described in the real domain using Hadamard products. We close by showing that this description also extends quite nicely to $`𝝈_3`$ dephasing as well as nonunital relaxation back towards a nonrandom equilibrium state (i.e. $`T_2`$ and $`T_1`$ relaxation in NMR parlance).

Consider the Choi matrix of the propagator which rotates a single qubit (Eqs. (2527) above):

$`𝑪𝒉𝒐𝒊\left(𝐏_{00}𝐏_{00}+\mathrm{sin}(𝝂t/2)𝓡_{\widehat{𝝂}}+(1\mathrm{cos}(𝝂t/2))𝓡_{\widehat{𝝂}}^{\mathrm{\hspace{0.17em}2}}\right)`$
$`=`$ $`\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 1\end{array}\right]+\mathrm{sin}(𝝂t/2)\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& \widehat{\nu }_{11}& \widehat{\nu }_{10}\\ 0& \widehat{\nu }_{11}& 0& \widehat{\nu }_{10}\\ 0& \widehat{\nu }_{01}& \widehat{\nu }_{01}& 0\end{array}\right]+\mathrm{}`$ (55)
$`\mathrm{}(1\mathrm{cos}(𝝂t/2))\left[\begin{array}{cccc}0& 0& 0& \widehat{\nu }_{10}^2\widehat{\nu }_{11}^2\\ 0& 0& \widehat{\nu }_{10}\widehat{\nu }_{01}& \widehat{\nu }_{01}\widehat{\nu }_{11}\\ 0& \widehat{\nu }_{10}\widehat{\nu }_{01}& 0& \widehat{\nu }_{01}\widehat{\nu }_{11}\\ \widehat{\nu }_{01}^2\widehat{\nu }_{11}^2& \widehat{\nu }_{10}\widehat{\nu }_{11}& \widehat{\nu }_{10}\widehat{\nu }_{11}& \widehat{\nu }_{01}^2\widehat{\nu }_{10}^2\end{array}\right],`$
where the “hat” on the $`\nu `$’s indicates normalization by $`𝝂`$. A general operator sum can be derived from this matrix (as well as by expanding the implied commutators via Eq. (23)), but since this is a bit involved we shall restrict ourselves to rotations by an angle $`\vartheta =𝝂t/2`$ about the $`𝗑`$, $`𝗒`$ or $`𝗓`$ coordinate axes. On substituting $`\widehat{\nu }_{10}=1`$ and $`\widehat{\nu }_{01}=\widehat{\nu }_{11}=0`$ into the above Choi matrix we obtain the following singular value decomposition:

$$\begin{array}{cc}\hfill \left[\begin{array}{cccc}1& 0& 0& \mathrm{cos}(\vartheta )\\ 0& 0& 0& \mathrm{sin}(\vartheta )\\ 0& 0& 0& \mathrm{sin}(\vartheta )\\ 1& 0& 0& \mathrm{cos}(\vartheta )\end{array}\right]=& \left[\begin{array}{cc}\mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\\ \mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\\ \mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\\ \mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\end{array}\right]\mathrm{}\hfill \\ & \mathrm{}\left[\begin{array}{cc}\mathrm{cos}(\vartheta /2)& 0\\ 0& \mathrm{sin}(\vartheta /2)\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 1\\ 1& 0& 0& 1\end{array}\right],\hfill \end{array}$$
(56)
as may be readily verified using the usual half-angle formulae. The corresponding operator sum for rotation by $`\vartheta `$ about the $`𝗑`$-axis is simply:

$$\begin{array}{cc}\hfill 𝒰\left(𝒆^{ı(\vartheta /2)𝐏_{10}}𝝆𝒆^{ı(\vartheta /2)𝐏_{10}}\right)=& \mathrm{cos}(\vartheta /2)\left[\begin{array}{cc}\mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\\ \mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]+\mathrm{}\hfill \\ & \mathrm{sin}(\vartheta /2)\left[\begin{array}{cc}\mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\\ \mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].\hfill \end{array}$$
(57)
In a similar fashion, it can be shown that the operator sum for a $`𝗒`$-rotation is:

$$\begin{array}{cc}\hfill 𝒰\left(𝒆^{ı(\vartheta /2)𝐏_{01}}𝝆𝒆^{ı(\vartheta /2)𝐏_{01}}\right)=& \mathrm{cos}(\vartheta /2)\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]\left[\begin{array}{cc}\mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\\ \mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\end{array}\right]+\mathrm{}\hfill \\ & \mathrm{sin}(\vartheta /2)\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]\left[\begin{array}{cc}\mathrm{sin}(\vartheta /2)& \mathrm{cos}(\vartheta /2)\\ \mathrm{cos}(\vartheta /2)& \mathrm{sin}(\vartheta /2)\end{array}\right].\hfill \end{array}$$
(58)
For a $`𝗓`$-rotation, on the other hand, the Choi matrix turns out to be rank $`4`$ with singular value decomposition:

$$\begin{array}{cc}& \left[\begin{array}{cccc}1& 0& 0& \mathrm{cos}(\vartheta )\\ 0& 0& \mathrm{sin}(\vartheta )& 0\\ 0& \mathrm{sin}(\vartheta )& 0& 0\\ \mathrm{cos}(\vartheta )& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 1& 0& 0& 1\end{array}\right]\mathrm{}\hfill \\ & \mathrm{}\left[\begin{array}{cccc}\mathrm{cos}^2(\vartheta /2)& 0& 0& 0\\ 0& \mathrm{cos}(\vartheta /2)\mathrm{sin}(\vartheta /2)& 0& 0\\ 0& 0& \mathrm{cos}(\vartheta /2)\mathrm{sin}(\vartheta /2)& 0\\ 0& 0& 0& \mathrm{sin}^2(\vartheta /2)\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 1& 0& 0& 1\end{array}\right].\hfill \end{array}$$
(59)
This corresponds to the operator sum

$$\begin{array}{cc}\hfill 𝒰\left(𝒆^{ı(\vartheta /2)𝐏_{11}}𝝆𝒆^{ı(\vartheta /2)𝐏_{11}}\right)=& \mathrm{cos}^2(\vartheta /2)𝝈_{}+\mathrm{sin}^2(\vartheta /2)𝐏_{11}𝝈_{}𝐏_{11}+\mathrm{}\hfill \\ & \mathrm{}ı\mathrm{cos}(\vartheta /2)\mathrm{sin}(\vartheta /2)\left(𝐏_{10}𝝈_{}𝐏_{01}+𝐏_{01}𝝈_{}𝐏_{10}\right),\hfill \end{array}$$
(60)
which has the pleasant feature that the trigonometric functions occur as scalar factors in each term and not embedded in the operators. This enables us to use Eq. (50) to rewrite it in terms of the Hadamard product as follows:

$$\begin{array}{cc}\hfill 𝒰\left(𝒆^{ı(\vartheta /2)𝐏_{11}}𝝆𝒆^{ı(\vartheta /2)𝐏_{11}}\right)& =\left[\begin{array}{cc}1& \mathrm{cos}(\vartheta )\\ \mathrm{cos}(\vartheta )& 1\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]+\mathrm{}\hfill \\ \hfill \mathrm{}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]& \left(\left[\begin{array}{cc}0& \mathrm{sin}(\vartheta )\\ \mathrm{sin}(\vartheta )& 0\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}\\ \sigma _{10}& \sigma _{11}\end{array}\right]\right)\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right].\hfill \end{array}$$
(61)
As our next example, wherein the Hadamard product enables even greater simplifications, consider an “Ising-type” interaction between two qubits of the form $`𝝈_3𝝈_3=𝐏_{33}`$, which is also known as “weak scalar coupling” in NMR Ernst et al. (1987). An operator sum expression for this could be obtained by expanding the general formula given in Eq. (39), but because this Hamiltonian is again diagonal in the $`𝝈_3`$ eigenbasis a simpler expression can be obtained directly starting from the diagonal matrix of the corresponding propagator, i.e.

$$\begin{array}{cc}\hfill |𝒆^{ı𝐏_{33}\pi Jt/2}𝝆𝒆^{ı𝐏_{33}\pi Jt/2}=& \left(𝒆^{ı𝐏_{33}\pi Jt/2}𝒆^{ı𝐏_{33}\pi Jt/2}\right)|𝝆\hfill \\ \hfill =& \mathrm{𝐃𝐢𝐚𝐠}\left(|𝐉(t)\right)|𝝆𝓙(t)|𝝆,\hfill \end{array}$$
(62)
where

$$𝐉(t)\left[\begin{array}{cccc}1& e^{ı\pi Jt}& e^{ı\pi Jt}& 1\\ e^{ı\pi Jt}& 1& 1& e^{ı\pi Jt}\\ e^{ı\pi Jt}& 1& 1& e^{ı\pi Jt}\\ 1& e^{ı\pi Jt}& e^{ı\pi Jt}& 1\end{array}\right].$$
(63)
On transforming this into the real domain via Eq. (47), we find that $`𝓙`$ has been converted into the sum of the diagonal matrix

$$𝓒(t)\mathrm{𝐃𝐢𝐚𝐠}\left(\text{}\left[\begin{array}{cccc}1& \mathrm{cos}(\pi Jt)& \mathrm{cos}(\pi Jt)& 1\\ \mathrm{cos}(\pi Jt)& 1& 1& \mathrm{cos}(\pi Jt)\\ \mathrm{cos}(\pi Jt)& 1& 1& \mathrm{cos}(\pi Jt)\\ 1& \mathrm{cos}(\pi Jt)& \mathrm{cos}(\pi Jt)& 1\end{array}\right]\text{}\right),$$
(64)
and the anti-diagonal matrix

$$𝐏_{15,0}𝓢(t)\mathrm{𝐃𝐢𝐚𝐠}\left(\text{}\left[\begin{array}{cccc}0& \mathrm{sin}(\pi Jt)& \mathrm{sin}(\pi Jt)& 0\\ \mathrm{sin}(\pi Jt)& 0& 0& \mathrm{sin}(\pi Jt)\\ \mathrm{sin}(\pi Jt)& 0& 0& \mathrm{sin}(\pi Jt)\\ 0& \mathrm{sin}(\pi Jt)& \mathrm{sin}(\pi Jt)& 0\end{array}\right]\text{}\right),$$
(65)
where the (self-inverse) left factor of $`𝐏_{15,0}=𝝈_1^4`$ simply reverses the order of the rows. The nonzero entries of the Choi matrix of $`𝓒(t)+𝓢(t)`$ turn out to comprise two $`4\times 4`$ blocks along the diagonal, which are exactly the two matrices above, i.e.

$$\begin{array}{cc}\hfill 𝓒(t)=& \mathrm{𝐃𝐢𝐚𝐠}\left(|𝓔_𝓒𝑪𝒉𝒐𝒊(𝓒(t))𝓔_𝓒^{}\right)\hfill \\ \hfill \text{and}𝐏_{15,0}𝓢(t)=& \mathrm{𝐃𝐢𝐚𝐠}\left(|𝓔_𝓢𝑪𝒉𝒐𝒊(𝓢(t))𝓔_𝓢^{}\right),\hfill \end{array}$$
(66)
wherein

$$𝓔_𝓒\underset{i=0}{\overset{3}{}}𝐞_i(𝐞_i𝐞_i)^{}\text{and}𝓔_𝓢\underset{i=0}{\overset{3}{}}𝐞_i(𝐞_i𝐞_{3i})^{}$$
(67)
project out the rows / columns of their respective blocks.

Thus we can obtain the desired operator sum representation by computing the eigenvalues and eigenvectors of the $`4\times 4`$ symmetric matrices

$$𝐂(t)𝓔_𝓒𝑪𝒉𝒐𝒊(𝓒(t))𝓔_𝓒^{}\text{and}𝐒(t)𝓔_𝓢𝑪𝒉𝒐𝒊(𝓢(t))𝓔_𝓢^{},$$
(68)
letting the operators’ matrices be the diagonal / anti-diagonal matrices formed from the entries of these eigenvectors, and multiplying each term in the sum by the corresponding eigenvalue. The results are

$$𝓒(t)|𝝈_{}=|\frac{1}{2}\left(1+\mathrm{cos}(\pi Jt)\right)𝝈_{}+\frac{1}{2}\left(1\mathrm{cos}(\pi Jt)\right)𝐏_{33}𝝈_{}𝐏_{33}$$
(69)
and

$$\begin{array}{cc}\hfill 𝓢(t)|𝝈_{}=\frac{1}{2}\mathrm{sin}(\pi Jt)|𝐏_{30}\mathrm{𝐃𝐢𝐚𝐠}\left([1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},1]\right)𝝈_{}\mathrm{𝐃𝐢𝐚𝐠}\left([1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},1]\right)𝐏_{30}& \mathrm{}\hfill \\ \hfill \mathrm{}𝐏_{30}\mathrm{𝐃𝐢𝐚𝐠}\left([1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1}]\right)𝝈_{}\mathrm{𝐃𝐢𝐚𝐠}\left([1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1}]\right)𝐏_{30}& .\hfill \end{array}$$
(70)
By using Eq. (50) to replace these operator sums by Hadamard products and taking advantage of the symmetry of $`𝐒(t)`$, however, we can obtain an even simpler expression, namely

$$\left(𝓒(t)+𝓢(t)\right)|𝝈_{}=|𝐂(t)𝝈_{}𝐒(t)(𝐏_{30}𝝈_{}𝐏_{30}).$$
(71)
Finally, we show how one can also use Hadamard products with the real density matrix to describe simple relaxation processes, in a manner similar to that described in Havel et al. (2001) for the usual Hermitian density matrix. For a single qubit undergoing $`T_1`$ (dissipation) and $`T_2`$ (decoherence) relaxation, the time derivative is given by:

$$_t𝝈_{}(t)=𝐑𝝈_{}(t)\left[\begin{array}{cc}0& 1/T_2\\ 1/T_2& 1/T_1\end{array}\right]\left[\begin{array}{cc}1& \sigma _{01}(t)\\ \sigma _{10}(t)& \sigma _{11}(t)\end{array}\right].$$
(72)
Assuming that these relaxation processes are uncorrelated, this can immediately be extended to any number of qubits using the fact that the Hadamard product satisfies the mixed product formula with the Kronecker product (Eq. (15)). In the case of two qubits relaxing with Hadamard relaxation matrices $`𝐑^1`$, $`𝐑^2`$, for example, we obtain

$$\begin{array}{cc}\hfill _t𝝈_{}(t)=& \left(𝐑^1(\mathrm{𝟏𝟏}^{})+(\mathrm{𝟏𝟏}^{})𝐑^2\right)𝝈_{}(t)\hfill \\ \hfill =& \left[\begin{array}{cccc}0& 1/T_2^2& 1/T_2^1& 1/T_2^1+1/T_2^2\\ 1/T_2^2& 1/T_1^2& 1/T_2^1+1/T_2^2& 1/T_2^1+1/T_1^2\\ 1/T_2^1& 1/T_2^1+1/T_2^2& 1/T_1^1& 1/T_1^1+1/T_2^2\\ 1/T_2^1+1/T_2^2& 1/T_2^1+1/T_1^2& 1/T_1^1+1/T_2^2& 1/T_1^1+1/T_1^2\end{array}\right]𝝈_{}(t).\hfill \end{array}$$
(73)
where $`\mathrm{𝟏}`$ is a $`4\times 1`$ vector of $`1`$’s. The fact that uncorrelated $`T_1`$ as well as $`T_2`$ relaxation can be extended so easily to multiple spins in this way is actually a significant advantage of the real density matrix over the Hermitian, since in the latter case the diagonal terms are mixtures of terms decaying at differing rates, substantially complicating their treatment via Hadamard products Havel et al. (2001).

When correlations are present, however, these advantages are largely lost, since then the off-diagonal entries of the real density matrix consist of mixtures of terms with differing decay rates (fortunately, $`T_1`$ relaxation is usually largely uncorrelated Ernst et al. (1987)). Let us work through the case of two qubits in detail, assuming for simplicity that the $`T_2`$ relaxation processes at the two qubits are totally correlated and have the same rate $`1/T_2`$, as for example in an NMR gradient-diffusion experiment Havel et al. (2001). In this case the Hadamard relaxation matrix for the Hermitian density matrix has the form

$$𝐑=\frac{1}{T_2}\left[\begin{array}{cccc}0& 1& 1& 4\\ 1& 0& 0& 1\\ 1& 0& 0& 1\\ 4& 1& 1& 0\end{array}\right],$$
(74)
and the corresponding $`16\times 16`$ diagonal relaxation superoperator

$$𝓡=\mathrm{𝐃𝐢𝐚𝐠}\left(|𝐑\right)$$
(75)
is easily exponentiated into a diagonal matrix of survival probabilities for the entries of the (traceless part of the) Hermitian density matrix. In this case, however, it turns out to be almost as easy, but more revealing, to convert $`𝓡`$ into the real domain and perform the integration there. The result of the first step is

$$𝓡\stackrel{𝒰}{}\frac{1}{T_2}\left[\begin{array}{cccccccccccccccc}0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& 2& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& 2& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}2}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}1}& \mathrm{\hspace{0.25em}0}\\ 0& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}& \mathrm{\hspace{0.25em}0}\end{array}\right],$$
(76)
which again has nonzero entries only on its diagonal and its anti-diagonal. Representing the action of this matrix on $`|𝝈_{}`$ as a sum of a Hadamard product and a Hadamard product coupled with row/column inversion as in Eq. (71), we obtain the real equation of motion:

$$T_2_t𝝈_{}(t)=\left[\begin{array}{cccc}0& 1& 1& 2\\ 1& 0& 2& 1\\ 1& 2& 0& 1\\ 2& 1& 1& 0\end{array}\right]𝝈_{}(t)+\left[\begin{array}{cccc}0& 0& 0& 2\\ 0& 0& 2& 0\\ 0& 2& 0& 0\\ 2& 0& 0& 0\end{array}\right]\left(𝐏_{30}𝝈_{}(t)𝐏_{30}\right)$$
(77)
It may readily be verified that the operations on $`𝝈_{}`$ which occur in the two terms of this expression commute, and hence this equation can be integrated by exponentiating them separately. With the first term this leads to a simple Hadamard (entrywise) exponential Havel et al. (2001), namely

$$𝐃(t)\mathrm{𝐄𝐱𝐩}_{}\left(\frac{t}{T_2}\left[\begin{array}{cccc}0& 1& 1& 2\\ 1& 0& 2& 1\\ 1& 2& 0& 1\\ 2& 1& 1& 0\end{array}\right]\right)=\left[\begin{array}{cccc}1& e^{t/T_2}& e^{t/T_2}& e^{2t/T_2}\\ e^{t/T_2}& 1& e^{2t/T_2}& e^{t/T_2}\\ e^{t/T_2}& e^{2t/T_2}& 1& e^{t/T_2}\\ e^{2t/T_2}& e^{t/T_2}& e^{t/T_2}& 1\end{array}\right].$$
(78)
To exponentiate the second term, we note that the Hadamard product is with $`2𝐏_{03}=2𝝈_2𝝈_2`$ and resort briefly to superoperators in order to simplify the exponential as follows:

$$\begin{array}{cc}& \mathrm{𝐄𝐱𝐩}\left((2t/T_2)\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{03})(𝐏_{30}𝐏_{30})\right)\hfill \\ \hfill =& \underset{k=0}{\overset{\mathrm{}}{}}\left((2t/T_2)\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{03})(𝐏_{30}𝐏_{30})\right)^k/k!\hfill \\ \hfill =& \begin{array}{cc}\hfill 𝐏_{00}𝐏_{00}+\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{30})\underset{k=1}{\overset{\mathrm{}}{}}\frac{(2t/T_2)^{2k}}{2k!}\mathrm{}& \\ \hfill \mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{03})& (𝐏_{30}𝐏_{30})\underset{k=0}{\overset{\mathrm{}}{}}\frac{(2t/T_2)^{2k+1}}{(2k+1)!}\hfill \end{array}\hfill \\ \hfill =& \begin{array}{cc}\hfill 𝐏_{00}𝐏_{00}+\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{30})(\mathrm{cosh}(2t/T_2)\mathrm{\hspace{0.17em}1})\mathrm{}& \\ \hfill \mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{03})& (𝐏_{30}𝐏_{30})\mathrm{sinh}(2t/T_2).\hfill \end{array}\hfill \end{array}$$
(79)
This derivation relies on the facts that $`𝐏_{30}𝐏_{30}`$ squares to the identity, $`\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{03})`$ squares to $`\mathrm{𝐃𝐢𝐚𝐠}(|𝐏_{30})`$, and each commutes with the other.

Going back to operator sum notation and abbreviating $`𝝈_{}𝝈_{}(0)`$, we thus obtain in all:

$`𝝈_{}(t)=\begin{array}{cc}\hfill 𝐃(t)(𝝈_{}+(\mathrm{cosh}(2t/T_2)1)𝐏_{30}𝝈_{}\mathrm{}& \\ \hfill \mathrm{sinh}& (2t/T_2)𝐏_{03}(𝐏_{30}𝝈_{}𝐏_{30}))\hfill \end{array}`$
$`=`$ $`\left[\begin{array}{cccc}\sigma _{00}& \sigma _{01}e^{t/T_2}& \sigma _{02}e^{t/T_2}& \mathrm{\Delta }(t;\sigma _{03},\sigma _{30})\\ \sigma _{10}e^{t/T_2}& \sigma _{11}& \mathrm{\Delta }(t;\sigma _{12},\sigma _{21})& \sigma _{13}e^{t/T_2}\\ \sigma _{01}e^{t/T_2}& \mathrm{\Delta }(t;\sigma _{21},\sigma _{12})& \sigma _{22}& \sigma _{23}e^{t/T_2}\\ \mathrm{\Delta }(t;\sigma _{30},\sigma _{03})& \sigma _{31}e^{t/T_2}& \sigma _{32}e^{t/T_2}& \sigma _{33}\end{array}\right],`$ (80)
where $`\mathrm{\Delta }(t;x,y)(\mathrm{cosh}(2t/T_2)x+\mathrm{sinh}(2t/T_2)y))\mathrm{exp}(2t/T_2)`$. From this we see that the anti-diagonal entries decoher into mixtures with their symmetrically placed opposites in the real density matrix. These mixtures correspond to the real and imaginary parts of the $`\rho _{12}=\overline{\rho }_{21}`$ entries in the Hermitian density matrix, otherwise known as *zero-quantum coherences*, which are immune to correlated noise Ernst et al. (1987).

## VI Epilogue

We have seen that one can, with some effort, do pretty much everything with the real density matrix that one could with the usual Hermitian one. This may be useful as a didactic device, or in calculations with experimental (e.g. NMR) data where it is desirable to keep the experimentally measured values of the observables in sight at all times. This work is also a good demonstration of the power of Choi matrix decompositions as a means of finding operator sum representations of linear superoperators Havel (2003).

Although the Hermitian density matrix is expected to be better suited, by and large, for the purposes of numerical calculations, it is worth emphasizing that for theoretical and/or expository purposes the compact but lucid notation of geometric algebra offers significant advantages over any matrix formalism. In this regard, we point out that Havel et al. (2003) have recently introduced a *parity-even* (rather than reverse-even, aka Hermitian) multi-qubit density operator via geometric algebra, which generalizes the *multi-particle space-time algebra* introduced for isolated systems to open multi-qubit systems. It is our hope that in due course such a geometric formulation may provide new insights into some of the conceptual problems that underlie quantum physics.

The existence of the real density matrix is further of some theoretical interest, since it provides a coordinate ring within which one can study the issues of entanglement and decoherence via invariant theoretic methods Grassl et al. (1998); Makhlin (2002). There are intimate connections between invariant theory and geometric algebra, and it is often easier to automate symbolic computations in an invariant ring than it is at the more abstract level of geometric algebra Sturmfels (1993); Havel (1997, 2001).

###### Acknowledgements.

The author thanks Nicolas Boulant, David Cory and Chris Doran for useful discussions. This work was supported by ARO grants DAAD19-01-1-0519, DAAD19-01-1-0678, by DARPA grant MDA972-01-1-0003, and by a grant from the Cambridge-MIT Institute, Ltd.


format: text

text:
# Color Octet Contribution to High 𝑝_𝑇 J/Ψ Production in pp Collisions at √𝑠 = 500 and 200 GeV at RHIC

## Abstract

We compute $`\frac{d\sigma }{dp_T}`$ of the $`J/\psi `$ production in pp collisions at RHIC at $`\sqrt{s}`$ = 500 and 200 GeV by using both the color octet and singlet models in the framework of non-relativistic QCD. The $`J/\psi `$ we compute here includes the direct $`J/\psi `$ from the partonic fusion processes and the $`J/\psi `$ coming from the radiative decays of $`\chi _J`$’s both in the color octet and singlet channel. The high $`p_T`$ $`J/\psi `$ production cross section is computed within the PHENIX detector acceptance ranges: $`0.35<\eta <0.35`$ and $`1.2<|\eta |<2.4`$, the central electron and forward muon arms. It is found that the color octet contribution to $`J/\psi `$ production is dominant at RHIC energy in comparison to the color singlet contributions. We compare our results with the recent preliminary data obtained by PHENIX detector for the high $`p_T`$ $`J/\psi `$ measurements. While the color singlet model fails to explain the data completely the color octet model is in agreement with the single data point above 2 GeV transverse momentum. A measurement of $`J/\psi `$ production at RHIC in the next run with better statistics will allow us to determine the validity of the color octet model of $`J/\psi `$ production at RHIC energies. This is very important because it is necessary to know the exact mechanism for $`J/\psi `$ production in pp collisions at RHIC if one is to make predictions of $`J/\psi `$ suppression as a signature of quark-gluon plasma. These mechanisms also play an important role in determining the polarized spin structure function of the proton at RHIC.

Recently there has been much progress in understanding of heavy quarkonium production mechanism in high energy hadronic collisions. The two prominent models for heavy quarkonium production mechanism which are quite successful at different collider energies are 1) color singlet model sing  and 2) color octet model oct ; cho1 ; cho2 ; oct2 . The $`J/\psi `$ production measurement by CDF collaboration at Tevatron teva  has ruled out the color singlet model for high $`p_T`$ quarkonium production braaten  in that center of mass energy regime. The color singlet model gives a reasonable description of large $`p_T`$ $`J/\psi `$ production at ISR energies isr . In this context large $`p_T`$ $`J/\psi `$ production data from pp collisions at $`\sqrt{s}`$ = 500 and 200 GeV at RHIC might tell us the relative importance and validity of these two models in this new energy range. In particular a valid $`J/\psi `$ production mechanism at pp collisions will play an important role in determining the mechanism for $`J/\psi `$ suppression at Au-Au collisions which is suggested to be a prominent signature for quark-gluon plasma detection at RHIC nayak ; nayk . These mechanisms also play an important role in determining the polarized spin structure function of the proton.

In the color singlet model the quarkonium is formed as a non relativistic bound state of heavy quark-antiquark pair via static gluon exchange. In this model it is assumed that the Q$`\overline{Q}`$ is produced in the color singlet state at the production point with appropriate spin (S) and orbital angular momentum quantum number (L) to evolve to a bound state $`{}_{}{}^{2S+1}L_{J}^{}`$ with total angular momentum quantum number (J). The relative momentum of the $`Q\overline{Q}`$ pair inside the quarkonium is assumed to be small compared to the mass $`M_Q`$ of the heavy quark so that the $`Q`$ and $`\overline{Q}`$ will not fly apart to form heavy mesons. The non-relativistic non-perturbative wave functions and its derivatives appearing in the color singlet model calculation are either determined from the potential model or taken from experiment.

On the other hand in the color octet model relativistic effects are taken into account which are neglected in the color singlet model. In the color octet model using an effective field theory called non-relativistic QCD (NRQCD), the dynamical gluon enters into Fock state decompositions of the quarkonium states. In NRQCD the expansion is carried out in terms of the relative velocity $`v`$ ($`v^2`$ 0.23 for C$`\overline{C}`$ system and 0.1 for B$`\overline{B}`$ system) of the Q$`\overline{Q}`$ bound state. The NRQCD Lagrangian density is given by:

$$_{NRQCD}=_{light}+_{heavy}+_{correction}$$
(1)
where

$$_{light}=\frac{1}{4}F^{a\mu \nu }F_{\mu \nu }^a+\overline{q}\gamma _\mu D^\mu [A]q,$$
(2)
is the usual QCD Lagrangian for gluons and the light flavors and

$`_{heavy}=\psi ^{}(iD_t+{\displaystyle \frac{D^2}{2M}})\psi +\varphi ^{}(iD_t+{\displaystyle \frac{D^2}{2M}})\varphi `$ (3)
is the leading order heavy quark part and

$`_{correction}={\displaystyle \frac{1}{8M^3}}[\psi ^{}D^4\psi \varphi ^{}D^4\varphi ]`$
$`+{\displaystyle \frac{g}{8M^2}}[\psi ^{}(DEED)\psi +\varphi ^{}(DEED)\varphi ]`$
$`+{\displaystyle \frac{ig}{8M^2}}[\psi ^{}\sigma (D\times EE\times D)\psi +\varphi ^{}\sigma (D\times EE\times D)\varphi ]`$
$`+{\displaystyle \frac{g}{2M}}[\psi ^{}\sigma B\psi \varphi ^{}\sigma B\varphi ]+\mathrm{}`$ (4)
is the correction term. In the above equations $`\psi ,\varphi `$ are the two component Dirac spinors of the heavy quark. The various correction terms in the above equation come from the corrections to kinetic energy, electric field coupling, magnetic field coupling and from double electric scattering etc.. In NRQCD the dynamical gluon enter into the Fock state decompositions of different physical states. The wave functions of S-wave orthoquarkonia state $`|\psi _Q>`$ appear as:

$`|\psi _Q>=O(1)|Q\overline{Q}[^3S_1^{(1)}]>+O(v)|Q\overline{Q}[^3P_J^{(8)}]g>+O(v^2)|Q\overline{Q}[^3S_1^{(1,8)}]gg>+`$
$`O(v^2)|Q\overline{Q}[^1S_0^{(8)}]g>+O(v^2)|Q\overline{Q}[^3D_J^{(1,8)}]gg>+\mathrm{}..`$ (5)
and the wave functions of P-wave orthoquarkonium state $`|\chi _{QJ}>`$:

$`|\chi _{QJ}>=O(1)|Q\overline{Q}[^3P_J^{(1)}]>+O(v)|Q\overline{Q}[^3S_1^{(8)}]g>+\mathrm{}.`$ (6)
In the above equation (1,8) refers to singlet and octet state of the $`Q\overline{Q}`$ pair. After a Q$`\overline{Q}`$ is formed in its color octet state it may absorb a soft gluon and transform into $`|\chi _{QJ}>`$ via eq. (6) and then become $`J/\psi `$ by photon decay or the Q$`\overline{Q}`$ in the color octet state can emit two long wavelength gluons and become a $`J/\psi `$ via eq. (5) and so on. All these low energy interactions are negligible and the non-perturbative matrix elements can be fitted from the experiments or can be determined from the lattice calculations.

Using both mechanisms for heavy quarkonium production we compute the $`p_T`$ distribution of J/$`\psi `$ production at RHIC energy in this paper. The $`J/\psi `$ production includes both direct $`J/\psi `$ and the $`J/\psi `$ coming from the radiative decay of $`\chi _J`$’s: $`\chi _JJ/\psi +\gamma `$. We will calculate the $`J/\psi `$ production in this paper by considering all the parton fusion processes. In the color octet model the non-perturbative matrix elements in the parton fusion processes are extracted in cho2  by fitting CDF experimental data.

The PHENIX detector at RHIC has measured a few data points in the $`p_T`$ distribution of $`J/\psi `$ production at $`\sqrt{s}`$ = 200 GeV pp collisions. As these data are preliminary our comparison to this experimental data is not very conclusive. However, in the future the PHENIX detector at RHIC will measure high $`p_T`$ $`J/\psi `$ production up to $`p_T`$ $``$ 10 GeV with better statistics. The acceptance range of single electron for the PHENIX detector is $`0.35<\eta <0.35`$ and for single muon it is $`1.2<|\eta |<2.4`$. We compute the transverse momentum distribution of $`J/\psi `$ production by using both the color octet and the singlet model in the PHENIX detector acceptance range.

The $`22`$ parton fusion processes which contribute to the $`p_T`$ distribution of the $`J/\psi `$ and $`\chi _{cJ}`$ production are order $`\alpha _s^3`$. The differential cross section for the production of a quarkonium state ($`{}_{}{}^{2S+1}L_{J}^{}`$) with different quantum numbers (L,S,J) in the color singlet model is given by sing :

$`E{\displaystyle \frac{d^3\sigma }{d^3p}}(AB^{2S+1}L_J+X)`$ $`={\displaystyle \underset{i,j}{}}{\displaystyle 𝑑x_ix_if_{i/A}(x_i,Q^2)𝑑x_jx_jf_{j/B}(x_j,Q^2)}`$ (7)
$`\times {\displaystyle \frac{\widehat{s}}{\pi }}{\displaystyle \frac{d\widehat{\sigma }}{d\widehat{t}}}(ij(^{2S+1}L_J)k)\delta (\widehat{s}+\widehat{t}+\widehat{u}M^2).`$
Here $`x_i`$ and $`x_j`$ are the momentum fractions carried by partons $`i`$ and $`j`$ inside the hadrons A and B respectively. The $`f_{i/A}(x_i)`$ and $`f_{j/B}(x_j)`$ are the parton distribution function inside the hadron. $`\frac{d\widehat{\sigma }}{d\widehat{t}}(ij(^{2S+1}L_J)k)`$ is the partonic level differential cross section to form a quarkonium state ($`{}_{}{}^{2S+1}L_{J}^{}`$) which involves the radial wave function and its derivatives at the origin. Using the delta function integration in Eq. (7) one obtains:

$`{\displaystyle \frac{d\sigma }{dp_T}}(ABJ/\psi ,\chi _J+X)={\displaystyle \underset{i,j}{}}{\displaystyle 𝑑y𝑑x_ix_if_{i/A}(x_i,Q^2)x_j^{}f_{j/B}(x_j^{},Q^2)}`$
$`{\displaystyle \frac{2p_T}{x_i\frac{M_T}{\sqrt{s}}e^y}}{\displaystyle \frac{d\widehat{\sigma }}{d\widehat{t}}}(ij^{2S+1}L_Jk),`$ (8)
where

$$x_j^{}=\frac{1}{\sqrt{s}}\frac{x_j\sqrt{s}M_Te^yM^2}{x_j\sqrt{s}M_Te^y}.$$
(9)
In the above $`M`$ is the mass of the bound state quarkonium and $`M_T^2=\sqrt{p_T^2+M^2}`$. The partonic level differential cross section $`\frac{d\widehat{\sigma }}{d\widehat{t}}(ij(^{2S+1}L_J)k)`$ for the gluon fusion processes are obtained in gastmans  by using helicity decomposition method and the quark, gluon processes are derived in sing . The $`\frac{d\widehat{\sigma }}{d\widehat{t}}(ij(^{2S+1}L_J)k)`$ as given in sing  and gastmans  contain the non-relativistic wave function $`|R(0)|^2`$ (for direct $`J/\psi `$ process) and its derivatives $`|R^{}(0)|^2`$ (for $`\chi _J`$ processes) at the origin. For the non relativistic wave functions at the origin we take the Buchmuller-Type wave function with charm quark mass $`M_c`$=1.48 GeV. The numerical value is bm : $`|R(0)|^2`$=0.81 GeV3. For the derivative of the radial wave function at origin we use brat  $`\frac{9}{2\pi }\frac{|R^{}(0)|^2}{M_c^4}`$=15 MeV.

In the color octet model the differential cross section for heavy quarkonium production at order $`\alpha _s^3`$ is given by:

$`{\displaystyle \frac{d\sigma }{dp_T}}(AB\psi _Q(\chi _J)+X)={\displaystyle \underset{i,j}{}}{\displaystyle 𝑑y𝑑x_ix_if_{i/A}(x_i,Q^2)x_j^{}f_{j/B}(x_j^{},Q^2)}`$
$`\times {\displaystyle \frac{2p_T}{x_i\frac{M_T}{\sqrt{s}}e^y}}{\displaystyle \frac{d\widehat{\sigma }}{d\widehat{t}}}(ijC\overline{C}[^{2S+1}L_J^{(8)}]k\psi _Q(\chi _J)),`$ (10)
where

$`{\displaystyle \frac{d\widehat{\sigma }}{d\widehat{t}}}(ijC\overline{C}[^{2S+1}L_J^{(8)}]k\psi _Q(\chi _J))=`$
$`{\displaystyle \frac{1}{16\pi \widehat{s}^2}}\mathrm{\Sigma }|𝒜(ijC\overline{C}[^{2S+1}L_J^{(8)}]k)_{\mathrm{short}}|^2<0|𝒪_8^{\psi _Q(\chi _J)}(^{2S+1}L_J)|0>.`$ (11)
In this paper we include the contribution from the $`{}_{}{}^{3}P_{J}^{(8)}`$ and $`{}_{}{}^{1}S_{0}^{(8)}`$ processes which are present in the Fock state decomposition amplitude (5) in addition to $`{}_{}{}^{3}S_{1}^{(8)}`$ state of the color octet amplitude. Similar process is followed for the $`\chi _J`$ production by using (6) which takes into account the color octet $`{}_{}{}^{3}S_{1}^{(8)}`$ state. These $`\chi _J`$ formed in the color octet state again decay radiatively to give $`J/\psi `$. The total contribution to $`J/\psi `$ production in the color octet model is the sum of the direct $`J/\psi `$ and the contribution coming from $`\chi _J`$’s decay in the color octet model. In the above equation the short distance part of the partonic matrix element square in different quantum states for various partonic processes: $`\mathrm{\Sigma }|𝒜(ijC\overline{C}[^{2S+1}L_J^{(8)}]k)_{\mathrm{short}}|^2`$ are very lengthy and are given in cho2 . We use all the partonic processes in our calculation. We use the non-perturbative matrix elements: $`<0|𝒪_8^{\psi _Q(\chi _J)}(^{2S+1}L_J)|0>`$ from the tabulation given in cho2  which are extracted from the CDF data. In the extraction of the these matrix elements $`<0|𝒪_8^{\psi _Q}(^3P_0)|0>=M_Q^2<0|𝒪_8^{\psi _Q}(^1S_0)|0>`$ was set and then the values of the linear combination was extracted which is given by: $`\frac{<0|𝒪_8^{J/\psi }(^3P_0)|0>}{M_c^2}+\frac{<0|𝒪_8^{J/\psi }(^1S_0)|0>}{3}=(2.2\pm 0.5)\times 10^2`$ GeV3. For other matrix elements the values are given by: $`<0|𝒪_8^{J/\psi }(^3S_1)|0>=(6.6\pm 2.1)10^3`$ GeV3, $`<0|𝒪_8^{\chi _1}(^3S_1)|0>=(9.8\pm 1.3)10^3`$ GeV3. We use the central values of these matrix elements in our calculation.

We present our results for the rapidity range covered by the central electron arm and forward muon arm at the PHENIX detector at RHIC. The electron pseudo rapidity range: $`0.35<\eta <0.35`$ corresponds to the $`J/\psi `$ rapidity range: $`0.5<y<0.5`$. For the muon pseudo rapidity range: $`1.2<|\eta |<2.4`$ we present our results in the $`J/\psi `$ rapidity range $`1<y<2`$ which is well covered by the RHIC/PHENIX experiment. So the results we will present in this paper is in the $`J/\psi `$ rapidity range: $`0.5<y<0.5`$ and $`1<y<2`$. We will use GRV98 NLO grv98  parton distribution function throughout our calculation. The $`Q`$ scale in the structure function and the renormalization scale are chosen to be $`Q=M_T=\sqrt{M^2+p_T^2}`$, where $`M`$ and $`p_T`$ are the mass and the transverse momentum of the heavy quarkonium respectively. In extracting the NRQCD matrix elements, the MRSD0 parton distribution functions inside a nucleon was used in the color octet model calculation in cho2 . In this paper we are using, instead the GRV98 NLO parton distribution function for the parton distribution functions. We have checked that there is no significant difference in the results obtained by using MRSD0 and GRV98 in quantities calculated here.

In Fig. 1 we present the direct $`J/\mathrm{\Psi }`$ production cross section both from the color singlet and color octet models at $`\sqrt{s}`$= 500 and 200 GeV pp collisions in the $`J/\psi `$ rapidity range $`0.5<y<0.5`$. The solid line is the color octet contribution and dotted line is the color singlet contribution at $`\sqrt{s}`$=500 GeV. The dashed line is the color octet contribution and dot-dashed line is the color singlet contribution at $`\sqrt{s}`$=200 GeV. The direct $`J/\psi `$ production from color octet model at both the colliding energies are much higher than the color singlet counterparts. In what follows, we present the results of our computation for the $`p_T`$ range from 2 to 10 GeV. At smaller $`p_T`$, the color singlet and octet cross sections are corrupted by collinear divergences as discussed in cho1 ; cho2 . Since the intrinsic motion of the incident partons inside colliding protons renders the differential cross section uncertain for $`p_T`$ 2 GeV, as discussed in cho1 ; cho2  we have presented our results for $`p_T`$ 2 GeV in this paper. For low $`p_T`$ calculation one might follow the procedure where one might have to take intrinsic motion of the incident parton into account. This calculation is, however, beyond the scope of the present paper. In this paper we report the predictions for high $`p_T`$ $`J/\psi `$ production. High $`p_T`$ limit $``$ 10 GeV corresponds to the possible $`p_T`$ measurement range at PHENIX. In the very high $`p_T`$ limit ($`6GeV`$) the $`B`$ decay contribution $`BJ/\psi +X`$ and the gluon and charm quark fragmentation contribution might become important. However, we have not included these process here. In this paper we have computed only the prompt $`J/\psi `$ production from the parton fusion processes. These are the same contributions that were used in cho2  at Tevatron energies to extract NRQCD matrix elements.

In Fig. 2 we present the predictions for the $`J/\psi `$ production cross section coming from the radiative decay of the $`\chi _J`$’s both from the color octet and singlet model at $`\sqrt{s}`$= 500 and 200 GeV pp collisions in the $`J/\psi `$ rapidity range $`0.5<y<0.5`$. The solid line is the color octet contribution and dashed line is the color singlet contribution at $`\sqrt{s}`$=500 GeV. The dotted line is the color octet contribution and dot-dashed line is the color singlet contribution at $`\sqrt{s}`$=200 GeV. The $`\chi _{cJ}`$ contribution in each case in the figure is the sum of the contributions from $`\chi _0`$, $`\chi _1`$ and $`\chi _2`$. The $`\chi _{cJ}`$’s decay coming from the color singlet channel is slightly higher at small $`p_T`$ but very much lower at large $`p_T`$ in comparison to the color octet contribution. It can also be noted from Fig. 1 and Fig. 2 that the direct $`J/\psi `$ production in the color octet channel is larger than the contributions coming from all other processes. Thus at RHIC energies, we find that the dominant contribution to $`J/\psi `$ production comes from this direct color octet process.

In Fig.3 we present the results of our computation for the direct $`J/\psi `$ production for both the color singlet and octet models in the $`J/\psi `$ rapidity range $`1<y<2`$ at $`\sqrt{s}`$=500 and 200 GeV. The solid line is the contribution from the color octet process and dotted line is the contribution from the color singlet process at $`\sqrt{s}`$=500 GeV. The dashed line and dot-dashed lines are the similar curves at $`\sqrt{s}`$=200 GeV. In this rapidity regime we also find that the color octet contribution is much larger than the color singlet contribution. In Fig. 4 we present the $`J/\psi `$ production cross section coming from the $`\chi _J`$’s decay in the $`J/\psi `$ rapidity range $`1<y<2`$ at $`\sqrt{s}`$= 500 and 200 GeV pp collisions at RHIC. The solid line is the contribution from the color octet process and dashed line is the contribution from the color singlet process at $`\sqrt{s}`$=500 GeV. The dotted line and dot-dashed lines are the similar curves at $`\sqrt{s}`$=200 GeV. The total $`\chi _{cJ}`$ contributions in the figure includes the contribution from $`\chi _0`$, $`\chi _1`$ and $`\chi _2`$. As $`\psi ^{}`$ contribution to $`J/\psi `$ is small at this energy range we have not considered the $`\psi ^{}`$ contribution in our calculation. As can be seen again the dominant contribution to the $`J/\psi `$ production comes mainly from the direct $`J/\psi `$ production in the color octet channel.

The total $`J/\psi `$ $`p_T`$ distribution coming from combined color singlet model and color octet model contributions in all channels is plotted in Fig. 5 in the $`J/\psi `$ rapidity range $`0.5<y<0.5`$ at $`\sqrt{s}`$=500 and 200 GeV pp collisions. The solid line is the $`p_T`$ distribution of the total (octet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay) plus singlet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay)) $`J/\psi `$ production cross section at $`\sqrt{s}`$=500 GeV. The dashed line is the contribution from the color octet model and dot-dashed line is the contribution from the singlet model at the same energy. As pointed out in the above the total cross section is the sum of the cross section from direct $`J/\psi `$ and from the radiative decays of the $`\chi _J`$’s both in color octet and singlet channels. It can be seen that the total cross section is almost same as the cross section obtained from the color octet channels. The color singlet channel contribution is very small as can be seen from the dot-dashed line.

Similar results are obtained for $`\sqrt{s}`$=200 GeV pp collisions in the same figure in the $`J/\psi `$ rapidity range $`0.5<y<0.5`$. The dotted line is the total (octet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay) plus singlet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay)) $`J/\psi `$ production cross section at $`\sqrt{s}`$=200 GeV. Long-dashed line is the octet contribution and thin-solid line is the singlet contribution. It can be seen from the figure that the total $`J/\psi `$ production cross section is almost same as that of the color octet model contribution and the singlet contribution is much smaller. A measurement of the $`J/\psi `$ differential cross section at RHIC in the future (both at $`\sqrt{s}`$ = 500 and 200 GeV pp collisions) will be able to tell us whether the color octet model reproduces the RHIC data.

As PHENIX has measured few data points for $`p_T`$ distribution of $`J/\psi `$ in the $`\mu ^+\mu ^{}`$ rapidity range in pp collisions at $`\sqrt{s}`$= 200 GeV, we will compare our results with these available experimental data. As stated earlier we will compare the available data in the $`J/\psi `$ rapidity range $`1<y<2`$ as this range is well covered by the RHIC/PHENIX experiment in the muon channel. In Fig. 6 we present the $`p_T`$ distribution of total $`J/\psi `$ production cross section from color singlet and octet model in the $`J/\psi `$ rapidity range $`1<y<2`$ at $`\sqrt{s}`$=500 GeV. The solid line is the $`p_T`$ distribution of the total (octet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay) plus singlet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay)) $`J/\psi `$ production cross section at $`\sqrt{s}`$=500 GeV. The dashed line is the contribution from the color octet model and dot-dashed line is the contribution from the singlet model at the same energy. The $`p_T`$ distribution of total $`J/\psi `$ production cross section from color octet and singlet model in the $`J/\psi `$ rapidity range $`1<y<2`$ at $`\sqrt{s}`$=200 GeV pp collisions is plotted in Fig. 6. The dotted line is the total (octet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay) plus singlet (direct $`J/\psi `$ \+ $`\chi _J`$’s decay)) $`J/\psi `$ production cross section at $`\sqrt{s}`$=200 GeV. Long-dashed line is the octet contribution and thin-solid line is the singlet contribution.

Let us compare our results with the preliminary experimental data obtained at the PHENIX detector. As mentioned earlier, at small $`p_T`$ the color singlet and octet cross sections are corrupted by collinear divergences as discussed in cho2 . Since the intrinsic motion of the incident partons inside colliding protons renders the differential cross section uncertain for $`p_T`$ 2 GeV, as discussed in cho1 ; cho2  we will compare our results at $`p_T`$ 2 GeV with the PHENIX data qm02  in this paper. We compare our calculation with the only available data point ($``$ 2 GeV) at $`p_T`$ 2.5 GeV (see qm02 ). The data point is presented as ”+” in Fig. 6. It can be seen that the present experimental data fits the color octet model predictions of $`J/\psi `$ production and the color singlet model prediction is well below the data point. Hence the color singlet model is not the relevant mechanism at RHIC energy. As for the color octet model, the single data point matches with the color octet contribution, but many more data with high precision measurements are necessary to confirm the validity of the color octet model at RHIC energy. The next run results from PHENIX experiment will be able to shed light on this.

In summary, we have computed the $`J/\psi `$ production differential cross section in pp collisions at RHIC at $`\sqrt{s}`$ = 500 and 200 GeV by using both the color octet and singlet models in the framework of non-relativistic QCD. The $`J/\psi `$ production cross section we compute here includes the direct $`J/\psi `$ production from the partonic fusion processes and the $`J/\psi `$ coming from the radiative decays of $`\chi _J`$’s both in the color octet and singlet channel. The high $`p_T`$ $`J/\psi `$ production cross section is computed within the $`J/\psi `$ rapidity range: $`0.5<y<0.5`$ and $`1<|\eta |<2`$, in the the central electron and forward muon arms respectively. It is found that the color octet contribution to $`J/\psi `$ production is dominant at RHIC energy in comparison to the color singlet contributions. We have compared our results with the recent data obtained by PHENIX detector for $`p_T`$ 2 GeV at $`\sqrt{s}`$ =200 GeV pp collisions. While the color singlet model fails to explain the data completely the color octet model seems to reproduce the data. As the present high $`p_T`$ data of $`J/\psi `$ production at PHENIX detector at RHIC is very restricted, measurement of high $`p_T`$ $`J/\psi `$ production at RHIC in the next run with better statistics will be able to tell us the degree of validity of the color octet model of $`J/\psi `$ production at RHIC energies. It is necessary to know the exact mechanism of $`J/\psi `$ production in pp collisions at RHIC in order to make a prediction of $`J/\psi `$ suppression as a signature of the quark-gluon plasma in Au-Au collisions. Knowledge of the $`J/\psi `$ production mechanism is also crucial for isolating the polarized spin structure functions of the proton from the data.

###### Acknowledgements.

This research is supported by the Department of Energy, under contract W-7405-ENG-36. We thank Rajiv Gavai, Mike Leitch, Pat McGaughey, Emil Mottola and Johann Rafelski for useful discussions.