%\documentstyle[psfig,12pt]{llncs}
\documentclass[12pt]{article}
\usepackage{psfig} 

% variant of \vec command from llncs,
% but more my style ...
\def\VEC#1{{\mathchoice%
{\mbox{\boldmath$\displaystyle #1$}}%
{\mbox{\boldmath$\textstyle #1$}}%
{\mbox{\boldmath$\scriptstyle #1$}}%
{\mbox{\boldmath$\scriptscriptstyle #1$}}}}

% miscellaneous math formatting and symbols
\newcommand{\OL}{\overline}
\newcommand{\UL}{\underline}
\newcommand{\FUN}[1]{{\rm #1}}
\newcommand{\MAT}[1]{{\bf #1}}
\newcommand{\LAB}[1]{{\sf #1}}
\newcommand{\EMB}[1]{{\VEC #1}}
\newcommand{\FRAC}[2]{\mbox{$\frac{#1}{#2}$}}
\newcommand{\SQRT}[1]{\mbox{$\sqrt{#1}$}}
\newcommand{\KET}[1]{|#1\rangle}
\newcommand{\BRA}[1]{\langle#1|}
\newcommand{\AVG}[1]{\langle#1\rangle}
\newcommand{\HALF}{{\FRAC{1}{2}}}
\newcommand{\RHO}{{\EMB\rho}}
\newcommand{\VRHO}{{\EMB\varrho}}
\newcommand{\EP}{\VEC{E}_{+}}
\newcommand{\EM}{\VEC{E}_{-}}
\newcommand{\FP}{\VEC{F}_{+}}
\newcommand{\FM}{\VEC{F}_{-}}
\newcommand{\GP}{\VEC{G}_{+}}
\newcommand{\GM}{\VEC{G}_{-}}
\newcommand{\EPM}{\VEC{E}_{\pm}}
\newcommand{\EMP}{\VEC{E}_{\mp}}
\newcommand{\FPM}{\VEC{F}_{\pm}}
\newcommand{\FMP}{\VEC{F}_{\mp}}
\newcommand{\GPM}{\VEC{G}_{\pm}}
\newcommand{\GMP}{\VEC{G}_{\mp}}
\newcommand{\IP}{\VEC{I}_{+}}
\newcommand{\IM}{\VEC{I}_{-}}
\newcommand{\IPM}{\VEC{I}_{\pm}}
\newcommand{\IMP}{\VEC{I}_{\mp}}
\newcommand{\IX}{\VEC{I}_\LAB{x}}
\newcommand{\IY}{\VEC{I}_\LAB{y}}
\newcommand{\IZ}{\VEC{I}_\LAB{z}}

%\addtolength{\footnotesep}{6pt}
\renewcommand{\footnoterule}{\vspace{3pt}%
\kern-3pt\hrule width 2 true cm\kern 2.6pt}
\renewcommand{\arraystretch}{1.2}
\renewcommand{\baselinestretch}{1.2}
\newcommand{\skiplinehalf}{\vspace{6pt}}

\pagestyle{plain}
\flushbottom	% sloppy spacing vs. my time
\setlength{\parskip}{0pt plus 0.5\baselineskip}
\setlength{\topsep}{0.125\baselineskip plus 0.125\baselineskip}
\setlength{\itemsep}{0.125\baselineskip plus 0.125\baselineskip}

\begin{document}% comes before \maketitle in this style

\title{Principles and Demonstrations of Quantum
Information Processing by NMR Spectroscopy\thanks{
Portions of this survey were presented at the AeroSense
Workshop on Photonic Quantum Computing II, held in Orlando,
Florida on April 16, 1998, at the Dagstuhl Seminar
on Quantum Algorithms, held in Schloss Dagstuhl,
Germany on May 10 -- 15, 1998, and at the Workshop
on Quantum Information, Decoherence and Chaos, held
on Heron Island, Australia September 21 -- 25, 1998;
this paper is an updated and extended version of one
published in the proceedings of the AeroSense meeting,
available as vol.\ 3385 from the International Society for
Optical Engineering, 1000 20th St., Bellingham, WA 98225, USA.}}

%\author{T.~F.~Havel\inst{1}, S.~S.~Somaroo\inst{1},
%C.-H.~Tseng\inst{2} and D.~G.~Cory\inst{3}}

%\institute{BCMP, Harvard Medical School, Boston, MA 02115
%\and Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138
%\and Nuclear Engineering, MIT, Cambridge, MA 02139}

\maketitle

\begin{abstract}
This paper surveys our recent research on quantum information
processing by nuclear magnetic resonance (NMR) spectroscopy.
We begin with a brief introduction to the product operator
formalism, on which the theory of NMR spectroscopy is based,
and use it throughout the rest of the paper to show
how it provides an concise framework within which
to analyze quantum computations and decoherence.
The implementation of quantum algorithms by NMR depends
upon the availability of special kinds of mixed states,
called pseudo-pure states, and we consider a number
of different methods for preparing pseudo-pure states,
along with what is known about how they scale with the number of spins.
The quantum-mechanical nature of processes involving such
macroscopic pseudo-pure states also is a matter of debate,
and we attempt to make this debate more concrete by presenting
the results of NMR experiments which validate Hardy's paradox,
subject to certain assumptions that we explicitly state.
Finally, a detailed product operator description is
given of recent NMR experiments which demonstrate the
principles behind a three-bit quantum error correcting code.
\end{abstract}

\eject
\section{INTRODUCTION}
An {\em ensemble quantum computer\/} consists of a massive
collection of independent and identical quantum computers,
which all operate on the same coherent superposition in lock-step,
and for which one can determine the sum of their one-bit
observables over the ensemble in a single measurement.
This enables it to use classical ensemble parallelism to
estimate the {\em expectation values\/} of the observables,
rather than random eigenvalues thereof \cite{CorFahHav:97}.
Although we do not expect any asymptotic performance gains
with an ensemble of fixed size, an ensemble quantum computer
retains all the essential features of a quantum computer,
including entanglement and interference,
and hence can solve the same problems in polynomial time.
The further ability to estimate expectation values can
potentially reduce the amount of time that would be
needed on an equivalent quantum computer by a very large,
albeit constant factor, for example when using Grover's
algorithm to solve search problems \cite{BoBrHoTa:98,Grover:97a}.

The main reason an ensemble quantum computer is interesting,
however, is that it can be implemented with surprising ease,
at least on a limited scale, by NMR spectroscopy
on ordinary liquids at room temperature and pressure.
This implementation takes advantage of several
very favorable features of NMR spectroscopy on
the spin $= 1/2$ nuclei of molecules in a liquid,
in particular their seconds-long decoherence times.
The implementation also depends on the availability
of special kinds of mixed states, whose behavior
under coherent Hamiltonians is identical to that of
true pure states \cite{CorFahHav:97,GershChuan:97}.
The general methods that are presently available
for preparing such {\em pseudo-pure\/} states from
the equilibrium state of the spins in a magnetic
field result in a rapid loss of observable
magnetization with increasing number of spins.
This in turn restricts such implementations to
ca.\ 8 -- 10 spins in the foreseeable future.
Although this is certainly a severe restriction,
NMR has now provided direct experimental demonstrations
of many of the basic elements of quantum information
processing \cite{ChGeKuLe:98,CorFahHav:96,CorPriHav:98},
including simple quantum algorithms
\cite{ChVaZhLeLl:98,ChuGerKub:98,JonMosHan:98,JonesMosca:98},
teleportation \cite{NieKniLaf:98},
and error correcting codes \cite{CMPKLZHS:98}.

This paper will survey our recent research on
quantum information processing by NMR spectroscopy.
We begin with a brief description of the mathematical
framework within which NMR experiments are widely analyzed,
known as the {\em product operator formalism\/}, and show
how it applies to quantum computing (cf.~\cite{SomCorHav:98}).
Next, we give an overview of the basic ideas behind
ensemble quantum computing by NMR spectroscopy,
with emphasis on pseudo-pure state preparation and scaling.
We then consider some ways in which quantum correlations
can be manifest even in weakly polarized spin ensembles,
and illustrate this with the results of NMR experiments which
(with certain assumptions) validate Hardy's paradox \cite{Hardy:92}.
Finally, the utility of NMR and its associated
product operator formalism as a means of studying
decoherence will be demonstrated by an analysis
of our recent experiments with a three-bit quantum
error correcting code \cite{CMPKLZHS:98}.
The reader is assumed throughout to be familiar with
the basic notions of quantum information processing,
as presented in e.g.\ Refs.\ \cite{Preskill:98,Steane:98,WilliClear:98}.
A more introductory account of our work on
ensemble quantum computing by NMR spectroscopy
may be found in Ref.\ \cite{CorPriHav:98}.

\section{THE PRODUCT OPERATOR FORMALISM}
The basic observables in the NMR spectroscopy are
the $\LAB{x}$- and $\LAB{y}$-components of the spin
angular momentum (traditionally in units of $\hbar$).
These give rise to absorptive and dispersive peaks,
respectively, centered on the resonance frequencies of the
spin in the applied static magnetic field \cite{Slichter:90}.
The operators for the angular momentum components
of a spin are denoted by $\IX$, $\IY$ and $\IZ$,
which for the spin $= 1/2$ nuclei dealt with in this paper
are typically represented by one half the Pauli matrices.
Although products of the angular momentum components
of different spins cannot be observed directly,
they evolve naturally into single components,
and hence it is also natural to
define operators for such products.
This can be done using the {\em mixed product formula\/}
between the operator composition product and tensor products:
\begin{equation}
(\VEC{A} \otimes \VEC{B})(\VEC{C} \otimes \VEC{D})
~=~ (\VEC{AC}) \otimes (\VEC{BD})
\end{equation}

Thus we have, for example, the following operator
for the product of the $\LAB{x}$ component of the
first spin with the $\LAB{y}$ component of a second:
\begin{equation} \begin{array}{rcl}
\BRA{\phi^1}\IX^1\KET{\phi^1} \BRA{\psi^2}\IY^2\KET{\psi^2}
&~=~& \BRA{\phi}\IX\KET{\phi} \otimes \BRA{\psi}\IY\KET{\psi}
\\ &~=~&
(\BRA{\phi}\otimes\BRA{\psi}) (\IX \otimes \IY)
(\KET{\phi}\otimes\KET{\psi})
\\ &~=~&
\BRA{\phi^1\psi^2} (\IX \otimes \VEC{1})
(\VEC{1} \otimes \IY) \KET{\phi^1\psi^2}
\\ &~=~&
\BRA{\phi^1\psi^2} \IX^1 \IY^2 \KET{\phi^1\psi^2}
\end{array} \end{equation}
Operators like $\IX^1\IY^2$ are called {\em product operators\/}
\cite{BoulaRance:94a,ErnBodWok:87,SomCorHav:98,SoEiLeBoEr:83,vdVenHilbe:83}.
The mixed product formula,
together with the multiplication rules inherited from
the Pauli matrix algebra, gives us a well-defined
multiplication between any pair of product operators.
Henceforth, we shall identify the identity
operator $\VEC{1}$ with the scalar identity $1$.

Product operators enable us to represent not only
the observables of NMR, but also the spin states and
their transformations within a single algebraic system.
To explain this, we first note that the state of the spins
in an NMR sample at room temperature is necessarily highly mixed,
and so must be described by a density operator \cite{Blum:81}.
We shall also use density operators to describe pure states.
For example, the density operator of two spins
antiparallel along the $\LAB{z}$-axis is
\begin{equation}
\KET{01}\BRA{01} ~=~
(\KET{0}\BRA{0})\otimes(\KET{1}\BRA{1})
~=~ \HALF(1 + 2\IZ^1) \HALF(1 - 2\IZ^2)
~\equiv~ \EP^1 \EM^2 ~.
\end{equation}
It is easily verified that the operators
$\EPM^k = \HALF(1 \pm 2\IZ^k)$ of the
individual spins are {\em idempotent\/},
meaning they are equal to their squares,
and that they are mutually {\em annihilating\/},
i.e.\ $\EP^k\EM^k = 0 = \EM^k\EP^k$.
A general pure state corresponds to a
{\em primitive\/} idempotent, meaning it cannot
be expressed as a sum of two nonzero idempotents,
and the pure state is {\em uncorrelated\/} if it can
be factorized into one-spin idempotents (as above).
A general mixed state is given by a sum of
(not necessarily factorizable) primitive idempotents
with nonnegative coefficients $p_i$ summing to one,
which are the probabilities of the corresponding
pure states in the thermodynamic ensemble:
\begin{equation}
\RHO ~\equiv~ \overline{\KET{\psi}\BRA{\psi}}
~=~ {\sum}_i \, p_i \, \KET{\psi_i}\BRA{\psi_i}
\end{equation}

This representation of quantum states
by density operators makes it clear that,
if $\VEC{U}$ is a unitary transformation
which acts on state vectors by left-multiplication
$\KET{\psi} \rightarrow \VEC{U}\KET{\psi}$,
then it acts on the corresponding
density operator by conjugation, i.e.
\begin{equation}
\RHO ~=~ \overline{\KET{\psi}\BRA{\psi}}
\quad\longrightarrow\quad
\overline{(\VEC{U} \KET{\psi}) {({\VEC{U} \KET{\psi}})}^{\sim}}
~=~ \VEC{U} \overline{\KET{\psi}\BRA{\psi}} \tilde\VEC{U}
~=~ \VEC U \RHO\, \tilde\VEC{U}
\end{equation}
(where $\tilde\VEC{U} \equiv \VEC{U}^{\sim}$
denotes the Hermitian conjugate).
It is also easily seen that the ensemble-average
expectation value of any observable $\VEC A$ is given
by the trace of its product with the density operator:
\begin{equation}
\overline{\BRA{\psi}\VEC A\KET{\psi}} ~=~
\FUN{Tr}(\overline{\BRA{\psi}\VEC A\KET{\psi}})
~=~ \FUN{Tr}(\VEC A\overline{\KET{\psi}\BRA{\psi}})
~\equiv~ \FUN{Tr}(\VEC A\RHO)
\end{equation}
This in turn is just $2^N$ times the (real) {\em scalar part\/}
$\langle\VEC A\RHO\rangle \equiv 2^{-N}\,\FUN{Tr}(\VEC A\RHO)$
of the product.\footnote{
The factor of $2^{-N}$ in the scalar
part is not as unnatural as it may seem,
since as we shall see it also appears
in the partition function. }
For example, using the relations $2\AVG{\EPM^k} = 1$
and $2 \EPM^k \IZ^k = \pm \EPM^k$, we can show that the
$\LAB{z}$-component of the total angular momentum depends
on the polarizations $p$ and $q$ of two uncorrelated spins as
\begin{equation} \begin{array}{rcl}
&& 4 \left\langle ((1-p)\EP^1+p\,\EM^1)
((1-q)\EP^2+q\,\EM^2) (\IZ^1+\IZ^2) \right\rangle
\\ &~=~&
4 \left\langle (1-p)(1-q)\EP^1\EP^2 - pq\,\EM^1\EM^2 \right\rangle
\\ &~=~&
(1-p)(1-q) - pq ~=~ 1 - p - q ~.
\end{array} \end{equation}

Turning now to the transformations,
we first note that NMR spectroscopists
typically use an interaction representation,
wherein all the transformations and states
are referred to a right-handed coordinate frame
which co-rotates at the receiver frequency about
the axis of precession $\LAB{z}$ \cite{Slichter:90}.
The unitary transformation that corresponds to a
right-hand rotation of the $k$-th spin by an angle
$\theta$ about the $\LAB{x}$ axis in this frame
is given by
\begin{equation} \begin{array}{rcl}
e^{-\imath\theta\IX^k} &~=~& 1 -
\imath \left(\FRAC{\theta}{2}\right) 2\IX^k -
\HALF \left(\FRAC{\theta}{2}\right)^2 +
\FRAC{\imath}{6} \left(\FRAC{\theta}{2}\right)^3
\,2 \IX^k + \cdots \\ &~=~&
\cos\left(\FRAC{\theta}{2}\right) -
\imath \sin\left(\FRAC{\theta}{2}\right) 2\IX^1 ~,
\end{array} \end{equation}
where we have used $(2\IX^k)(2\IX^k) = 1$.
In NMR, this operation is implemented by
an RF (radio-frequency) pulse whose frequency
range spans the resonances of the $k$-th spin.
We shall also encounter {\em correlated rotations\/},
in particular that induced by the bilinear
interaction known in NMR as {\em scalar coupling\/},
\begin{equation}
e^{-\imath \,t\, \VEC{H}_\LAB{J}} ~=~
e^{-\imath \,t\, 2\pi J^{12} \IZ^1\IZ^2} ~=~
\cos\left(\FRAC{\pi}{2}J^{12}t\right) - \imath
\sin\left(\FRAC{\pi}{2}J^{12}t\right) 4\IZ^1\IZ^2 ~,
\end{equation}
where $\VEC{H}_\LAB{J}$ is the (weak) coupling
Hamiltonian and $J^{12}$ the coupling constant in Hz.
As we shall see momentarily, this interaction permits
us to implement general quantum logic gates by NMR.

In quantum computing, the NOT operation on the $k$-th spin
simply rotates it by $\pi$; according to the above formula:
\begin{equation} \begin{array}{rcl}
e^{-\imath\pi\IX^k} \EP^k e^{\imath\pi\IX^k}
&~=~& (-2\imath\IX^k) \EP^k (2\imath\IX^k)
~=~ \HALF (1 + 8 \IX^k \IZ^k \IX^k) \\
&~=~& \HALF (1 - 8 \IX^k \IX^k \IZ^k)
~=~ \HALF (1 - 2\IZ^k) ~=~ \EM^k 
\end{array} \end{equation}
The c-NOT (controlled-NOT) operation is a $\pi$
rotation of e.g.\ the first spin {\em conditional\/}
on the polarization of a second.
Using the relation $\EPM^k \EMP^k = 0$,
we can easily show that
\begin{equation}
(-2\imath\IX^1 \EM^2 + \EP^2) (\VEC{E}_\epsilon^1 \EPM^2)
(2\imath\IX^1 \EM^2 + \EP^2) ~=~ \VEC{E}_{\pm\epsilon}^1 \EPM^2
\end{equation}
($\epsilon \in \{\pm{}\}$).
The phase factor $\imath$ multiplying $\IX^1$
complicates the action of the c-NOT on a superposition,
but can be eliminated by a conditional phase shift.
Using $\EM^k + \EP^k = 1$, this phase-corrected
c-NOT gate is seen to be $\VEC{S}^{1|2} \equiv$
\begin{equation} \begin{array}{rcl}
2\IX^1 \EM^2 + \EP^2 &~=~&
(-\imath \EM^2 + \EP^2) (2\imath\IX^1 \EM^2 + \EP^2)
\\ &~=~&
\left( 1 + (e^{-\imath\frac\pi2} - 1) \EM^2 \rule{0pt}{10pt} \right)
\left( 1 + (e^{\imath\pi\IX^1} - 1) \EM^2 \right)
\\ &~=~&
e^{-\imath\frac{\pi}{2}\EM^2} e^{\imath\pi\IX^1\EM^2}
~=~ e^{-\imath\frac{\pi}{2}(1 - 2\IX^1)\EM^2} ~,
\end{array} \end{equation}
and hence the idempotents $\EPM^k$ also give
us an algebraic description of the c-NOT gate,
in addition to the density operators of pure states.
It is well-known that single spin rotations,
together with the c-NOT, are sufficient to
implement any quantum logic gate \cite{BBCDMSSSW:95}.

It is also well-known how these basic operations
can be implemented by NMR \cite{CorPriHav:98}.
The individual spins can be rotated by selective
radio-frequency pulses, as mentioned above.
Conditional rotations, and in particular the c-NOT,
can be implemented by combining selective
rotations with scalar coupling.
Noting that
\begin{equation}
e^{\imath \frac{\pi}{2} \IX^1}
e^{-\imath \frac{\pi}{2} \IY^1} ~=~
e^{\imath \frac{\pi}{2} \IZ^1}
e^{\imath \frac{\pi}{2} \IX^1} ~,
\end{equation}
we obtain:
\begin{equation} \begin{array}{rcl}
\VEC S^{1|2} &~=~&
e^{-\imath\frac{\pi}{2}\EM^2}
e^{\imath\pi\IX^1\EM^2}
\\ &~=~&
e^{-\imath\frac{\pi}{4}(1-2\IZ^2)}
e^{\imath\frac{\pi}{2}\IX^1(1-2\IZ^2)}
\\ &~=~&
e^{-\imath\frac{\pi}{4}} e^{\imath\frac{\pi}{2}\IZ^2}
e^{\imath\frac{\pi}{2}\IX^1} e^{-\imath\pi\IX^1\IZ^2}
\\ &~=~&
e^{-\imath\frac{\pi}{4}} e^{\imath\frac{\pi}{2}\IZ^2}
e^{\imath\frac{\pi}{2}\IX^1} e^{-\imath\frac{\pi}{2}\IY^1}
e^{-\imath\pi\IZ^1\IZ^2} e^{\imath\frac{\pi}{2}\IY^1}
\\ &~=~&
e^{-\imath\frac{\pi}{4}} e^{\imath\frac{\pi}{2}\IZ^2}
e^{\imath\frac{\pi}{2}\IZ^1} e^{\imath\frac{\pi}{2}\IX^1}
e^{-\imath\pi\IZ^1\IZ^2} e^{\imath\frac{\pi}{2}\IY^1}
\end{array} \end{equation}
Since the overall phase of the transformation
has no effect on the density operator, it follows that
the c-NOT $\VEC{S}^{1|2}$ can be implemented in NMR by
applying the following sequence of effective Hamiltonians:
\begin{equation}
[-\FRAC{\pi}{2}\IY^1] \rightarrow [\pi\IZ^1\IZ^2] \rightarrow
[-\FRAC{\pi}{2}\IX^1] \rightarrow [-\FRAC{\pi}{2}(\IZ^1+\IZ^2)]
\end{equation}
In practice, the effective Hamiltonian $[\pi\IZ^1\IZ^2]$
is obtained by applying a $\pi$-pulse to both spins in the
middle and at the end of a $1/(2J^{12})$ evolution period,
to refocus the Zeeman evolution \cite{CorPriHav:98}.
The $[-\FRAC{\pi}{2}\IY^1]$ and $[-\FRAC{\pi}{2}\IX^1]$
Hamiltonians are implemented by RF pulses as above,
while the $[-\FRAC{\pi}{2}(\IZ^1+\IZ^2)]$ transformation
can most easily be implemented by letting one spin
evolve while refocusing the other, then vice versa,
and finally realigning the transmitter's phase with the spins'.

In closing this section, we note that the
algebra of a single spin is isomorphic to the
{\em Clifford algebra\/} of space \cite{SomCorHav:98}.
This is naturally contained within the Dirac
matrix algebra, which in turn is isomorphic
to the Clifford algebra of space-time.
A relativistic multiparticle version of
the latter is given by the Clifford algebra
of a direct sum of multiple copies of space-time,
and contains an algebra isomorphic to the product
operator algebra \cite{DorLasGul:93,DoLaGuSoCh:96}.
Thus, the product operator algebra is not just
a specialized tool for NMR spectroscopists,
but plays a fundamental role in physics.

\section{PSEUDO-PURE STATE PREPARATION AND SCALING}
In liquids, the dipole-dipole interaction
among the spins is averaged to zero by the
rapid rotational diffusion of the molecules.
This effectively makes each molecule
into an independent quantum computer.
It also enables the density operator to
be factorized into a product of density
operators for the individual molecules,
\begin{equation}
\VRHO ~=~ \VRHO^1 \cdots \VRHO^M ~=~
\RHO^1 \otimes\cdots\otimes \RHO^M ~,
\end{equation}
where $\VRHO^k \equiv \MAT 1 \otimes \cdots
\otimes \RHO^k \otimes \cdots \otimes \MAT 1$.
In a pure liquid (or if we are looking
at just one component of a solution),
all the molecules are equivalent so that
all these density operators are the same.
It follows that we can work with a {\em reduced\/}
density operator $\EMB\rho ~\equiv~ \RHO^1 ~=\cdots=~
\RHO^M$ of dimension $2^N$, rather than $2^{MN}$.

The equilibrium density operator $\RHO_\LAB{eq}$ of
an $N$-spin system in a magnetic field is given by
the Boltzman operator determined by its Hamiltonian,
$\exp(-\VEC{H}/k_\LAB{B}T)$, divided by the partition
function $Z_\LAB{eq} = \FUN{Tr}(\exp(-\VEC{H}/k_\LAB{B}T))$.
At ``high temperatures'' (which, given the gyromagnetic
ratios of nuclear spins, means more than a few $mK$),
$\VEC{H}$ is very small compared to $k_\LAB{B}T$,
so that a linear approximation is quite accurate:
\begin{equation}
\RHO_\LAB{eq} ~\approx~ (1 - \VEC{H} / k_\LAB{B}T) / \FUN{Tr}(
1 - \VEC{H} / k_\LAB{B}T) ~=~ (1 - \VEC{H} / k_\LAB{B}T) / 2^N
\end{equation}
The identity component $1 / 2^N$
gives rise to no net magnetization,
since it represents a state in which
each spin is up as often as it is down.
The Hamiltonian $\VEC{H}$, on the other hand,
is well-approximated by the dominant Zeeman precession
term $\VEC{H}_\LAB{Z}$, which can be written as
\begin{equation}
\VEC{H}_\LAB{Z} ~\equiv~ -\hbar\, B_\LAB{z}
(\gamma^1 \IZ^1 +\cdots+ \gamma^N \IZ^N)
\end{equation}
where $\gamma^k$ is the gyromagnetic ratio of the $k$-th spin,
and $B_\LAB{z}$ is the $\LAB{z}$-component of the magnetic field.
In practice, NMR spectroscopists
usually drop the unobservable identity
component from the density operator.
Moreover, in homonuclear systems
they also leave out both the constant
factor $\hbar B_\LAB{z} \gamma^k$ as
well as the partition function $2^{-N}$ ---
though we shall carefully avoid doing the latter!
In these terms, the equilibrium density
operator of a two-spin system,
and its matrix representation in
the standard $\IZ$-eigenbasis, is
\renewcommand{\arraystretch}{1.0}
\begin{equation} \begin{array}{rcl}
&& \hat{\RHO}_\LAB{eq} ~=~ \FRAC{1}{4} (\IZ^1 + \IZ^2)
~=~ \FRAC{1}{4} (\KET{00}\BRA{00} - \KET{11}\BRA{11})
\vspace*{3pt} \\ &~\leftrightarrow~&
\FRAC{1}{4}\, \MAT{Diag}( 1, 0, 0, -1 )
~\equiv~ {\displaystyle\frac{1}{4}}
\left[ \begin{array}{rrrr}
1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\ 0 & ~~0 & ~~0 & -1
\end{array} \right] ~,
\end{array} \end{equation}
\renewcommand{\arraystretch}{1.2}
where the ``hat'' on $\RHO_\LAB{eq}$
signifies its traceless part.

In contrast, the density operator of two
spins in their pseudo-pure ground state is
\begin{equation} \begin{array}{rcl}
\hat{\RHO}_{00} &~\equiv~& \pm\FRAC{1}{6}
(\IZ^1 + \IZ^2 + 2 \IZ^1 \IZ^2) ~=~
\pm\FRAC{1}{3} \left( \EP^1 \EP^2 - \FRAC{1}{4} \right)
\\ &~=~&
\pm\FRAC{1}{3} \left( \KET{00}\BRA{00} - \FRAC{1}{4} \right)
~\leftrightarrow~ \pm\FRAC{1}{12}\, \MAT{Diag}( 3, -1, -1, -1 ) ~.
\end{array} \end{equation}
The overall sign depends on whether we have a
population excess or deficit in the ground state;
for consistency, we shall generally assume the former.
Observe that a unitary transformation of the density
operator induces a transformation of the corresponding
state vector just as it does for true pure states, since
\begin{equation}
\VEC{U} \hat{\RHO}_{00}\, \tilde\VEC{U} ~=~
\FRAC{1}{3} \left( \left( \VEC{U} \KET{00} \right)
\left( \VEC{U} \KET{00} \right)^{\sim} - \FRAC{1}{4} \right) ~.
\end{equation}
Similarly, because the NMR observables $\VEC A = \IX^k, \IY^k$
are traceless, the ensemble-average expectation value relative
to a pseudo-pure density operator yields the ordinary
expectation value versus the corresponding state vector:
\begin{equation}
\FUN{Tr}( \VEC{A} \hat{\RHO}_{00} ) ~=~ \FRAC{1}{3} \left( \FUN{Tr}
( \VEC{A} \KET{00}\BRA{00} ) - \FRAC{1}{4} \FUN{Tr}( \VEC{A} ) \right)
~=~ \FRAC{1}{3} \BRA{00} \VEC{A} \KET{00}
\end{equation}
The general form of a pseudo-pure density operator is
\begin{equation} \label{eq:ppform}
\hat{\RHO}_\LAB{\psi} ~=~ \FRAC{N/2}{\rule{0pt}{6pt}2^N-1}
\left( \KET{\psi}\BRA{\psi} - 2^{-N} \right) ~,
\end{equation}
where $\KET{\psi}$ is a normalized $N$-spin state vector,
and the prefactor has been chosen so as to keep
the maximum eigenvalue $\| \hat{\RHO}_{\psi} \|_2$ equal
to that of the $N$-spin equilibrium density operator.

Even though we have defined them
to have the same maximum eigenvalue,
for $N > 1$ the remaining eigenvalues of
$\hat{\RHO}_\LAB{eq}$ and $\hat{\RHO}_{\psi}$ differ,
and hence there is no unitary
transformation taking one to the other.
There are nevertheless a number
of nonunitary processes by which
one can prepare pseudo-pure states.
The most direct is to generate a spatially
varying distribution of states across the sample,
such that the ensemble average is pseudo-pure.
This can be done by using a {\em field gradient\/}
along the $\LAB{z}$-axis to create a
position-dependent phase whose average is zero,
thereby in effect setting the transverse ($\LAB{xy}$)
components of the density operator to zero.\footnote{
In the homonuclear case, the zero-quantum coherences
are not rapidly dephased by a $\LAB{z}$-gradient,
so a slightly more complicated procedure is necessary. }
For example, it is straightforward to show that the sequence
\begin{equation}
[\FRAC{\pi}{4}(\IX^1+\IX^2)] \rightarrow [\pi\IZ^1\IZ^2]
\rightarrow [-\FRAC{\pi}{6}(\IY^1+\IY^2)]
\end{equation}
applied to the two-spin equilibrium state $\hat{\RHO}_\LAB{eq}$ yields
\begin{equation}
2^{-\frac52} \left( \sqrt3\, \left( \EP^1\EP^2 - \FRAC14 -
\IX^1\IX^2 \right) - \IX^1\EM^2 - \EM^1\IX^2 \right) ~,
\end{equation}
which is reduced by a $\LAB{z}$-gradient
to $(3/32)^{\frac12}(\EP^1\EP^2 - 1/4)$.
Further RF and gradient pulse sequences
which convert the equilibrium state of
two and three spin systems to pseudo-pure states
may be found in Ref.~\cite{CorPriHav:98},
but extending them to larger numbers
of spins is not straightforward.

An alterative proposed by E.~Knill {\em et al.\/}
\cite{KniChuLaf:98} is to ``time-average''
the results of several separate experiments.
In the simple case of a two spin system,
for example, the average of the three states
\begin{equation} \begin{array}{rcl}
\hat{\RHO}_{123} ~~\equiv~& \FRAC{1}{4} (\IZ^1 + \IZ^2)
&\leftrightarrow~~ \FRAC{1}{4}\, \MAT{Diag}( 1, 0, 0, -1 ) \\ 
\hat{\RHO}_{231} ~~\equiv~& \FRAC{1}{4} (\IZ^1 + 2 \IZ^1 \IZ^2)
&\leftrightarrow~~ \FRAC{1}{4}\, \MAT{Diag}( 1, 0, -1, 0 ) \\ 
\hat{\RHO}_{312} ~~\equiv~& \FRAC{1}{4} (2 \IZ^1 \IZ^2 + \IZ^2)
&\leftrightarrow~~ \FRAC{1}{4}\, \MAT{Diag}( 1, -1, 0, 0 )
\end{array} \end{equation}
is the pseudo-pure state
\begin{equation} \begin{array}{rcl} \label{eq:cyclic_avg}
\FRAC{1}{3} (\hat{\RHO}_{123} + \hat{\RHO}_{231} + \hat{\RHO}_{312})
&~=~& \FRAC{1}{12} (2 \IZ^1 + 2 \IZ^2 + 4 \IZ^1 \IZ^2)
\\ &~\leftrightarrow~&
\FRAC{1}{12}\, \MAT{Diag}( 3, -1, -1, -1 ) ~.
\end{array} \end{equation}
More generally, one can obtain the same
results that one would get on a pseudo-pure
state by averaging the results of the calculations
over all $2^N-1$ cyclic permutations of the nonground
state populations of the equilibrium density operator.
Although this naive approach is not efficient,
Knill {\em et al.\/} have shown that one can average over
smaller groups in time $O(N^3)$ with much the same effect.

A fundamentally different approach, first proposed
by N.~Gershenfeld \& \linebreak I.~Chuang \cite{GershChuan:97},
involves working with subpopulations of the molecules
distinguished by the states of additional {\em ancilla\/} spins.
Gershenfeld \& Chuang \cite{ChGeKuLe:98} have given an
example of a two-spin {\em conditional pseudo-pure state\/}
(as we call it), which is obtained by row/column
permutation of the diagonal equilibrium density matrix
$\MAT{Diag}( 3, 1, 1, -1, 1, -1, -1, -3 ) / 16$
of a three-spin system including one ancilla, namely
\begin{equation} \begin{array}{rcl} \label{eq:condpure3a}
&&\FRAC{1}{16}\, \MAT{Diag}( 3, -1, -1, -1, -3, 1, 1, 1 )
\\ &~\leftrightarrow~&
\FRAC{1}{16} (2\IZ^1 (2\IZ^2 + 2\IZ^3 + 4\IZ^2\IZ^3))
\\ &~=~&
\FRAC{1}{4} (\EP^1 - \EM^1) (\EP^2\EP^3 - \FRAC{1}{4}) ~.
\end{array} \end{equation}
The last form makes it clear that in
the subpopulation with the first spin ``up'',
which is labeled by $\EP^1$,
and in the subpopulation with it ``down'',
which is labeled by $\EM^1$,
spins $2$ and $3$ are in the
pseudo-pure state $\EP^2\EP^3 - 1/4$.
Since the spectrum of spins $2$ and $3$ is
antiphase with respect to the ancilla spin $1$,
one can select the subpopulations just by
keeping only either positive or negative peaks.
Although the situation is considerably
more complicated with more spins,
Gershenfeld and Chuang have shown that conditional
pure states can be obtained (with some loss of signal)
using as few as $O(\log(N))$ ancillae;
to our knowledge, the time-space-signal
trade-offs of the various possibilities
have not been fully worked out.

An alternative to conditional pure states,
which we call {\em relative pseudo-pure states\/},
can be obtained via the partial trace operation
(in NMR, decoupling \cite{Slichter:90}),
rather than peak selection as above.
For example, a two-spin relative pseudo-pure state is given
by the partial trace over the ancilla spins 1 \& 2 in
\begin{equation} \begin{array}{rcl} \label{eq:relpure4}
&& \FRAC{1}{32}\, \MAT{Diag}( 4, 2, 2, 0, 2, 0,
-2, 0, 0, -2, 0, 2, 0, -2, -2, -4 )
\\ &\leftrightarrow& \FRAC{1}{16}
\left( \EP^1\EP^2 (\EP^3 + \EP^4) +
\EP^1\EM^2 (\EP^3\EP^4 - \EM^3\EP^4)
\right. \\ &&  ~+ \left.
\EM^1\EP^2 (\EM^3\EM^4 - \EP^3\EM^4) -
\EM^1\EM^2 (\EM^3 + \EM^4) \right) ~,
\end{array} \end{equation}
which is again a permutation of the
diagonal elements of $\hat{\RHO}_\LAB{eq}$.
This can be seen by adding up the $4\times 4$ blocks of
the matrix, obtaining $\MAT{Diag}( 6, -2, -2, -2 ) / 32$.
Alternatively, since the partial trace in the product
operator formalism corresponds to simply eliminating
those terms depending on either spins $1$ or $2$ and
multiplying the remaining terms by $4$ \cite{SomCorHav:98},
we need only add up the multipliers of
$\EP^1\EP^2, \ldots, \EM^1\EM^2$, which yields
\begin{equation} \begin{array}{rcl}
&& \FRAC{1}{16} ((1 + \EP^3 - \EM^3) (1 + \EP^4 - \EM^4) - 1)
\\ &~=~&
\FRAC{1}{4} \left( \EP^3 \EP^4 - \FRAC{1}{4} \right)
~\leftrightarrow~
\FRAC{1}{32}\,\MAT{Diag}( 6, -2, -2, -2 ) ~.
\end{array} \end{equation}

\begin{figure}
\begin{picture}(300,200)
\put(5,-5){ \psfig{file=polarization.eps,width=4.5in}
} \end{picture} \caption{
Negative base ten logarithm of the polarization
$P$ as a function of the logarithm of the ratio
of the energy level spacing to $k_\LAB{B}T$ for
a one, four, sixteen and sixty-four spin
pseudo-pure state obtained by cyclic averaging.
For protons in a standard 500 MHz.\ spectrometer,
$P \sim 10^{-5}$ at equilibrium. }
\end{figure}

We now briefly consider the SNR (signal-to-noise ratio)
of these methods of creating pseudo-pure states.
It has been argued that since the equilibrium population of the
ground state falls off exponentially with the number of spins,
and all these methods are aimed in some fashion at
isolating the signal from the ground state population,
the SNR of all these methods must likewise
decline exponentially with $N$ \cite{Warren:97}.
Although this argument carries considerable weight,
we shall see that the number and variety of the available
methods renders the actual situation more complex.
The standard to which the signal strength must be compared
is that of a single spin in its equilibrium state, namely
\begin{equation}
\hat{\RHO}_\LAB{eq} ~=~ \HALF \IZ^1 ~=~ \FRAC{1}{4}
(\KET{0}\BRA{0} - \KET{1}\BRA{1}) ~.
\end{equation}
The maximum eigenvalue $\| \hat{\RHO}_\LAB{eq} \|_2 = 1/4$
is what we will use as the standard signal strength for
spins of like gyromagnetic ratio (as assumed throughout).
We shall therefore calculate the SNR of a pseudo-pure
state by transforming it to the corresponding ground state
$\KET{0\cdots 0}\BRA{0\cdots 0} - 2^{-N}$ (if need be),
taking the partial trace over all but one of the spins,
and multiplying the maximum eigenvalue of the result by $4$.

The maximum eigenvalue of the partial trace over all
but one of the spins in a pseudo-pure state obtained
by cyclic averaging, as in Eq.~(\ref{eq:cyclic_avg}),
is easily seen to be $N/(4(2^N-1))$, which decays
almost exponentially with the number of spins $N$.
There is an additional factor of $\sqrt{2^N-1}$
which comes from averaging over $2^N-1$ experiments,
and gives a net SNR of $N/(4\sqrt{2^N-1})$ for the average.
The exponential time requirements of cyclic averaging
will nonetheless force one to average over smaller groups,
with consequently smaller improvements in the SNR.
In any case, the SNR declines superpolynomially with $N$.
Figure 1 shows how the signal strength changes as a function
of the ratio of the energy level spacing to $k_\LAB{B}T$,
relative to the signal in a perfectly polarized sample,
when the pseudo-pure state is obtained by cyclic averaging,
for varying numbers of spins.

Because of the many possible variations on the
ideas and the difficulty of analyzing all of them,
it is not practical to present simple formulae for the
SNR of the other methods of preparing pseudo-pure states.
Further complexity is added to the situation by
the ability to combine the various methods above.
A number of such combinations are given in Knill
{\em et al.\/} \cite{KniChuLaf:98},
along with bounds on the SNR for each.
In our laboratory we are developing a new
method, again based on field gradients,
which enables the sample to be divided into
discrete volumes and separate unitary
transformations to be applied to each.
In principle, this permits multiple
experiments to be performed, and their
results added, in a single experiment,
thereby performing an average over
multiple experiments in constant time.
This new method could also be used in a
variety of combinations with existing methods.

It is nevertheless encouraging to observe that
the SNR of the two-spin conditional and relative
pseudo-pure states given in Eqs.~(\ref{eq:condpure3a})
and (\ref{eq:relpure4}) is $1/2$ in both cases;
this is exactly the decline in the ground state
population of a two-spin system compared to a one-spin.
In Eq.~(\ref{eq:condpure3a}), we attain this ``theoretical limit''
because the expansion of the density operator consists of
a single term conditioned on the state of a single ancilla;
it is not possible to do as well with more spins.
In Eq.~(\ref{eq:relpure4}), however, it is
because such permutations are able to transfer
polarization from one set of spins to another.
We have found this makes it possible to derive
a two-spin pseudo-pure state from a six-spin
equilibrium state with {\em no\/} loss of SNR,
whereas a simplistic ground-state population
argument implies we should lose at least $1/2$.
This may be seen by adding up the rows in
the rearrangement of $\hat{\RHO}_\LAB{eq}$
shown in Eq.~(\ref{eq:relpure6}) below,
which corresponds taking the traces of the
four $16\times 16$ blocks along the diagonal,
and yields $\MAT{Diag}( 48, -16, -16, -16 )$.
%\renewcommand{\arraystretch}{1.0}
\begin{equation} \label{eq:relpure6}
\MAT{Diag} \begin{array}[t]{rrrrrrrrrrrrrrrrr} (
~6, &  4, &  4, &  4, &  4, &  4, &  4, &  2, & 
 2, &  2, &  2, &  2, &  2, &  2, &  2, &  2, & \\
 0, &  0, &  0, &  0, &  0, &  0, &  0, &  0, & 
 0, &  0, & -2, & -2, & -2, & -2, & -4, & -4, & \\
 0, &  0, &  0, &  0, &  0, &  0, &  0, &  0, & 
 0, &  0, & -2, & -2, & -2, & -2, & -4, & -4, & \\
~2, & ~2, & ~2, & ~2, & ~2, & ~2, & -2, & -2, & 
-2, & -2, & -2, & -2, & -2, & -4, & -4, & -6\,  &
~) \end{array}
\end{equation}
%\renewcommand{\arraystretch}{1.2}%
The partial trace over one of the two remaining
spins then gives $\MAT{Diag}( 32, -32 )$,
which when divided by $128$ (twice the partition
function) yields $\FRAC{1}{2} \IZ$ as claimed.

A general algorithm has recently been given by
Schulman \& Vazirani \cite{SchulVazir:98} whereby one
can ``distill'' an $M$-spin relative pure state from
an ensemble of molecules each containing $N$ spins.
Starting from a net polarization per spin of $\epsilon$,
this algorithm yields $M$ perfectly polarized
spins providing $M/N \sim O(\epsilon^2)$,
a result anticipated by earlier work in NMR
which showed that the amount of polarization
that can be transferred to a single spin is
at most $O(\sqrt{N})$ \cite{Sorensen:89}.
Unfortunately, given that $\epsilon \approx 10^{-5}$ for
protons at equilibrium in a standard 500 MHz.\ spectrometer,
a molecule with of order $10^{10}$ spins would be needed
to prepare a perfectly polarized state on a single spin
--- which is always in a pseudo-pure state!
The importance of Schulman \& Vazirani's algorithm
thus lies in the fact that it shows that it is at
least theoretically possible to perform fully quantum
computations with only polynomial overhead in space
and time, starting from a thermal equilibrium state.

In summary, although no general and practical method
of preparing pseudo-pure states with bounded SNR from
high-temperature equilibrium states is currently known,
it remains possible that such a method will be found.
It is further worth pointing out that physical
methods of ``refrigerating'' spins are available,
for example optical pumping \cite{NSRATP:96}.
These are presently confined to very simple systems,
but such a source of polarization could in principle be
used in conjunction with polarization transfer techniques
to produce (pseudo-)pure states on large numbers of spins.
Whether or not NMR proves to be a suitable means of building
a large-scale ensemble quantum computer in the future,
some very significant tests of quantum mechanics and
demonstrations of quantum computing are possible with
the small numbers of quantum bits we can achieve today.
To illustrate this, we will now present the results
of NMR experiments which (with certain assumptions)
disprove the existence of ``hidden variables''.

\section{MACROSCOPIC CONSEQUENCES OF QUANTUM CORRELATIONS}
Given the success of the purely classical Bloch
equations (and their multispin extensions) in
describing the NMR phenomena \cite{Slichter:90},
it is perhaps surprising that phenomena widely
regarded as uniquely ``quantum'' are nonetheless
manifest in macroscopic, room-temperature NMR spectroscopy.
For example, Seth Lloyd has recently proposed that
the nonclassical correlations in (Mermin's version of)
the GHZ state can be validated using NMR \cite{Lloyd:98}.
His approach involves using a fourth ``observer'' spin
to perform a nondemolition measurement on the three
spins in a GHZ state (or a pseudo-pure analogue thereof).
Here we shall consider another, rather different way in
which quantum correlations can be demonstrated by NMR.
We stress at the outset that, because all the
interacting spins involved reside in the same
molecule, neither Lloyd's approach nor ours can
establish the nonlocality of the correlations.

The approach taken here was inspired by an educational
paper published a few years ago, in which T.~F.~Jordan
has shown that the contradictions with hidden variables
implied by violations of Bell's inequalities as well as
by the GHZ and Hardy's paradox can be derived entirely by
consideration of the expectation values of product operators,
rather than by observations on single spins \cite{Jordan:94}.
This shows that, in principle, it is not necessary to use
nondemolition measurements with an observer spin in order to
perform experiments which demonstrate these contradictions by NMR;
it can be done directly from observations on ensembles of the spins
of interest, providing at least they are in a (pseudo-)pure state.
In a companion paper to Jordan's, N.~D.~Mermin points out that in
real-life experiments it is nevertheless not possible to perform
the measurements, either of single spins nor (by implication)
of expectation values, with sufficient precision to establish the
``perfect'' (total) correlations on which ``EPR'' arguments against
the existence of hidden variables are based \cite{Mermin:94}.
In that same paper, however, Mermin shows that Hardy's paradox
is a consequence of the Clauser-Horne form of Bell's inequality.
This enables Hardy's paradox \cite{Branning:97,Hardy:92} to be
extended to an open set in the Hilbert space of only two spins,
to which the practically available experimental data can confine us.

In the following, we present the results of NMR experiments
which implement the specific example of Hardy's paradox
presented by Mermin in an Appendix to his paper \cite{Mermin:94}.
Let us map the ``red'' and ``green'' eigenstates $\KET{\LAB{1G}}$
and $\KET{\LAB{1R}}$ of Mermin's measurement $1$ to
the spin states $\KET{0}$ and $\KET{1}$, respectively.
It will be clearer here to relabel this measurement as ``$\LAB{A}$'',
and to use $\KET{\alpha_\LAB{G}} \equiv \KET{0}$ and
$\KET{\alpha_\LAB{R}} \equiv \KET{1}$ as synonyms for its eigenbasis.
Correspondingly, we will relabel Mermin's measurement
$2$ as ``$\LAB{B}$'', and denote its the eigenbasis by
\begin{equation}
\KET{\beta_\LAB{G}} ~\equiv~ \SQRT{\FRAC{3}{5}} \KET{0}
- \SQRT{\FRAC{2}{5}} \KET{1} \quad\mbox{and}\quad
\KET{\beta_\LAB{R}} ~\equiv~ \SQRT{\FRAC{2}{5}} \KET{0}
+ \SQRT{\FRAC{3}{5}} \KET{1} ~.
\end{equation}
Then the state which Mermin has shown leads to a
near-maximum violation of Bell's inequality while
also providing an example of Hardy's paradox is
\begin{equation} \label{eq:sigh}
\KET{\psi} ~\equiv~ \HALF \KET{00} + \SQRT{\FRAC{3}{8}}
\KET{01} + \SQRT{\FRAC{3}{8}} \KET{10} ~.
\end{equation}

To translate this into the context of NMR,
we first note that the observable whose expectation
value is the probability that measurement $\LAB{A}$
yields the state $\KET{\alpha_\LAB{G}}$ is given by
$\VEC A \equiv \EP \equiv \HALF(1+2\IZ)$
(we drop the usual spin index because the measurements
$\LAB{A}$ \& $\LAB{B}$ are assumed the same for both spins).
Similarly the observable which gives the probability that
measurement $\LAB{B}$ yields $\KET{\beta_\LAB{G}}$ is
\begin{equation} \begin{array}{rcl}
\VEC B ~\equiv~ \KET{\beta_\LAB{G}}\BRA{\beta_\LAB{G}}
&~=~& \FRAC{3}{5} \KET{0}\BRA{0} + \FRAC{2}{5} \KET{1}\BRA{1} -
\SQRT{\FRAC{6}{25}} (\KET{0}\BRA{1} + \KET{1}\BRA{0}) \\ 
&~=~& \HALF + \FRAC{1}{5} \IZ - \SQRT{\FRAC{24}{25}} \IX ~.
\end{array} \end{equation}
In addition, the density operator
(including the identity) of Mermin's state is
\begin{equation} \begin{array}{rcl} \label{eq:Sigh}
\EMB\Psi ~\equiv~ \KET{\psi}\BRA{\psi}
&~=~& \FRAC{1}{4} + \FRAC{1}{8} \left(\IZ^1 + \IZ^2\right)
- \HALF \IZ^1\IZ^2 \\ 
&& +\, \SQRT{\FRAC{3}{8}} \left(\IX^1(1 + 2\IZ^2) +
(1 + 2\IZ^1)\IX^2\right) \\ 
&& +\, \FRAC{3}{4}\left( \IX^1\IX^2 + \IY^1\IY^2 \right) ~.
\end{array} \end{equation}

The state $\KET{01}$ is obviously related
to $\KET{00}$ by a rotation of spin $2$,
while $\KET{00}$ can likewise by rotated to
$\KET{10}$, but without affecting $\KET{01}$,
by a {\em conditional\/} rotation of spin $1$.
We shall denote these by
\begin{equation}
\VEC P(\phi) ~\equiv~ e^{-\imath\phi\IY^2}
\quad\mbox{and}\quad
\VEC Q(\theta) ~\equiv~ e^{-\imath\theta\IY^1\EP^2} ~,
\end{equation}
respectively.
They act consecutively on the ground state to yield
\begin{equation}
\BRA{00} \tilde{\VEC P}(\phi) \tilde{\VEC Q}(\theta) ~=~
\left[ \,\cos(\theta/2) \cos(\phi/2),\, \sin(\theta/2)
\cos(\phi/2),\, \sin(\phi/2),\, 0\, \right] ~,
\end{equation}
which is easily verified to equal %the desired vector
$\BRA{\psi} = [ 1/2, \sqrt{3/8}, \sqrt{3/8},\, 0 ]$ when
\begin{equation}
\phi ~=~ 2 \arctan( \SQRT{3/5} )
\quad\mbox{and}\quad
\theta ~=~ 2 \arctan( \SQRT{3/2} ) ~.
\end{equation}
Using the product operator techniques presented in section 2,
these transformations are readily implemented by NMR pulse sequences.

The next thing to notice is that if we take expectation values
with the usual idempotents $\EP^1\EP^2, \ldots, \EM^1\EM^2$, we get
\begin{equation} \begin{array}{rcl}
\FRAC{1}{4} ~=~ & 4\left\langle \EMB\Psi \EP^1\EP^2 \right\rangle
& ~\equiv~ 4\left\langle \EMB\Psi \VEC A^1 \VEC A^2 \right\rangle
\\ 
\FRAC{3}{8} ~=~ & 4\left\langle \EMB\Psi \EP^1\EM^2 \right\rangle
& ~\equiv~ 4\left\langle \EMB\Psi \VEC A^1 (1 - \VEC A^2) \right\rangle
\\ 
\FRAC{3}{8} ~=~ & 4\left\langle \EMB\Psi \EM^1\EP^2 \right\rangle
& ~\equiv~ 4\left\langle \EMB\Psi (1 - \VEC A^1) \VEC A^2 \right\rangle
\\ 
0 ~=~ & 4\left\langle \EMB\Psi \EM^1\EM^2 \right\rangle
& ~\equiv~ 4\left\langle \EMB\Psi (1 - \VEC A^1) (1 - \VEC A^2)
\right\rangle ~.
\end{array} \end{equation}
These correspond to the diagonal of the
density matrix in the usual $\IZ$ basis,
\begin{equation}
{\bf diag}(\EMB\Psi) ~=~
[ \FRAC{1}{4}, \FRAC{3}{8}, \FRAC{3}{8}, 0 ]
\qquad\mbox{($\LAB{A}$ on 1, $\LAB{A}$ on 2)} ~,
\end{equation}
which contains the probabilities of the four possible outcomes
of performing measurement $\LAB{A}$ on both spins (as shown).

The product operator form of $\VEC B$ immediately
makes clear that measurement $\LAB{B}$ is nothing more
or less than the measurement of the magnetization of the
spin along an axis inclined at an angle of $\zeta \equiv
\arctan(\sqrt{24}) = \pi - \theta$ to the $z$-axis in the $xz$-plane.
Letting $\VEC R \equiv \exp(-\imath \zeta \IY)$,
it follows that the probability that measurement $\LAB{B}$ on
spin $1$ yields ``$\LAB{G}$'' (i.e.~$\KET{\beta_\LAB{G}}$) is
\newcommand{\tildeR}% define \tilde{\VEC R} the height of \VEC{R}
{{\rule[0pt]{0pt}{7pt}\smash{\tilde\VEC{R}}}}
\begin{equation}
4 \AVG{ \EMB\Psi \VEC B^1 } ~=~
4 \AVG{ \EMB\Psi \tildeR^1
\VEC A^1 \VEC{R}^1 } ~=~
4 \AVG{ \VEC{R}^1 \EMB\Psi \tildeR^1 \VEC A^1 } ~.
\end{equation}
with a similar expression for spin $2$.
More generally, the probabilities of the outcomes
of the other combinations of measurements are given
by the diagonals of the transformed density matrices:
\begin{equation} \begin{array}{rcl}
{\bf diag}\left( \VEC{R}^2 \EMB\Psi \tildeR^2 \right)
&~=~& [ 0, \FRAC{5}{8}, \FRAC{9}{40}, \FRAC{3}{20} ]
\qquad\qquad\mbox{($\LAB{A}$ on 1, $\LAB{B}$ on 2)} \\
{\bf diag}\left( \VEC{R}^1 \EMB\Psi \tildeR^1 \right)
&~=~& [ 0, \FRAC{9}{40}, \FRAC{5}{8}, \FRAC{3}{20} ]
\qquad\qquad\mbox{($\LAB{B}$ on 1, $\LAB{A}$ on 2)} \\
{\bf diag}\left( \VEC{R}^1\VEC{R}^2
\EMB\Psi \tildeR^2 \tildeR^1 \right)
&~=~& [ \FRAC{9}{100}, \FRAC{27}{200}, \FRAC{27}{200}, \FRAC{16}{25} ]
\qquad\mbox{($\LAB{B}$ on 1, $\LAB{B}$ on 2)}
\end{array} \end{equation}
For compactness, let us denote these probabilities by $\Psi_{kl}^{ij}$,
where $i,j \in \{\LAB{A},\LAB{B}\}$ are the measurements and
$k,l \in \{\LAB{G},\LAB{R}\}$ are the corresponding outcomes,
e.g.~$\Psi_\LAB{GR}^\LAB{AB} = 4\AVG{\EMB\Psi \VEC A^1(1 - \VEC B^2)}$.

We may translate Mermin's proof \cite{Mermin:94}
that these probabilities are incompatible
with hidden variables associated with the
individual spins into this context as follows:
First, since $\Psi_\LAB{GG}^\LAB{AB} = \Psi_\LAB{GG}^\LAB{BA} = 0$,
in any molecule wherein one of the spins is parallel
to the $z$-axis the other must be antiparallel to
the axis of measurement $\LAB{B}$ and vice versa.
Hence, since $\Psi_\LAB{GG}^\LAB{BB}$ is nonzero,
in some molecules ($9$\%, to be precise) both
spins must be antiparallel to the $z$-axis.
But this contradicts the fact that $\Psi_\LAB{RR}^\LAB{AA} = 0$.
More generally, Mermin has shown that
\begin{equation}
\Psi_\LAB{GG}^\LAB{BB} ~\le~ \Psi_\LAB{GG}^\LAB{AB}
+ \Psi_\LAB{RR}^\LAB{AA} + \Psi_\LAB{GG}^\LAB{BA}
\end{equation}
is an example of the Clauser-Horne form of Bell's inequality.
Thus to disprove the existence of such one-particle hidden variables
it is sufficient to determine these probabilities to $\pm2$\% or so.

At this point we encounter a significant complication,
which is that the ``strong'' (von Neumann) measurements
assumed in their analyses by Jordan and
Mermin {\em cannot\/} be implemented by NMR;
we can only perform ``weak'' (nonperturbing) measurements
of the {\em relative population differences\/} between
states connected by single spin flips \cite{CorPriHav:98}.
This is done by applying a magnetic field gradient
along the $\LAB{z}$-axis, which (as previously described)
dephases any transverse components in the density operator.
Thereafter, a pair of ``soft'' $\pi/2$ readout pulses,
each tuned to the frequency of just one of the two spins,
produces a pair of spectra each with two peaks whose
heights are proportional to the population differences
between pairs of states connected by flips of that spin.
The factor relating the peak heights to the
corresponding differences in the {\em probabilities\/}
of the states can be determined from spectra
collected on the pseudo-pure ground state,
after which it is straightforward to convert
the differences into the corresponding
absolute probabilities by linear least squares,
subject to the constraint that their sum is unity.
We shall encounter field gradients again in the next section,
when we show how they can also be used to implement
precisely controlled decoherence models.

Thus the overall experiment consists of
collecting ten spectra, as follows:
\begin{enumerate}
\item Prepare the state $\EMB\Psi$, by first
preparing the pseudo-pure ground state $\KET{00}$
using one of the previously described methods,
and then transforming it by $\VEC Q(\theta) \VEC P(\phi)$.
\item Use a selective radio-frequency pulse to
apply the rotation $\VEC R$ to those spins on
which measurement $\LAB{B}$ is to be performed.
\item Use a $\LAB{z}$-gradient to dephase the transverse
components of the resulting density operator.
\item Apply a readout pulse to one of the
spins, and collect the corresponding spectrum;
repeat steps 1 -- 3 and then do the same for the other spin.
\item Repeat steps 1 - 4 for each of the four combinations of
measurements $\LAB{AA}$, $\LAB{AB}$, $\LAB{BA}$ and $\LAB{BB}$
on the two spins.
\item Collect two additional calibration spectra by applying soft
readout pulses to each spin in the pseudo-pure ground state.
\end{enumerate}
These experiments were performed on a Bruker
400 MHz.~spectrometer using the two spin $\HALF$
nuclei in ${}^{13}\LAB{C}$-labeled chloroform.

\begin{figure}
\begin{picture}(324,275)
\put(0,144){ \psfig{file=tseng_aa.eps,width=2.25in} }
\put(75,255){ $\LAB{A^CA^H}$ }
\put(170,144){ \psfig{file=tseng_ab.eps,width=2.25in} }
\put(245,255){ $\LAB{A^CB^H}$ }
\put(0,0){ \psfig{file=tseng_ba.eps,width=2.25in} }
\put(75,110){ $\LAB{B^CA^H}$ }
\put(170,0){ \psfig{file=tseng_bb.eps,width=2.25in} }
\put(245,110){ $\LAB{B^CB^H}$ }
\end{picture} \caption{
The pairs of ${}^{13}\LAB{C}$-labeled chloroform
spectra (carbon left, proton right) obtained by
performing the four combinations of measurements
$\LAB{A^CA^H}$, $\LAB{A^CB^H}$, $\LAB{B^CA^H}$ and $\LAB{B^CB^H}$
on the pseudo-pure form of Mermin's state $\EMB\Psi$.
The spectra have been normalized by the height of
the peak of the corresponding spin in the pseudo-pure
ground state, and the horizontal axis is in kHz.
The transitions of the peaks, from left to right, are
$\KET{0^\LAB{C}0^\LAB{H}} \leftrightarrow \KET{1^\LAB{C}0^\LAB{H}}$,
$\KET{0^\LAB{C}1^\LAB{H}} \leftrightarrow \KET{1^\LAB{C}1^\LAB{H}}$,
$\KET{0^\LAB{C}0^\LAB{H}} \leftrightarrow \KET{0^\LAB{C}1^\LAB{H}}$,
$\KET{1^\LAB{C}0^\LAB{H}} \leftrightarrow \KET{1^\LAB{C}1^\LAB{H}}$.}
\label{fig:hardy}
\end{figure}

The spectra obtained from steps 1 -- 5 of this
experiment are shown in Fig.~\ref{fig:hardy}.
The probabilities derived from these peak heights,
and the residual associated with each,
are shown in the table below.
%
\begin{table} \begin{center}
%\renewcommand{\arraystretch}{1.0}
\begin{tabular}{|c|rrrr|c|} \hline
\multicolumn{6}{|c|}{\parbox[t]{3.625in}{\bf
Probabilities of Outcomes {\boldmath$\LAB{G}$}
\& {\boldmath$\LAB{R}$} for the Measurements
{\boldmath$\LAB{A}$} \& {\boldmath$\LAB{B}$}
(Carbon, Proton) Demonstrating Hardy's Paradox,
as Derived from the Chloroform NMR Spectra in
Fig.\ \ref{fig:hardy} \vspace{3pt}}
} \\ \hline\hline
~Measurements~ & $\quad(\LAB{G},\LAB{G})$ & $\qquad(\LAB{G},\LAB{R})$ &
$\qquad(\LAB{R},\LAB{G})$ & $\qquad(\LAB{R},\LAB{R})\quad$ & ~Residuals~
\\ \hline
$(\LAB{A^C},\LAB{A^H})$ & $0.253$ & $0.380$ & $0.366$ & $0.001\quad$
& $0.008$  \\
$(\LAB{A^C},\LAB{B^H})$ & $0.029$ & $0.609$ & $0.217$ & $0.145\quad$
& $0.018$  \\
$(\LAB{B^C},\LAB{A^H})$ & $-0.002$ & $0.230$ & $0.614$ & $0.159\quad$
& $0.005$ \\
$(\LAB{B^C},\LAB{B^H})$ & $0.097$ & $0.125$ & $0.156$ & $0.622\quad$
& $0.021$ \\ \hline
\end{tabular}
%\renewcommand{\arraystretch}{1.2}
\vspace{-10pt}
\end{center} \end{table}
%
It follows that Bell's inequality is violated by
\begin{equation} \begin{array}{rl}
& \Psi_\LAB{GG}^\LAB{AB} + \Psi_\LAB{RR}^\LAB{AA} +
\Psi_\LAB{GG}^\LAB{BA} - \Psi_\LAB{GG}^\LAB{BB} \\
=~ & 0.029 + 0.001 - 0.002 - 0.097 ~=~ -0.069 ~.
\end{array} \end{equation}
A rigorous error analysis is not possible,
because the dominant errors in NMR spectra
(e.g.~RF field inhomogeneity) are not statistical.
If we nevertheless take the mean RMS residual
(half the total residuals shown in the table)
of $0.0065$ as an estimate of the errors and
assume they are independent between spectra,
the expected error in the sum of these four
numbers is only $0.013$, so that this violation
of Bell's inequality appears significant.

It is important to point out that,
because these measurements are performed
on a weakly polarized pseudo-pure state
rather than on a true pure state,
they are subject to interpretation.
If one were to pull the molecules out of the
sample one at a time and perform the measurements
$\LAB{A}$ and $\LAB{B}$ with a Stern-Gerlach apparatus,
the resulting probabilities would all be very close to $1/4$,
which would {\em not\/} violate Bell's inequality.
This is because the vast majority of the density
operator is contained in its identity component,
which however does {\em not\/} contribute to the NMR spectra.

Indeed, a recent preprint by Braustein, Caves,
Jozsa, Linden, Popescu and Shack claims that
since the total density operator of such a
weakly polarized sample {\em can\/} be explained
by ensemble of unentangled spin systems,
such experiments should be regarded as
``simulations'' of true quantum systems and
not instances thereof \cite{BrausteinEtAl:98}.
Nevertheless, this is not the only valid
definition of mixed state entanglement
(as Braustein {\em et al.\/} also point out).
For example, it is well-known that one can
teleport an unknown mixed state with better
fidelity than could be achieved by local
operations and classical communication,
even though the state does not violate any Bell
inequality (see exercise 5.5 in \cite{Preskill:98}).

In contrast to Braustein {\em et al.}'s point of view,
our claim to have obtained a violation of Bell's
inequality is based on the {\em assumption\/}
that the spins of each molecule in a pseudo-pure
ensemble are in a definite, albeit unknown,
pure state, and that the molecules contributing
to the spectra of such an ensemble are all
in the {\em same\/} corresponding pure state.
The arbitrary choice of ensemble for a density operator
is in many respects similar to the arbitrary choice of
a coordinate system in physical problems more generally,
and our ensemble interpretation for that part of a
pseudo-pure density operator which can be observed
by NMR is arguably the simplest, in that it involves
the smallest possible number of pure states: one.
It should further be clear that, because the
unobserved molecules are in a different state
from those which contribute to the spectra,
this interpretation of our experiments is
quite different from the ``loophole'' in
the Aspect experiments \cite{AspectEtAl:82},
where the unobserved photon pairs are expected to
follow the same statistics as the observed ones.

This microscopic interpretation of pseudo-pure density
operators, which in one way or another underlies
all existing claims that it is possible to observe
quantum phenomena by macroscopic NMR spectroscopy,
will be further justified in future publications.
For now, we will simply point out that the interpretation
is at least supported by the otherwise incredible
ability of NMR spectroscopy on pseudo-pure states
to reproduce all manner of quantum phenomena.
To further emphasize this, we will now present
experimental results which demonstrate the principles
of quantum error correction from an NMR perspective.

\section{QUANTUM ERROR CORRECTION BY NMR SPECTROSCOPY}
The error correcting code we have chosen to illustrate
by NMR is well-known in the field \cite{KnillLafla:97},
and uses two ancilla (labeled $2$ \& $3$) to
encode the state of a data spin (labeled $1$).
Letting $\VEC{S}^{2|1}$ and $\VEC{S}^{3|1}$ be
c-NOT's, and $\VEC{R}_{90}^{123} \equiv
\exp(-\imath\FRAC{\pi}{2}(\IY^1+\IY^2+\IY^3))$,
the encoding operation proceeds as follows:
\begin{equation} \begin{array}{rcl}
&& (\alpha\KET{0} + \beta\KET{1}) \KET{00} ~
\stackrel{\VEC{S}^{2|1}}{\longrightarrow}
\stackrel{\VEC{S}^{3|1}}{\longrightarrow}
\stackrel{\VEC{R}_{90}^{123}}{\longrightarrow}
~ \alpha \KET{\mbox{$+$$+$$+$}} + \beta\KET{\mbox{$-$$-$$-$}}
\\ &&
\left( \mbox{where}~ \KET{\mbox{$\pm$$\pm$$\pm$}} \equiv
(\KET{0}\pm\KET{1}) (\KET{0}\pm\KET{1})(\KET{0}\pm\KET{1}) \right)
\end{array} \end{equation}
Decoding consists of applying the
inverse operations in the reverse order,
which acts on the states obtained by
single sign-flip errors as follows:
\begin{equation} \begin{array}{rcl}
&&
\alpha \KET{\mbox{$+$$+$$-$}} + \beta\KET{\mbox{$-$$-$$+$}}
\stackrel{\VEC{R}_{-90}^{123}}{\longrightarrow}
\stackrel{\VEC{S}^{3|1}}{\longrightarrow}
\stackrel{\VEC{S}^{2|1}}{\longrightarrow}
(\alpha\KET{0} + \beta\KET{1})\KET{01}
\\ &&
\alpha \KET{\mbox{$+$$-$$+$}} + \beta\KET{\mbox{$-$$+$$-$}}
\stackrel{\VEC{R}_{-90}^{123}}{\longrightarrow}
\stackrel{\VEC{S}^{3|1}}{\longrightarrow}
\stackrel{\VEC{S}^{2|1}}{\longrightarrow}
(\alpha\KET{0} + \beta\KET{1})\KET{10}
\\ &&
\alpha \KET{\mbox{$-$$+$$+$}} + \beta\KET{\mbox{$+$$-$$-$}}
\stackrel{\VEC{R}_{-90}^{123}}{\longrightarrow}
\stackrel{\VEC{S}^{3|1}}{\longrightarrow}
\stackrel{\VEC{S}^{2|1}}{\longrightarrow}
(\alpha\KET{1} + \beta\KET{0})\KET{11}
\end{array} \end{equation}
It follows that a Toffli gate $\VEC{T}^{1|23}$,
which flips the data spin conditional on
the ancillae being in the state $\KET{11}$,
will correct a sign-flip error in the
data spin and leave it alone otherwise,
even if an error occurs in the ancillae.

In practice, errors in quantum computers
are not expected to be single sign-flips,
but rather small random phase errors
which cumulatively result in decoherence.
Nevertheless, we can show that the ability to
correct sign-flips implies the ability to cancel
the effect of such phase errors to first order.
Random phase errors correspond to the propagator
$\exp(-\imath(\chi^1\IZ^1+\chi^2\IZ^2+\chi^3\IZ^3))$,
which acts to first order on the encoded state as:
\begin{equation} \begin{array}{rcl} &&
\exp(-\imath(\chi^1\IZ^1+\chi^2\IZ^2+\chi^3\IZ^3) \left(
\alpha \KET{\mbox{$+$$+$$+$}} + \beta\KET{\mbox{$-$$-$$-$}}
\right) \\ &~\approx~& \left(
\alpha \KET{\mbox{$+$$+$$+$}} + \beta\KET{\mbox{$-$$-$$-$}}
\right) - \imath\chi^1 \left(
\alpha \KET{\mbox{$-$$+$$+$}} + \beta\KET{\mbox{$+$$-$$-$}}
\right) \\ && -\, \imath\chi^2 \left(
\alpha \KET{\mbox{$+$$-$$+$}} + \beta\KET{\mbox{$-$$+$$-$}}
\right) - \imath\chi^3 \left(
\alpha \KET{\mbox{$+$$+$$-$}} + \beta\KET{\mbox{$-$$-$$+$}}
\right)
\end{array} \end{equation}
Since decoding and the error-correcting
Toffoli gate are likewise linear,
it follows that the first-order effects
of phase errors are cancelled as claimed.
Note this argument makes no assumptions
concerning the correlations among the errors!

Experimental results demonstrating these expectations
have recently been published \cite{CMPKLZHS:98}.
In the following, we shall present a more
detailed explanation of how the error correction
works using the product operator formalism,
along with selected experimental data
illustrating and validating this explanation.
We shall assume that the data spin is in one of the
states $1$ (unpolarized), $\IX^1$, $\IY^1$ or $\IZ^1$.
Although these are mixed states,
each consists of an incoherent sum of pure states,
e.g.\ $2\IZ^1 = \KET{0}\BRA{0} - \KET{1}\BRA{1}$,
so if error correction works on these pure states,
by linearity it will also work on the mixtures (and vice versa).
In these terms, a complete set of initial states
$\hat{\RHO}_\LAB{A}$ for error correction are:
\begin{equation} \begin{array}{rcl} \label{eq:rhoA}
&& \left. \begin{array}{l}
\EP^2\EP^3 \\ \IX^1 \EP^2\EP^3 \\
\IY^1 \EP^2\EP^3 \\ \IZ^1 \EP^2\EP^3
\end{array} \right\}
~\equiv~ \hat{\RHO}_\LAB{A}^1 \EP^2 \EP^3 ~=~ \hat{\RHO}_\LAB{A}
\end{array} \end{equation}
The corresponding states $\hat{\RHO}_\LAB{B}$
to which they are mapped by encoding are:
\begin{equation} \begin{array}{rcl} \label{eq:rhoB}
\hat{\RHO}_\LAB{B} &~\equiv~& \left\{ \begin{array}{l}
\FRAC14 + \IX^1\IX^2 + \IX^1\IX^3 + \IX^2\IX^3 \\
\IZ^1\IY^2\IY^3 + \IY^1\IZ^2\IY^3 + \IY^1\IY^2\IZ^3 - \IZ^1\IZ^2\IZ^3 \\
\IZ^1\IZ^2\IY^3 + \IZ^1\IY^2\IZ^3 + \IY^1\IZ^2\IZ^3 - \IY^1\IY^2\IY^3 \\
\FRAC{1}{4} (\IX^1 + \IX^2 + \IX^3) + \IX^1\IX^2\IX^3
\end{array} \right.
\end{array} \end{equation}
We note the last three states in Eq.~(\ref{eq:rhoA})
can be prepared (with a 50\% loss of polarization)
from the average of twice Eq.~(\ref{eq:condpure3a}) with
\begin{equation} \begin{array}{rcl} \label{eq:condpure3b}
&& \FRAC{1}{16}\, \MAT{Diag}( 3, 1, 1, 1, -3, -1, -1, -1 )
\\ &~\leftrightarrow~&
\FRAC{1}{16} (\IZ^1 (3 + 2\IZ^2 + 2\IZ^3 + 4\IZ^2\IZ^3))
\\ &~=~&
\FRAC{1}{4} (\EP^1 - \EM^1) (\EP^2\EP^3 + \FRAC{1}{2}) ~.
\end{array} \end{equation}

In NMR, decoherence occurs principally through
the randomly fluctuating external magnetic fields
at each spin $k$ along the $\LAB{z}$-axis, $B_\LAB{z}^k$.
The effect of these fields is most simply described
in the {\em spherical\/} operator basis $1$,
$\IZ^k$ and $\IPM^k \equiv \IX^k \pm \imath \IY^k$,
as opposed to the {\em Cartesian\/} basis used up to now.
The products of these basis elements can be shown
to decay exponentially at rates proportional to the
mean-square field $\OL{(B_\LAB{z}^k)^2}$ for $\IPM^k$
(as well as $\IPM^k \IZ^\ell$, $\IPM^k \IZ^\ell \IZ^m$),
and to
\begin{equation} \label{eq:fields}
\begin{array}[t]{rll}
& \OL{(B_\LAB{z}^k - B_\LAB{z}^\ell)^2} &
\quad\mbox{for $\IP^k\IM^\ell$ \& $\IM^k\IP^\ell$,} \\
& \OL{(B_\LAB{z}^k + B_\LAB{z}^\ell)^2} &
\quad\mbox{for $\IP^k\IP^\ell$ \& $\IM^k\IM^\ell$,} \\ 
& \OL{(B_\LAB{z}^k + B_\LAB{z}^\ell - B_\LAB{z}^m)^2} &
\quad\mbox{for $\IP^k\IP^\ell\IM^m$ \& $\IM^k\IM^\ell\IP^m$,~~etc.,} \\
\mbox{and}\quad & \OL{(B_\LAB{z}^k + B_\LAB{z}^\ell + B_\LAB{z}^m)^2} &
\quad\mbox{for $\IP^k\IP^\ell\IP^m$ \& $\IM^k\IM^\ell\IM^m$.}
\end{array}
\end{equation}
These products are referred to as single (SQC1: $\IPM^k$),
zero (ZQC: $\IPM^k\IMP^\ell$), double (DQC: $\IPM^k\IPM^\ell$),
three-spin single (SQC3: $\IPM^k\IPM^\ell\IMP^m$, etc.) and triple
(TQC: $\IPM^k\IPM^\ell\IPM^m$) quantum coherences, respectively.

We shall consider two extreme forms of decoherence.
In the first, the fields at the different spins are
uncorrelated, and hence the random variables $\chi^k$ can
be assumed to be identically distributed and independent.
In the second, they are assumed to be totally correlated.
By Eq.\ (\ref{eq:fields}), the relative rates
of decoherence in these two cases are:
\begin{equation}
\begin{array}{lccccc}
& ~\mbox{\sf ZQC}~ & ~\mbox{\sf SQC1}~
& ~\mbox{\sf SQC3}~ & ~\mbox{\sf DQC}~
& ~\mbox{\sf TQC}~ \\
\mbox{\sf Uncorrelated:}
& 2 & 1 & 3 & 2 & 3 \\
\mbox{\sf Totally Correlated:}
& 0 & 1 & 1 & 4 & 9
\end{array}
\end{equation}
Decomposing $\hat{\RHO}_\LAB{B}$ into a spherical basis,
multiplying by decaying exponentials with the above
rates normalized by the SQC1 decay rate $\tau$,
and returning to the Cartesian basis gives
\begin{equation}
\hat{\RHO}_\LAB{C} ~\equiv~ \left\{ \begin{array}{l}
\FRAC14 + (\IX^1\IX^2 + \IX^1\IX^3 + \IX^2\IX^3) e^{-2t/\tau}
\skiplinehalf \\
(\IZ^1\IY^2\IY^3 + \IY^1\IZ^2\IY^3 + \IY^1\IY^2\IZ^3) e^{-2t/\tau}
- \IZ^1\IZ^2\IZ^3
\skiplinehalf \\
(\IZ^1\IZ^2\IY^3 + \IZ^1\IY^2\IZ^3 + \IY^1\IZ^2\IZ^3) e^{-t/\tau}
- \IY^1\IY^2\IY^3 e^{-3t/\tau}
\skiplinehalf \\
\FRAC{1}{4} (\IX^1 + \IX^2 + \IX^3) e^{-t/\tau}
+ \IX^1\IX^2\IX^3 e^{-3t/\tau}
\end{array} \right.
\end{equation}
in the uncorrelated case, and
\begin{equation}
\hat{\RHO}_\LAB{C} ~\equiv~
\left\{ \begin{array}{l}
\FRAC14 + \HALF (\IX^1\IX^2 + \IX^1\IX^3 + \IX^2\IX^3 +
\IY^1\IY^2 + \IY^1\IY^3 + \IY^2\IY^3) \\ +\,
\HALF (\IX^1\IX^2 + \IX^1\IX^3 + \IX^2\IX^3 - \IY^1\IY^2
- \IY^1\IY^3 - \IY^2\IY^3) e^{-4t/\tau}
\skiplinehalf \\
\HALF (\IZ^1(\IX^2\IX^3+\IY^2\IY^3) +
\IZ^2(\IX^1\IX^3+\IY^1\IY^3) +
\IZ^3(\IX^1\IX^2+\IY^1\IY^2)) \\ -\,
\HALF (\IZ^1(\IX^2\IX^3-\IY^2\IY^3) +
\IZ^2(\IX^1\IX^3-\IY^1\IY^3) +
\IZ^3(\IX^1\IX^2-\IY^1\IY^2)) \\ \qquad
e^{-4t/\tau} - \IZ^1\IZ^2\IZ^3
\skiplinehalf \\
(\IZ^1\IZ^2\IY^3 + \IZ^1\IY^2\IZ^3 + \IY^1\IZ^2\IZ^3)
e^{-t/\tau} \\ +\, \FRAC14
(\IX^1\IX^2\IY^3 + \IX^1\IY^2\IX^3 + \IY^1\IX^2\IX^3 + 3\IY^1\IY^2\IY^3)
e^{-t/\tau} \\ -\, \FRAC14
(\IX^1\IX^2\IY^3 + \IX^1\IY^2\IX^3 + \IY^1\IX^2\IX^3 - \IY^1\IY^2\IY^3)
e^{-9t/\tau}
\skiplinehalf \\
\FRAC14 (\IX^1 + \IX^2 + \IX^3)
e^{-t/\tau} \\ +\, \FRAC14
(3\IX^1\IX^2\IX^3 + \IY^1\IY^2\IX^3 + \IY^1\IX^2\IY^3 + \IX^1\IY^2\IY^3)
e^{-t/\tau} \\ +\, \FRAC14
(\IX^1\IX^2\IX^3 - \IY^1\IY^2\IX^3 - \IY^1\IX^2\IY^3 - \IX^1\IY^2\IY^3)
e^{-9t/\tau}
\end{array} \right.
\end{equation}
in the totally correlated.
The decoding operation converts this to
\begin{equation}
\hat{\RHO}_\LAB{D} ~\equiv~ \left\{ \begin{array}{l}
\FRAC14 + (\IZ^2 + \IZ^3 + \IZ^2\IZ^3) e^{-2t/\tau}
\skiplinehalf \\
(\HALF \IX^1\IZ^2 + \HALF \IX^1\IZ^3 + \IX^1\IZ^2\IZ^3) e^{-2t/\tau}
+ \FRAC{1}{4} \IX^1
\skiplinehalf \\
(\FRAC{1}{4} \IY^1 + \HALF \IY^1\IZ^2 + \HALF \IY^1\IZ^3) e^{-t/\tau}
+ \IY^1\IZ^2\IZ^3 e^{-3t/\tau}
\skiplinehalf \\
(\FRAC{1}{4} \IZ^1 + \HALF \IZ^1\IZ^2 + \HALF \IZ^1\IZ^3) e^{-t/\tau}
+ \IZ^1\IZ^2\IZ^3 e^{-3t/\tau}
\end{array} \right.
\end{equation}
in the uncorrelated case, and
\begin{equation}
\hat{\RHO}_\LAB{D} ~\equiv~
\left\{ \begin{array}{l}
\FRAC14 + \HALF (\HALF\IZ^2 + \HALF\IZ^3 + \IZ^2\IZ^3 -
2\IX^1\IZ^2\IY^3 - 2\IX^1\IY^2\IZ^3 + \IY^2\IY^3) \\ +\,
\HALF (\HALF\IZ^2 + \HALF\IZ^3 + \IZ^2\IZ^3 +
2\IX^1\IZ^2\IY^3 + 2\IX^1\IY^2\IZ^3 - \IY^2\IY^3)
e^{-4t/\tau} 
\skiplinehalf \\
\HALF (\IX^1(\IY^2\IY^3+\IZ^2\IZ^3) +
\HALF \IZ^3(\IX^1-\IX^2) +
\HALF \IZ^2(\IX^1-\IX^3))
+ \FRAC{1}{4} \IX^1 \\ +\,
\HALF (\IX^1(\IY^2\IY^3-\IZ^2\IZ^3) +
\HALF \IZ^3(\IX^1+\IX^2) +
\HALF \IZ^2(\IX^1+\IX^3))
e^{-4t/\tau}
\skiplinehalf \\
(\HALF\IY^1\IZ^3 + \HALF\IY^1\IZ^2 + \FRAC{1}{4}\IY^1)
e^{-t/\tau} \\ +\, \FRAC14
(\IZ^1\IZ^2\IY^3 + \IZ^1\IY^2\IZ^3 - \IY^1\IY^2\IY^3 - 3\IY^1\IZ^2\IZ^3)
e^{-t/\tau} \\ -\, \FRAC14
(\IZ^1\IZ^2\IY^3 + \IZ^1\IY^2\IZ^3 - \IY^1\IY^2\IY^3 + \IY^1\IZ^2\IZ^3)
e^{-9t/\tau}
\skiplinehalf \\
\FRAC14 (\IZ^1 + 2 \IZ^1\IZ^2 + 2 \IZ^1\IZ^3)
e^{-t/\tau} \\ +\, \FRAC14
(3\IZ^1\IZ^2\IZ^3 + \IY^1\IY^2\IZ^3 + \IY^1\IZ^2\IY^3 + \IZ^1\IY^2\IY^3)
e^{-t/\tau} \\ +\, \FRAC14
(\IZ^1\IZ^2\IZ^3 - \IY^1\IY^2\IZ^3 - \IY^1\IZ^2\IY^3 - \IZ^1\IY^2\IY^3)
e^{-9t/\tau}
\end{array} \right.
\end{equation}
in the totally correlated.
This is clearly getting a little messy,
and it gets much worse after the Toffoli gate!
Therefore, we shall only present the partial trace
over the ancillae after applying the Toffoli, which is
\begin{equation}
\hat{\RHO}_\LAB{E}^1 ~\equiv~ \left\{ \begin{array}{l}
1 \\
\IX^1 \\
\IY^1 \left( \FRAC{3}{2} e^{-t/\tau} - \HALF e^{-3t/\tau} \right)
\\
\IZ^1 \left( \FRAC{3}{2} e^{-t/\tau} - \HALF e^{-3t/\tau} \right)
\end{array} \right.
\end{equation}
in the uncorrelated case, and
\begin{equation}
\RHO_\LAB{E}^1 ~\equiv~
\left\{ \begin{array}{l}
1 \\
\IX^1 \\
\IY^1 \left( \FRAC{3}{2} e^{-t/\tau} - \FRAC{3}{8} e^{-t/\tau}
- \FRAC{1}{8} e^{-9t/\tau} \right)
\\
\IZ^1 \left( \FRAC{3}{2} e^{-t/\tau} - \FRAC{3}{8} e^{-t/\tau}
- \FRAC{1}{8} e^{-9t/\tau} \right)
\end{array} \right.
\end{equation}
in the correlated.
The slope of these curves at $t = 0$ is zero in all cases, as expected.

In order to demonstrate these results by NMR solution-state spectroscopy,
a precise implementation of the above decoherence models is needed.
This was achieved by using {\em gradient-diffusion\/} methods.
In these methods, a magnetic field gradient
is created along the $\LAB{z}$-axis;
as previously described, this dephases
the transverse ($\LAB{xy}$) magnetization.
More precisely, a field gradient causes the
transverse magnetization to precess at rates
which depend linearly on its $\LAB{z}$-coordinate,
thereby winding it into a spiral about the $\LAB{z}$-axis
whose average transverse magnetization is essentially zero.
The gradient is turned off for a given time interval $t$, during
which diffusion of the molecules along $\LAB{z}$ blurs the spiral.
The gradient is then reversed, causing the
magnetization to refocus and so create an ``echo''.
Because those molecules which have moved now precess at
a different rate, their magnetization is not refocussed,
so the magnitude of the echo decays exponentially with $t$.
Because all the spins in each molecule are subject
to the same change in field, this constitutes a
true implementation of the totally correlated model.
Using radio-freqency gradients, it is also possible to dephase
each spin separately, thereby implementing the uncorrelated model.
At this time, however, we have collected and processed data only for
the $\hat{\RHO}_\LAB{A}^1 = \IZ^1$ state with the totally correlated model.

\begin{figure}
\begin{picture}(300,230)
\put(5,-5){ \psfig{file=experiment.eps,width=4.5in}%,height=4.5in}
} \end{picture} \caption{
Experimental NMR data illustrating the decay of each of the product
% EXPANDING THE FOLLOWING MACROS PRODUCES LINE TOO LONG IN AUX FILE (DAMN IT!)
% operators $\IZ^1$, $2\IZ^1\IZ^2$, $2\IZ^1\IZ^3$ and $4\IZ^1\IZ^2\IZ^3$
operators $\mbox{\boldmath$I$}_{\sf z}^1$,
$2\mbox{\boldmath$I$}_{\sf z}^1\mbox{\boldmath$I$}_{\sf z}^2$,
$2\mbox{\boldmath$I$}_{\sf z}^1\mbox{\boldmath$I$}_{\sf z}^3$,
$4\mbox{\boldmath$I$}_{\sf z}^1\mbox{\boldmath$I$}_{\sf z}^2
\mbox{\boldmath$I$}_{\sf z}^3$,
as functions of the time allowed for gradient diffusion (see text),
together with the least-squares fits to their logarithms.
The single and triple quantum coherences in
% $4\IZ^1\IZ^2\IZ^3$	<-- DITTO!
$4\mbox{\boldmath$I$}_{\sf z}^1\mbox{\boldmath$I$}_{\sf z}^2
\mbox{\boldmath$I$}_{\sf z}^3$,
(negative curves) have been plotted and fit separately.
The sum of these data and the fits are also shown (topmost curve),
which illustrates that error correction cancels the decay
of the encoded state to first order as expected.
\vspace*{-12pt} }
\end{figure}

Although it is possible to prepare the
state $\IZ^1\EP^2\EP^3$ as noted above,
we have chosen to illustrate the above analysis by
preparing the states $\IZ^1$, $2\IZ^1\IZ^2$, $2\IZ^1\IZ^3$
and $4\IZ^1\IZ^2\IZ^3$ in four separate experiments,
each using sixteen different decoherence times $t$.
Because the SQC and TQC contributing to $4\IZ^1\IZ^2\IZ^3$
refocussed at different times, this further
enabled us to follow their evolutions separately.
The results of these experiments are
plotted against the time $t$ in Figure 2,
along with the corresponding logarithmic fits.
It may be seen that the sum of the data
and the fits thereto (also shown) do
indeed exhibit a near-zero initial slope,
in accord with the above calculations.
Our published report \cite{CMPKLZHS:98}
includes the results of further experiments
(performed by E.~Knill and R.~Laflamme)
with the natural and far more complicated
decoherence processes that occur in solution.
These are more difficult to interpret,
but are nevertheless consistent with the state
preservation expected from error correction.
Additional experiments and more detailed
calculations are in progress.

While a method of inhibiting decoherence ($T_2$ relaxation)
during NMR pulse sequences would be highly desirable,
there are strong reasons to doubt that quantum
error correction will be useful in this regard.
First, the ancillae must be placed in a pseudo-pure state,
which as we have shown above entails a loss
of 50\% of the signal for each ``data'' spin;
this is more than is recovered by error correction.
In addition, the ancillae must be returned to a pseudo-pure
state uncorrelated with the state of the data spin(s), or else
``fresh'' ancillae in such a state must be continuously available,
in order to inhibit decoherence over an appreciable
period of time by the repeated correction of errors.
Nevertheless, we feel that the basic idea underlying
error correction of preparing multiple quantum coherences,
allowing them to decoher, and then mixing them so
as to determine their relative rates of relaxation,
may be of considerable use in NMR studies
of the statistics of molecular motion.
This in turn is one of the most important
applications of NMR spectroscopy.
Conversely, whereas NMR spectroscopists have
previously used their methods solely to unravel
the secrets of naturally occurring systems,
it now appears possible to use these same
methods to engineer artificial systems in
which the basic principles of quantum mechanics
and the emergence of the classical world through
decoherence can be studied in unprecedented detail
\cite{GuiliniEtAl:96}.

\section*{ACKNOWLEDGEMENTS}
We thank E.~Knill and R.~Laflamme of
Los Alamos National Labs for teaching
us about quantum error correcting codes,
and S.~Braunstein and R.~Jozsa for useful
discussions on mixed-state correlation.
This work was supported by the U.~S.~Army Research
Office under grant number DAAG 55-97-1-0342
from the DARPA Ultrascale Computing Program.

\bibliography{../../math,../../csci,../../nmr,../../phys,../../self}
\bibliographystyle{plain}
\end{document}