Spin qubits for quantum simulations

Abstract

The investigation of quantum mechanical systems mostly concentrates on single elementary particles. If we combine such particles into a composite quantum system, the number of degrees of freedom of the combined system grows exponentially with the number of particles. This is a major difficulty when we try to describe the dynamics of such a system, since the computational resources required for this task also grow exponentially. In the context of quantum information processing, this difficulty becomes the main source of power: in some situations, information processors based in quantum mechanics can process information exponentially faster than classical systems. From the perspective of a physicist, one of the most interesting applications of this type of information processing is the simulation of quantum systems. We call a quantum information processor that simulates other quantum systems a quantum simulator.

This review discusses a specific type of quantum simulator, based on nuclear spin qubits, and using nuclear magnetic resonance for processing. We review the basics of quantum information processing by nuclear magnetic resonance (NMR) as well as the fundamentals of quantum simulation and describe some simple applications that can readily be realized by today’s quantum computers. In particular, we discuss the simulation of quantum phase transitions: the qualitative changes that the ground states of some quantum mechanical systems exhibit when some parameters in their Hamiltonians change through some critical points. As specific examples, we consider quantum phase transitions where the relevant ground states are entangled. Chains of spins coupled by Heisenberg interactions represent an ideal system for studying these effects: depending on the type and strength of interactions, the ground states can be product states or they can be maximally entangled states representing different types of entanglement.

Graphical abstract

Keywords

quantum simulation / quantum computation / quantum information / quantum phase transition / nuclear magnetic resonance (NMR)

1 Introduction

Storing, exchanging and processing of information has always been a defining aspect of human culture. With the rise of semiconductor technology and microelectronics, it became possible to encode virtually every type of information electronically. By now, we have become accustomed to having digitally encoded information available at our fingertips, with continuously increasing speed of communication and processing, and virtually limitless storage capacity. On the level of the material basis, this is the result of great advances in the capacity and speed of the computer hardware.

While these improvements have continued at a steady pace [1], it is now becoming clear that there are limits to this process that are starting to make it ever more difficult to continue this trend [2]. Some of these limits are technological in nature and a large number of ingenious engineers are working hard to find ways around these obstacles; examples are speed limits that may be overcome by adjustments to the manufacturing process, e.g., switching from silicon to silicon-germanium or using new types of insulator materials.

A second type of problems are fundamental physical problems, such as the limitations on the size of electronic circuits by the finite size of atoms or the amount of heat being dissipated by the electronic circuits [3]. These limitations are much more fundamental and cannot be overcome by simply switching the material basis of the electronic circuit technology. The third type of limitations are software: there are a number of computational problems where no efficient solution is known. Here, examples include optimization problems like the “traveling salesman problem” or the factorization of large numbers. Again, for these examples no efficient solutions are known, and it is not even known if they exist.

One possible way around both these fundamental problems was sketched out in 1982, at an international conference on the “Physics of Computation”, organized at the Massachusetts Institute of Technology by Ed Fredkin, Rolf Landauer, and Tom Toffoli. At this conference, several participants explored the fundamental limits that physics imposes on the speed of computation [4], such as the minimal amount of energy that computational machinery must dissipate. One of the participants was Richard Feynman from Caltech. In his lecture, he discussed the requirements for numerical simulations of quantum mechanical systems and showed that this is also a “hard” problem, i.e. one for which no algorithm is known that scales less than exponentially with the size of the system [5]. He suggested, however, that it might be possible to get around this obstacle by using quantum mechanical systems as the carriers of information and construct computing devices that operate according to Schrödinger’s equation. That quantum mechanical systems can really be used for storing and processing information was shown by Benioff [6] by constructing “Quantum mechanical Hamiltonian models of Turing machines”, which are universal computers. Since the time evolution of quantum mechanical systems is a unitary operation and therefore inherently reversible, quantum systems are a specific example of reversible computation, which demonstrates that computation can be performed without any energy dissipation. At that time, it was not known if these devices, now called “quantum computers” can actually be more powerful than classical computers. This conjecture was proven in 1993 by Bernstein and Vazirani [7], and shortly thereafter, the first algorithms were devised that operate exponentially faster on quantum computers than on classical computers [8, 9].

Following this idea, Lloyd [10] showed in 1996 that universal quantum simulators can be built from quantum mechanical systems. Subsequently, quantum simulation algorithms were developed for physical problems [11–17], such as the calculation of eigenstates and eigenvectors for given Hamiltonians [13], the energy spectrum [14], the spectrum of observables or relevant correlation functions [11] and the dynamics of many-body Fermi systems [18] as well as bosonic systems [12]. Besides these typical physical problems, quantum simulators should also be able to solve chemical problems. Some polynomial-time quantum algorithms have been proposed to perform the simulation for the calculations of the molecular energies [19, 20], the thermal rate constant of chemical reactions [21], chemical reactions using quantum dots [22] and chemical dynamics [23]. An attractive application of quantum simulations would be to calculate the folding of a protein under conditions of slow cooling.

Quantum simulators may turn out to be the field where quantum computers first successfully challenged classical computers. While algorithms such as Shors factoring algorithm require quantum computers with at least a few thousand qubits, it may be possible to perform useful quantum simulations with quantum computers that incorporate as few as 20–50 qubits. They could be used for a variety of problems in physics, chemistry and biology.

A universal quantum computer is particularly appealing, since it can efficiently simulate the dynamics of any quantum system. However, from the experimental perspective, it may be more efficient to first realize a special-purpose quantum simulator that has less stringent requirements, e.g., on the type of gate operations that can be applied. Some successful examples are quantum simulation of the Bose–Hubbard model using cold atoms in optical lattices [24, 25] and that of the Anderson model, the Kondo lattice and the multichannel Kondo model using superconducting qubits [26].

In comparison to this vast body of theoretical work, relatively little experimental work has been published. The first few examples used liquid-state NMR [27–31] to simulate, e.g., three spins that are coupled by a 3-particle interaction [27], the dynamics of truncated quantum harmonic and anharmonic oscillators [28], or the fermionic Fano–Anderson model for many-body problems [31]. In the simulations, the average Hamiltonian theory (AHT) [32] is widely adopted to design appropriate control sequences for certain Hamiltonians to be simulated. Building a small-scale, special-purpose NMR quantum simulator is an interesting and significant work. On the technological side, the developments of operational demonstrations of efficient control of dynamics of few-qubit quantum systems that will be employed for performing quantum simulations will provide the groundwork for large-scale universal quantum simulators based on other physical systems.

In this paper, we first give a brief introduction to NMR quantum information processing (QIP). Based on this, we describe and experimentally demonstrate the quantum simulation of quantum phase transitions (QPTs, an important physical phenomenon in condensed-matter physics) for certain spin chains with many-body interaction on few-qubit liquid-state NMR simulators. Incorporating adiabatic evolution into the Hamiltonian simulation allows one to drive the system through quantum phase transitions and study the ground state properties in the different parameter ranges. In particular, we show the relationship between the ground-state entanglement and QPTs in the critical systems.

The paper is organized as follows: In Section 2, we introduce two quantum computation models: the standard circuit model and adiabatic quantum computation (AQC) model, in which the principles of implementation of quantum computation by NMR QIP are explained. In Section 3, we summarize the general scheme of quantum simulation by NMR QIP. We then describe two methods to complete Hamiltonian simulation, which is the key step in quantum simulation in Section 4. As an important indicator for some quantum phase transitions (QPTs), we explain how to characterize entangled ground states and introduce the concept of quantum phase transitions in Section 5. In Sections 6 and 7, we illustrate these concepts using simple examples in two- and three-qubit systems.

2 Computing with quanta

2.1 Fundamentals

Storing and processing of information are the main tasks of computers. Conventionally, the information is decomposed into bits, i.e., binary digits, which can assume the values 0 or 1. These units of information can be stored in different physical systems that can assume 2 distinct states, e.g., a capacitor with different amounts of charge or the magnetization of a spot on a magnetic disc, which can be oriented in different directions. Processing this information implies that a string of bits is changed into a different string, according to a finite set of instructions called the algorithm. In the case of a quantum computer,the system that represents the information is a quantum mechanical system, and the logical operations correspond to the temporal evolution of this system [33, 34].

The main difference between the quantum mechanical and the classical computational models is that in classical computers, the information being processed always corresponds to a single string of bits. In quantum mechanical information processors, however, the information being processed can also be a superposition of all the possible strings of bits. The logical operations are then applied to this superposition. Due to the linearity of the Schrödinger equation, the superposition of all possible input states evolves into a superposition of all states that would be obtained by processing the initial states individually. Algorithms taking advantage of this possibility can therefore increase the speed of computation dramatically: in some cases, the speed increases exponentially with the number of qubits of the processor. This possibility of processing superpositions of basis states is often referred to as “quantum parallelism”.

Different computational models exist for converting a given input into the desired solution. The model that has been most thoroughly investigated so far is the so-called circuit or network model [33, 34], which is closely related to the processing of classical information.

Figure 1 illustrates the basics of the network model of quantum information processing [35, 36]. The quantum information is stored in a string of qubits called the quantum register. Before the computation can start, this register has to be initialized into an initial state containing the input and, possibly, additional qubits required by the algorithm. Logical operations are then implemented as unitary transformations acting on the qubits. After the last operation, the result has to be read out by measuring suitable observables of some (possibly all) qubits.

2.2 Quantum bits

To store a quantum bit, we need a quantum system (i.e., a physical system that can be described by the laws of quantum mechanics) with two distinct (i.e., orthogonal) states, which we label |0⟩ and |1⟩. The only physical system that has exactly two orthogonal states is a spin with an angular momentum of I=\mathrm{\hslash }/2. In the context of this review, we will consider nuclear spins with I=1/2 (in units of \mathrm{\hslash }), which are easy to manipulate and store the quantum information for a relatively long time, since they are quite well isolated from other degrees of freedom.

While we are in principle free to choose any 2 orthogonal states as the basis states for the computation, one usually opts for the states parallel and anti-parallel to the direction of an external magnetic field as the basis states and labels them

\begin{array}{}& |0⟩=\left(\begin{array}{c}1\\ 0\end{array}\right)\text{ and }|1⟩=\left(\begin{array}{c}0\\ 1\end{array}\right)\end{array}

Each qubit can be in an arbitrary superposition of these states,

\begin{array}{}& |\psi ⟩=\alpha |0⟩+\beta |1⟩;{|\alpha |}^2+{|\beta |}^2=1\end{array}

Like in the classical case, we need a string of qubits to store more information. For instance, a pure state for a string of N qubits is represented by

\begin{array}{}& |{\psi }_N⟩={\displaystyle \sum _{x=00\mathrm{\dots }0}^{11\mathrm{\dots }1}}{c}_x|x⟩\end{array}

where the c_x are probability amplitudes with

\begin{array}{}& {\displaystyle \sum _{x=00\mathrm{\dots }0}^{11\mathrm{\dots }1}}{|c_x|}^2=1\end{array}

The quantum information is then contained in a state that lives in the Hilbert space spanned by the 2^N computational basis states from |00\mathrm{\dots }0⟩ to |11\mathrm{\dots }1⟩.

If the qubits are realized by nuclear spins, it is usually advantageous to choose a spin-1/2 system. Suitable nuclei are ¹H, ¹³C, ¹⁹F, ¹⁵N, ³¹P and ²⁹Si. While a spin I>1/2 can store more (quantum-)information than a qubit, and it is possible to exploit this, most algorithms assume that qubits are used. Furthermore, spins I>1/2 have a nuclear quadrupole moment, which is a rather effective mechanism for coupling the spin to its environment and therefore dissipate the quantum information by relaxation.

2.3 Quantum gates for single qubits

Logical operations being applied to the qubits can be described in terms of operators acting on the states. For a single qubit, every operator can be written as a 2\times 2 complex matrix and expanded in terms of a basis of four orthogonal operators. Suitable operators are, e.g., the unit operator \mathbf{1} together with the three spin operators I_x, I_y and I_z:

\begin{array}{}\text{(1)}& \begin{array}{c}\mathbf{1}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),I_x={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\hfill \\ I_y={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}0& -\mathrm{i}\\ \mathrm{i}& 0\end{array}\right),I_z={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\hfill \end{array}\end{array}

An equivalent alternative consists of the set \mathbf{1}, I_{\pm }=I_x\pm I_y, and I_z.

A logical operation is then represented by a unitary operator

\begin{array}{}\text{(1)}& U={\mathrm{e}}^{-\mathrm{i}H\tau }\end{array}

where H is a Hermitian operator describing a physical interaction and τ the duration for which this interaction is active.

All logical operations on information encoded in quantum states can be decomposed into operations on single qubits (single-qubit operations) and operations acting on two qubits (two-qubit operations). The single-qubit operations correspond to rotations of the (pseudo-)spins. In the context of spin-qubits, it is convenient to distinguish between rotations around the z-axis, which corresponds to precession in the static magnetic field, and rotations around axes in the xy-plane, which are induced by radio-frequency pulses.

The z-rotations require the Zeeman interaction with a static magnetic field. This interaction is represented by the Hamiltonian

\begin{array}{}\text{(1)}& H_Z={\Omega }_z{I}_z\end{array}

where the Larmor frequency {\Omega }_z=-\gamma B_0is proportional to the strength of the magnetic field B₀ and the gyromagnetic ratio γ of the spin.

The rotations around axes in the xy-plane are generated by pulsed radio-frequency (RF) fields. These operations are best described in the so-called rotating frame, which rotates around the direction of the static magnetic field with a frequency equal to the radio-frequency of the applied field. In this reference frame, the interaction with the RF field becomes static and can be written as:

\begin{array}{}\text{(1)}& H_{\mathrm{RF}}={\omega }_1{I}_{\eta }\end{array}

where {\omega }_1 is the strength of the RF field and \eta is a direction in the xy-plane that can be chosen by the phase of the radio-frequency.

Using these interactions, it is possible to perform arbitrary logical operations on individual qubits. Every such operation corresponds to a rotation of the spins. Mathematically, it is described by a unitary transformation:

\begin{array}{}\text{(1)}& \begin{array}{cc}U(\delta ,\alpha ,\beta ,\theta )\hfill & \\ ={\mathrm{e}}^{\mathrm{i}\delta }\left(\begin{array}{cc}{\mathrm{e}}^{\mathrm{i}(\alpha +\beta )/2}\cos {\displaystyle \frac{\theta }{2}}& {\mathrm{e}}^{\mathrm{i}(\alpha -\beta )/2}\sin {\displaystyle \frac{\theta }{2}}\\ {\mathrm{e}}^{-\mathrm{i}(\alpha -\beta )/2}\sin {\displaystyle \frac{\theta }{2}}& {\mathrm{e}}^{-\mathrm{i}(\alpha +\beta )/2}\cos {\displaystyle \frac{\theta }{2}}\end{array}\right)\hfill & \\ ={\mathrm{e}}^{\mathrm{i}\delta }{R}_z(\alpha )R_y(\theta )R_z(\beta )\hfill & \end{array}\end{array}

Here, δ is an overall phase factor, while R_{\eta }(\alpha ) is a rotation by an angle α around the η axis. Rotations around axes in the xy-plane can be implemented by suitable radio-frequency pulses, where the rotation angle is given by the product of the nutation frequency {\omega }_1 and the pulse duration τ, \alpha ={\omega }_1\tau. Similarly, a rotation around the z-axis corresponds to an evolution in the static magnetic field for a suitable period τ.

2.4 Two-qubit operations

Any nontrivial algorithm includes conditional operations. In standard programming languages they are represented by operations like if, while or for. Clearly, such operations operate on more than one (qu)bit: at least one that controls the operation (the control qubit) and at least one on which the operation is performed (the target qubit). It is possible to show that all multi-qubit operations can be decomposed into one- and two-qubit operations [33, 34]. Therefore, we discuss here only two-qubit operations.

To implement controlled operations in a physical device, information must be exchanged between the qubits participating in the operation. This requires an interaction between them. In the case of spins, this interaction is a spin–spin coupling. In the simplest case (to which we restrict ourselves here), this interaction is represented by the Hamiltonian:

\begin{array}{}\text{(2)}& H_J=2J_{ij}{I}_z^i{I}_z^j\end{array}

The interaction can be turned into an operation by letting the system evolve under this Hamiltonian. In reality, the coupling operator (2) is only part of the complete system Hamiltonian. For two spins, a typical system Hamiltonian takes the form

\begin{array}{}\text{(3)}& H_{\mathrm{NMR}}={\Omega }_z^1{I}_z^1+{\Omega }_z^2{I}_z^2+H_J\end{array}

Evolution under this Hamiltonian will therefore result in a different logical operation than under the pure coupling operator. To turn this into the desired evolution

\begin{array}{}\text{(3)}& R_{zizj}(\varphi )={\mathrm{e}}^{-\mathrm{i}\varphi I_z^i{I}_z^j}\end{array}

we have to refocus the Zeeman-terms using π-pulses (rotations by an angle π around an axis in the xy-plane) in the middle and at the end of the evolution period. The π rotation of both spins around the x-axis corresponds to the operation

\begin{array}{}\text{(3)}& {\mathrm{e}}^{-\mathrm{i}\pi (I_x^i+I_x^j)}\end{array}

and the free evolution of duration τ to

\begin{array}{}\text{(3)}& {\mathrm{e}}^{-\mathrm{i}(2J_{ij}{I}_z^i{I}_z^j+{\Omega }_z^i{I}_z^i+{\Omega }_z^j{I}_z^j)\tau }\end{array}

The sequence of a free evolution period, a π rotation, another free evolution period and another π rotation results therefore in the following evolution:

\begin{array}{}\text{(4)}& \begin{array}{ccc}{\mathrm{e}}^{-\mathrm{i}\pi (I_x^i+I_x^j)}{\mathrm{e}}^{-\mathrm{i}\tau (2J_{ij}{I}_z^i{I}_z^j+{\Omega }_z^i{I}_z^i+{\Omega }_z^j{I}_z^j)}\hfill & & \\ \cdot {\mathrm{e}}^{-\mathrm{i}\pi (I_x^i+I_x^j)}{\mathrm{e}}^{-\mathrm{i}\tau (2J_{ij}{I}_z^i{I}_z^j+{\Omega }_z^i{I}_z^i+{\Omega }_z^j{I}_z^j)}\hfill & & \\ ={\mathrm{e}}^{-\mathrm{i}\tau (2J_{ij}{I}_z^i{I}_z^j-{\Omega }_z^i{I}_z^i-{\Omega }_z^j{I}_z^j)}{\mathrm{e}}^{-\mathrm{i}\tau (2J_{ij}{I}_z^i{I}_z^j+{\Omega }_z^i{I}_z^i+{\Omega }_z^j{I}_z^j)}\hfill & & \\ ={\mathrm{e}}^{-\mathrm{i4}{J}_{ij}\tau I_z^i{I}_z^j}\hfill & & \end{array}\end{array}

as required. Depending on the state of the system, the second π-pulse can be omitted.

In the following, we will use the graphical representation of Fig. 2 to represent the single-qubit operation R_{\eta }(\alpha ) and the two-qubit operation R_{zizj}(\varphi ) in a quantum network.

Interactions are important because their effect on a quantum state depends not only on the state of a single qubit, but on the state of several (usually 2) qubits. This is essential for implementing multiple qubit operations. Perhaps the most important multi-qubit operation is the controlled-NOT (CNOT) operation: a NOT operation on one qubit (the target qubit) that is only executed if the second qubit (the control qubit) is 1. If the control qubit is zero, the CNOT operation becomes a unity (=NOP) operation. In the 2-qubit basis \{|00⟩,|01⟩,|10⟩,|11⟩\}, this operation can be expressed in matrixform as:

\begin{array}{}\text{(4)}& \mathrm{CNOT}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)\end{array}

The first qubit is the control-bit: it is never changed by the CNOT operation. The target-bit is inverted if the control-bit has the value 1 and left unchanged if the control-bit is 0. The overall effect is therefore an exchange of the two states |10⟩ and |11⟩. This operation can be implemented by single-qubit operations and the two-qubit operation R_{zizj}(\varphi ) as:

\begin{array}{}\text{(5)}& \begin{array}{cc}\mathrm{CNOT}=\hfill & {\mathrm{e}}^{-\mathrm{i}\delta }{R}_{-y}^1\left({\displaystyle \frac{\pi }{2}}\right)R_{zizj}\left(\pi \right)R_{-x}^1\left({\displaystyle \frac{\pi }{2}}\right)\hfill \\ & \cdot R_{-z}^2\left({\displaystyle \frac{\pi }{2}}\right)R_{-y}^2\left({\displaystyle \frac{\pi }{2}}\right)\hfill \end{array}\end{array}

Each rotation R_{\eta }^k acts only on a single qubit, indicated by the superscript k. The subscript gives the rotation axis, the argument in parentheses is the rotation angle.

2.5 Multi-qubit quantum gate operations

Universal quantum gate operations on arbitrary numbers of qubits can be decomposed into one- and two-qubit operations. These operations must be applied selectively, i.e., only to those qubits (or pairs of qubits) whose operations are required by the algorithm. In NMR quantum computing, the selective addressing is implemented via resonant excitation. The resonance frequencies of the different spins must therefore be sufficiently different. This is ideally satisfied if the different qubits are represented by different spin species (e.g., ¹H, ¹³C, ¹⁵N etc.). However, in most cases, it is necessary to use more than one spin of the same spin species. In this case, the distinction must be accomplished by different chemical shifts.

The internal Hamiltonian of the spins (Zeeman + spin-spin coupling) is time-independent, while gate operations require the Hamiltonian to be switched on and off. This is straightforward for the single-qubit operations, but the two-qubit operations require us to switch the effect of couplings on and off. As we have discussed in the case of a two-qubit system, the two-qubit gate operation can be implemented by a free precession period with refocusing of the chemical shift. The same principle can be used for multi-qubit systems in which couplings are present between all quibts. To implement a two-qubit operation between qubits i and k only, we let the system evolve freely, but apply refocusing pulses to spins i and k. This refocuses their chemical shift as well as their couplings to all other spins, but leaves the coupling between i and k invariant. In most cases, it will be necessary to refocus also the couplings between the other spins and their chemical shifts. We refer to the literature for more details on how to implement these operations [37–40].

The scheme that we have discussed so far uses free precession R_{zizj} under the coupling operator to realize two-qubit operations. Alternatively, it is also possible to implement them by selective irradiation of individual transitions in coupled multi-qubit systems (see, e.g., Refs. [41–53]). As an example, a transition-selective π-pulse applied to the transition between the states |10⟩\leftrightarrow |11⟩ implements the operation [50]:

\begin{array}{}\text{(6)}& \left[{\pi }_y^{|10⟩\leftrightarrow |11⟩}\right]=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& -1\\ 0& 0& 1& 0\end{array}\right)\end{array}

which is closely related to the CNOT operation. “Transition-selective” means in this context that only the two states |10⟩ and |11⟩ may be affected by the pulse. This implies that the RF pulse must be long compared to the inverse of the coupling constant or, equivalently, that {\omega }_1\ll |J_{ij}|.

This technique of implementing multi-qubit operations by selective RF pulses can also be extended to three-qubit operations (and beyond). As an example, the operation

\begin{array}{}\text{(7)}& \left[{\pi }_y^{|110⟩\leftrightarrow |111⟩}\right]=\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -1\\ 0& 0& 0& 0& 0& 0& 1& 0\end{array}\right)\end{array}

which closely resembles the quantum TOFFOLI-gate (=CCNOT), has been implemented by transition-selective RF pulses [50].

2.6 Initialization

Before starting a calculation, the system must be initialized into a well-defined state. In many algorithms, it is assumed that the initial state is the state |00\mathrm{\dots }0⟩, which corresponds to the ground state of the spin system. Unfortunately, for nuclear spins in liquids, only a fraction of the spins is in this state. For an ensemble of spins I=1/2 in thermal equilibrium at a temperature T>1 K,Boltzmann’s law predicts that the states |\uparrow ⟩ and |\downarrow ⟩ are almost equally populated, with an occupation probability of p_i\approx 1/2\pm ϵ, where ϵ=\mathrm{\hslash }{\Omega }_z/(k_{\mathrm{B}}T)\sim {10}^{-5} at room temperature.

For an ensemble of systems, each consisting of n nuclear spins, the populations of the 2^n levels are

\begin{array}{}\text{(7)}& p_i\approx 2^{-n}\pm ϵ{\displaystyle \sum _{i=1}^n}{m}_i\end{array}

where m_i=\pm 1/2 are the magnetic quantum numbers of the n spins. This state is very close to the maximally mixed state, where ϵ=0. Nevertheless, it is possible to coax such a system into simulating systems that are in a pure state, including systems in the ground state. For this purpose, we use the fact that a system, in which initially all states are equally populated, does not contribute to an observable signal. Conceptually, this is easy to understand, since for every spin \uparrow, there is exactly one spin \downarrow, and the signal contributions from each such pair of spins cancel exactly. It is therefore not necessary to bring all spins into the ground state, but it is sufficient to have an excess of subsystems in the state |000\mathrm{\dots }0⟩, while all other states are equally populated. This system behaves exactly like a system with all spins in the ground state, except for a signal reduction by ε. This type of states are referred to as “pseudo-pure” or “effective pure” states. The procedure was originally proposed by Cory et al. [54] and by Gershenfeld and Chuang [55] in 1997.

The correct description of these states is in terms of a density operator

\begin{array}{}\text{(7)}& {\rho }_{pp}={\displaystyle \frac{1-ϵ}{2^n}}I+ϵ|\psi ⟩⟨\psi |\end{array}

where |\psi ⟩ describes a pure quantum state of n qubits.

These states cannot be prepared from the thermal equilibrium state by unitary operations. Instead, the preparation must include nonunitary processes, which can be generated by averaging over different unitary processes. Different averaging techniques have been proposed, including the temporal averaging [56–58], spatial averaging [45, 54, 59–61], logical labeling[44, 55, 56] spatially averaged logical labeling [45, 62, 63], using cat states [64], pairs of pseudopure states [65], and so on.

As an example, let us consider the preparation of a pseudo-pure state for a two-spin system by the spatial averaging method [54, 66]. If the first spin is a ¹³C and the second a ¹H nuclear spin, the thermal equilibrium state is given by the density operator:

\begin{array}{}\text{(7)}& {\rho }_{eq}\propto {\gamma }_C{I}_z^1+{\gamma }_H{I}_z^2\end{array}

where we have omitted the term proportional to the unity operator and {\gamma }_C and {\gamma }_H are, respectively, the gyromagnetic ratios for ¹³C and ¹H. Starting from this state, we first apply an rf pulse to the ¹H nuclear spin, thus converting polarization into transverse magnetization. The transverse magnetization can be destroyed by applying a pulsed magnetic field gradient, which dephases spins in different parts of the sample. The resulting state is thus

\begin{array}{}\text{(8)}& {\gamma }_C{I}_z^1+{\gamma }_H{I}_z^2\stackrel{R_y^2\left(\frac{\pi }{3}\right)-G_z}{⟶}{\gamma }_C(I_z^1+2I_z^2)\end{array}

where we have used {\gamma }_H\approx 4{\gamma }_C. By applying another set of two pulses with an intervening period of free precession, and a second gradient pulse, this state can be converted into

\begin{array}{}\text{(9)}& \stackrel{R_y^2\left(\frac{\pi }{4}\right)-R_{zizj}\left(\pi \right)-R_x^2\left(\frac{\pi }{4}\right)-G_z}{⟶}{\gamma }_C(I_z^1+I_z^2+2I_z^1{I}_z^2)\end{array}

In this state, the population of the ground state is 3{\gamma }_C, while all other states have the population -{\gamma }_C[163].

This procedure can be extended to larger systems. As an example, for a three-spin system consisting of ¹³C, ¹H and ¹⁹F, it is possible to generate a spatially averaged pseudopure state via the following sequence [54, 66]:

\begin{array}{}\text{(10)}& \begin{array}{ccc}\hfill {\gamma }_C{I}_z^1+{\gamma }_H{I}_z^2+{\gamma }_F{I}_z^3& & \\ \hfill \stackrel{R_y^3\left[a\cos (2{\gamma }_C/{\gamma }_F)\right]-G_z}{⟶}{\gamma }_C(I_z^1+4I_z^2+2I_z^3)& & \\ \hfill \stackrel{R_y^2\left(\frac{\pi }{2}\right)-R_{z2z3}\left(\pi \right)-R_x^2\left(\frac{\pi }{2}\right)}{⟶}{\gamma }_C(I_z^1+2I_z^2{I}_z^3+2I_z^3)& & \\ \hfill \stackrel{R_y^3\left(\frac{\pi }{4}\right)-R_{z1z3}\left(\pi \right)-R_x^3\left(\frac{\pi }{4}\right)-G_z^{\prime }}{⟶}{\gamma }_C[I_z^1+I_z^3+2I_z^1{I}_z^3& & \\ \hfill +4I_z^2(I_z^3+2I_z^1{I}_z^3)]& & \\ \hfill \stackrel{R_y^2\left(\frac{\pi }{4}\right)-R_{z2z3}\left(\pi \right)-R_x^2\left(\frac{\pi }{4}\right)-G_z^{\prime \prime }}{⟶}{\gamma }_C(I_z^1+I_z^2+I_z^3+& & \\ \hfill +2I_z^1{I}_z^2+2I_z^1{I}_z^3+2I_z^2{I}_z^3+4I_z^1{I}_z^2{I}_z^3)& & \end{array}\end{array}

The main difficulty with this approach is that the preparation eliminates a large fraction of the polarization, resulting in poor sensitivity [67]. Depending on the system, it is sometimes possible to counter this sensitivity decrease by techniques for polarization enhancement, such as optical pumping, dynamic nuclear polarization, the use of parahydrogen and algorithmic cooling (see Ref. [68] and references therein).

2.7 Readout

As the last step of the quantum computation, it is necessary to read out the result from the quantum state by performing an appropriate set of measurements. As NMR QIP uses an ensemble of nuclear spin systems, which are only weakly coupled to the detection apparatus, the readout differs from the idealized projective measurements that are usually considered in textbooks on quantum mechanics: in NMR, the quantum state does not collapse, but continues to evolve with very little perturbation. In addition, the results of the measurements are not populations of quantum states, but the expectation values of the transverse magnetization of the spins, I_x and I_y (which can be measured simultaneously).

For most algorithms, it is sufficient to measure the polarization of individual qubits independently. To improve the understanding of the relevant processes and to check individual computational steps, however, it is sometimes useful to determine the complete density operator. This can be achieved via quantum state tomography: converting those parts of the density operator that are not directly observable, into observable magnetization and subsequently measuring the NMR signal allows one to tomographically reconstruct the density operator [55]. In the case of a two-qubit system, quantum state tomography can be achieved with the help of four independent readout pulses [69], in the case of three qubits, the same result can be achieved with seven readout pulses [70]. For larger systems, the number of readout pulses required to reconstruct the complete density operator increases exponentially with the number of qubits, making this approach unsuitable for large quantum registers. Possible alternatives include a two-demensional Fourier transform technique [71], but again this approach is not scalable to very large systems.

2.8 Adiabatic quantum computing

While the “circuit” or “network” model has been used in most theoretical and experimental approaches to quantum information processing, it is sometimes useful to consider alternative approaches. The second approach that we consider here is the model of adiabatic quantum computing, which was proposed by Farhi et al. [72] and which has been shown to be equivalent (in terms of computational complexity) to the circuit model [73].

Adiabatic quantum computing differs radically from the network approach: the information is not “processed” in the conventional sense. The main step is the generation of a problem Hamiltonian H_P, whose ground state encodes the desired output. Possible examples include search algorithms, where the Hamiltonian represents a penalty function that vanishes in the ground state.

To find the solution, it is then necessary to bring the quantum register into the ground state of the problem Hamiltonian. This is obviously a nontrivial task since we do not know the state beforehand. Here, the adiabatic theorem of quantum mechanics comes to the rescue: if we prepare the system in the ground state of any suitable Hamiltonian and then slowly change the Hamiltonian until it becomes equal to the problem Hamiltonian, the quantum adiabatic theorem predicts that the system will remain in the ground state and therefore encodes the desired solution once the Hamiltonian has become equal to the problem Hamiltonian. The procedure starts at time t=0 with the system in the ground state {\psi }_g(0) of a Hamiltonian H₀. The Hamiltonian is then changed continuously into the problem Hamiltonian H_P over a time T:

\begin{array}{}\text{(11)}& H(t)=[1-s(t)]H_0+s(t)H_P\end{array}

where the function s(t) increases monotonically from 0 at t=0 to 1 at t=T. The solution of the optimization problem is then determined by measuring the final ground state {\psi }_g(T) of H_P.

For this procedure to succeed, three conditions must be met:

• The system is initially in the ground state of the initial Hamiltonian H(0).

• There is a gap between the two lowest energy eigenvalues of the Hamiltonian at all times, {\epsilon }_e(t)-{\epsilon }_g(t)>0.

• The Hamiltonian H(t) is varied sufficiently slowly, such that \left|⟨{\psi }_g(t)|{\dot{\psi }}_e(t)⟩/[{\epsilon }_e(t)-{\epsilon }_g(t)]\right|\ll 1. Here |{\psi }_g(t)⟩ and |{\psi }_e(t)⟩ refer to the instantaneous ground state and the first excited state, respectively, and {\epsilon }_g(t), {\epsilon }_e(t) are the corresponding energies.

Several adiabatic quantum algorithms have been developed for solving classical optimization problems. Examples include the 3-SAT algorithm [72, 74] and the unsorted databases search [75]. Some small-scale adiabatic algorithms have already been implemented experimentally, using NMR systems [76–78]. In comparison to the circuit model, AQC appears to offer some advantages:

• It eliminates the need to perform fast quantum logic operations and measurements, which is particularly troublesome in the context of fault-tolerant quantum computation.

• It is inherently robust against errors due to dephasing, environmental noise and some unitary control errors [79, 80] as long as the system is kept at low temperature.

On the other hand, the scaling properties of some adiabatic quantum algorithms, especially for nondeterministic polynomial (NP)-complete problems [81], are still not completely clear. Certain specialized examples have been analyzed numerically. In some cases, these tests indicated exponential complexity [82], in others polynomial [72, 83]. It has been believed that the worst case performance is exponential [84] for NP-complete problems.

3 Basics of quantum simulation

For the purpose of this review, we will concentrate on one particular application of quantum information processing, the simulation of other quantum systems on a quantum computer. This was one of the original motivations for looking into the possibility of using quantum systems for information processing [5] and may well be the first application where quantum computers become sufficiently powerful to successfully challenge classical computers.

3.1 Mapping the system

Simulating another quantum system on a quantum computer [27, 28] requires first of all a quantum register that has at least as many degrees of freedom as the system to be simulated. Furthermore, the time evolution of the system must be simulated on the target physical system. In most cases, it is not possible to implement the Hamiltonian of the simulated system directly. In this common case, one has to implement a time-dependent Hamiltonian in such a way that after a total evolution time T, the overall evolution under the time-dependent Hamiltonian becomes the same as it would have been under the real Hamiltonian H_S.

Figure 3 summarizes the requirements for a quantum simulation: the left-hand column shows the system whose evolution we wish to study by simulating it on a quantum computer (right-hand column). The simulated system |s⟩ is mapped into the physical system |p⟩=\varphi |s⟩ by the unitary operation ϕ.

3.2 Evolution

If the simulated system evolves under a (time-independent) Hamiltonian H_S, and we require the relation between the simulated and the physical system to be time-independent, i.e.,

\begin{array}{}\text{(11)}& |p(t)⟩=\varphi |s(t)⟩,0⩽t⩽T\end{array}

the propagator for the physical system must be

\begin{array}{}\text{(11)}& V_T=\varphi U{\varphi }^{-1}\end{array}

In most cases, it will not be possible to directly implement a time-independent Hamiltonian H_P such that the desired evolution is obtained directly, V_T={\mathrm{e}}^{-\mathrm{i}{H}_{P}T}. Like in other quantum computational problems, it is then necessary to instead implement the overall propagator by a time-dependent Hamiltonian. In those cases, where the Hamiltonian is piecewise constant, we can write

\begin{array}{}\text{(12)}& V_T={\displaystyle \prod _k}{\mathrm{e}}^{-\mathrm{i}{H}^{(k)}{\tau }_k}\end{array}

where {\tau }_k is the time period during which the physical system evolves under the Hamiltonian H^{(k)}. Writing also

\begin{array}{}\text{(12)}& V_T={\mathrm{e}}^{-\mathrm{i}{\overline{H}}_{P}T}\end{array}

we obtain an effective Hamiltonian {\overline{H}}_P, which formally generates the desired time-evolution and satisfies the simulation condition

\begin{array}{}\text{(12)}& {\overline{H}}_P=\varphi H_S{\varphi }^{-1}\end{array}

For our purpose, the physical system P corresponds to a nuclear spin system. Its internal Hamiltonian is

where {\Omega }_z^i represents the Larmor frequency of spin i, J_{ij} the strength of the coupling between spins i and j, n is the number of spins, and we have assumed weak coupling.

The external Hamiltonian describing the effect of the RF pulses is

\begin{array}{}\text{(14)}& H_{\mathrm{ext}}={\displaystyle \sum _{i=1}^n}{\mathrm{e}}^{-\mathrm{i}({\omega }_{\mathrm{RF}}^i{t}_i+{\phi }_i)I_z^i}(-{\omega }_1^i{I}_x^i){\mathrm{e}}^{\mathrm{i}({\omega }_{\mathrm{RF}}^i{t}_i+{\phi }_i)I_z^i}\end{array}

where {\omega }_{\mathrm{RF}}^i is the angular frequency, {\phi }_i the initial phase, t_i the duration, and {\omega }_1^i the amplitude of the rf pulse acting on spin i. This is the most general case, where each spin interacts with its own rf field. In practice, we often work with several qubits corresponding to the same spin species. In those cases, the fields acting on different qubits may be identical.

To implement a quantum simulation, we have to choose the correct sequence of phases, frequencies, amplitudes and pulse durations, such that the resulting sequence of Hamiltonians

\begin{array}{}\text{(14)}& H^{(k)}=H_{\mathrm{NMR}}+H_{\mathrm{ext}}^{(k)}\end{array}

fulfills Eq. (12). This condition can always be fulfilled if the individual qubits can be distinguished by their frequencies and they are part of a single coupling network.

4 Generating effective Hamiltonians

4.1 Mathematical basics

The evolution that generates the required unitary transformation V_T can be generated in different ways. Here, we will always assume that the complete Hamiltonian of the physical system P consists of an internal part H_{\mathrm{NMR}}, which is time-independent, and a control part H_{\mathrm{ext}}, which describes the interaction with experimentally adjustable control fields that will be taken to be piecewise constant.

Depending on the sequence of the control fields, the overall propagator can be generated in different ways. In simple cases, it is possible to calculate the propagator for each segment and multiply them to obtain the overall propagator. We will refer to this approach as “exact decomposition”. As a simple example, we generate a Hamiltonian

\begin{array}{}\text{(14)}& H_{XX}=d{I}_x^i{I}_x^k\end{array}

when the available Hamiltonian includes single qubit operations and

\begin{array}{}\text{(14)}& H_{ZZ}=d{I}_z^i{I}_z^k\end{array}

We may then convert this Hamiltonian into the desired one by “sandwiching” it between two {(\pi /2)}_y rotations:

\begin{array}{}\text{(14)}& {\mathrm{e}}^{-\mathrm{i}{H}_{XX}\tau }={\mathrm{e}}^{-\mathrm{i}\frac{\pi }{2}(I_y^i+I_y^k)}{\mathrm{e}}^{-\mathrm{i}{H}_{ZZ}\tau }{\mathrm{e}}^{\mathrm{i}\frac{\pi }{2}(I_y^i+I_y^k)}\end{array}

This approach provides a simple, intuitive and efficient way to generate the required Hamiltonian. However, when the system size increases, it becomes difficult and often impossible to find an analytical solution for the propagators in the individual period and the numerical approach (which would correspond to a classical simulation) scales poorly with system size. It may thus be necessary to use a perturbation expansion, where each period is short and the overall propagator can be calculated by a perturbation expansion in the duration of each pulse cycle.

The mathematical basis for this perturbation theory is the Baker–Campbell–Hausdorff formula:

\begin{array}{}\text{(14)}& \begin{array}{cc}{\mathrm{e}}^B{\mathrm{e}}^A=\hfill & \exp \{A+B+{\displaystyle \frac{1}{2}}[B,A]\hfill \\ & +{\displaystyle \frac{1}{12}}([B,[B,A]]+[[B,A],A])+\mathrm{\dots }\}\hfill \end{array}\end{array}

In the context of NMR, it is often evaluated with the help of the Magnus expansion [32, 85]:

\begin{array}{}\text{(14)}& \overline{H}(\mathrm{\Delta }t)={\overline{H}}^{(0)}+{\overline{H}}^{(1)}+{\overline{H}}^{(2)}+\mathrm{\dots }\end{array}

with

\begin{array}{}\text{(14)}& \begin{array}{cc}{\overline{H}}^{(0)}=\hfill & {\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle {\int }_0^{\mathrm{\Delta }t}}d{t}_{1}H(t_1)\hfill \\ {\overline{H}}^{(1)}=\hfill & -{\displaystyle \frac{i}{2\mathrm{\Delta }t}}{\displaystyle {\int }_0^{\mathrm{\Delta }t}}d{t}_2{\displaystyle {\int }_0^{t_2}}d{t}_1[H(t_2),H(t_1)]\hfill \\ {\overline{H}}^{(2)}=\hfill & -{\displaystyle \frac{1}{6\mathrm{\Delta }t}}{\displaystyle {\int }_0^{\mathrm{\Delta }t}}d{t}_3{\displaystyle {\int }_0^{t_3}}d{t}_2\hfill \\ \hfill & \cdot {\displaystyle {\int }_0^{t_2}}\mathrm{d}{t}_1\{[H(t_3),[H(t_2),H(t_1)]]\hfill \\ \hfill & +[[H(t_3),H(t_2)],H(t_1)]\}\hfill \end{array}\end{array}

In many cases, Trotter’s formula is also helpful:

\begin{array}{}\text{(14)}& {\mathrm{e}}^{(A+B)\mathrm{\Delta }t}={\mathrm{e}}^{A\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{B\mathrm{\Delta }t}{\mathrm{e}}^{A\frac{\mathrm{\Delta }t}{2}}+O(\mathrm{\Delta }{t}^3)\end{array}

This symmetrisation in time makes the expression accurate to second order in \mathrm{\Delta }t. It can be extended to Hamiltonians with more than two components,

\begin{array}{}\text{(14)}& H={\displaystyle \sum _{k=1}^L}{H}_k\end{array}

For the present purpose, it is sufficient to assume that the Hamiltonian is time-independent and that the decomposition is such that each term H_k can be realized. We may then expand the propagator in a symmetric manner as:

\begin{array}{}\text{(15)}& \begin{array}{cc}{\mathrm{e}}^{-\mathrm{i}H\mathrm{\Delta }t}=\hfill & [{\mathrm{e}}^{-\mathrm{i}{H}_1\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_2\frac{\mathrm{\Delta }t}{2}}\mathrm{\dots }{\mathrm{e}}^{-\mathrm{i}{H}_L\frac{\mathrm{\Delta }t}{2}}]\hfill \\ & \cdot [{\mathrm{e}}^{-\mathrm{i}{H}_L\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_{L-1}\frac{\mathrm{\Delta }t}{2}}\mathrm{\dots }{\mathrm{e}}^{-\mathrm{i}{H}_1\frac{\mathrm{\Delta }t}{2}}]+O(\mathrm{\Delta }{t}^3)\hfill \end{array}\end{array}

which is again exact to second order in \mathrm{\Delta }t. By repeating this period with short evolution times \mathrm{\Delta }t, it is then possible to generate arbitrary evolutions. Higher-order decompositions are also possible [86].

4.2 Heisenberg Hamiltonian for 2 spins

As a specific example, we discuss the simulation of two spins interacting with local magnetic fields and with each other through a Heisenberg interaction. The system Hamiltonian is

\begin{array}{}\text{(16)}& \begin{array}{ccc}{H}_{XYZ}=\hfill & {\Omega }_z^1{I}_z^1+{\Omega }_z^2{I}_z^2+\hfill & \\ & J_x{I}_x^1{I}_x^2+J_y{I}_y^1{I}_y^2+J_z{I}_z^1{I}_z^2\hfill & \end{array}\end{array}

where the {\Omega }_z^k are the strengths of the external magnetic fields (along the z axis) acting on qubit k, and J_{\nu }(\nu =x,y,z) the Heisenberg exchange constants. Here {\Omega }_z^k and J_{\nu } are in angular frequency units. J_{\nu }>0 and correspond to the antiferromagnetic and ferromagnetic cases, respectively. For arbitrary J_{\nu }, this is often called the anisotropic Heisenberg XYZ model. Some special cases are:

• XXX (or isotropic Heisenberg): J_x=J_y=J_z

• XXZ: J_x=J_y\ne J_z

• XY: J_z=0

• XX: J_x=J_y,J_z=0

• Heisenberg–Ising: J_x=J_y=0

The system Hamiltonian (13) of the NMR system is an example of a Heisenberg–Ising Hamiltonian. Here, we discuss how to use it for simulating the behavior of a system with the more general Heisenberg Hamiltonian (16) by applying a suitable sequence of RF pulses and free evolutions.

4.3 Exact decomposition

The Hamiltonian (16) can be diagonalized by the transformation:

\begin{array}{}\text{(16)}& H_{XYZ}=R{H}_{diag}{R}^{†}\end{array}

where

\begin{array}{}\text{(17)}& R=\left(\begin{array}{cccc}\cos {\displaystyle \frac{{\theta }_1}{2}}\hfill & & & -\sin {\displaystyle \frac{{\theta }_1}{2}}\hfill \\ & \cos {\displaystyle \frac{{\theta }_2}{2}}\hfill & -\sin {\displaystyle \frac{{\theta }_2}{2}}\hfill & \\ & \sin {\displaystyle \frac{{\theta }_2}{2}}\hfill & \cos {\displaystyle \frac{{\theta }_2}{2}}\hfill & \\ \sin {\displaystyle \frac{{\theta }_1}{2}}\hfill & & & \cos {\displaystyle \frac{{\theta }_1}{2}}\hfill \end{array}\right)\end{array}

and

\begin{array}{}\text{(18)}& \begin{array}{c}\cos {\displaystyle \frac{{\theta }_1}{2}}=\sqrt{{\displaystyle \frac{{\eta }_1+{\Omega }_z^1+{\Omega }_z^2}{2{\eta }_1}}}\hfill \\ \cos {\displaystyle \frac{{\theta }_2}{2}}=\sqrt{{\displaystyle \frac{{\eta }_2+{\Omega }_z^1-{\Omega }_z^2}{2{\eta }_2}}}\hfill \end{array}\end{array}

with

\begin{array}{}\text{(18)}& \begin{array}{c}{\eta }_1=\sqrt{{({\Omega }_z^1+{\Omega }_z^2)}^2+{\left({\displaystyle \frac{J_x-J_y}{2}}\right)}^2}\hfill \\ {\eta }_2=\sqrt{{({\Omega }_z^1-{\Omega }_z^2)}^2+{\left({\displaystyle \frac{J_x+J_y}{2}}\right)}^2}\hfill \end{array}\end{array}

In terms of spin operators, the transformation R is

\begin{array}{}\text{(18)}& R={\mathrm{e}}^{-\mathrm{i}({\theta }_1+{\theta }_2)I_y^1{I}_x^2}{\mathrm{e}}^{-\mathrm{i}({\theta }_1-{\theta }_2)I_x^1{I}_y^2}\end{array}

This can be generated from the natural spin-spin coupling Hamiltonian by rotating the two spin operators from the z to the x and y axes, i.e., using the two relations:

\begin{array}{}\text{(19)}& I_y^1{I}_x^2=R_{-x}^1\left({\displaystyle \frac{\pi }{2}}\right)R_y^2\left({\displaystyle \frac{\pi }{2}}\right)I_z^1{I}_z^2{R}_{-y}^2\left({\displaystyle \frac{\pi }{2}}\right)R_x^1\left({\displaystyle \frac{\pi }{2}}\right)\end{array}

and

\begin{array}{}\text{(20)}& I_x^1{I}_y^2=R_y^1\left({\displaystyle \frac{\pi }{2}}\right)R_{-x}^2\left({\displaystyle \frac{\pi }{2}}\right)I_z^1{I}_z^2{R}_{-y}^1\left({\displaystyle \frac{\pi }{2}}\right)R_x^2\left({\displaystyle \frac{\pi }{2}}\right)\end{array}

The diagonal Hamiltonian H_{diag} is a special case of an NMR Hamiltonian (13):

\begin{array}{}\text{(20)}& H_{diag}={\displaystyle \frac{1}{2}}({\eta }_1+{\eta }_2)I_z^1+{\displaystyle \frac{1}{2}}({\eta }_1-{\eta }_2)I_z^2+J_z{I}_z^1{I}_z^2\end{array}

An evolution {\mathrm{e}}^{-\mathrm{i}{H}_{diag}\tau } can thus be implemented as

\begin{array}{}\text{(21)}& {\mathrm{e}}^{-\mathrm{i}{H}_{diag}\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{\mathrm{NMR}}{\tau }_1}\end{array}

provided we choose the experimentally accessible parameters as:

\begin{array}{}\text{(21)}& {\tau }_1=\tau {\displaystyle \frac{2J_z}{J_{12}}},{\Omega }_L^1=({\eta }_1+{\eta }_2){\displaystyle \frac{\tau }{2{\tau }_1}},{\Omega }_L^2=({\eta }_1-{\eta }_2){\displaystyle \frac{\tau }{2{\tau }_1}}\end{array}

Combining Eqs. (19)–(21), we can thus realize the propagator

\begin{array}{}\text{(21)}& U_{XYZ}(\tau )={\mathrm{e}}^{-\mathrm{i}{H}_{XYZ}\tau }=R{\mathrm{e}}^{-\mathrm{i}{H}_{diag}\tau }{R}^{†}\end{array}

by the pulse sequence shown in Fig. 4.

4.4 Perturbation expansion

The same evolution U_{XYZ}(\tau ) can also be realized in a different way. For this purpose, we use the Trotter approximation to write the Heisenberg Hamiltonian as:

\begin{array}{}\text{(21)}& H_{XYZ}=H_z+H_{xy}\end{array}

where

\begin{array}{}\text{(21)}& H_z={\Omega }_z^1{I}_z^1+{\Omega }_z^2{I}_z^2+J_z{I}_z^1{I}_z^2\end{array}

and

\begin{array}{}\text{(21)}& H_{xy}=J_x{I}_x^1{I}_x^2+J_y{I}_y^1{I}_y^2\end{array}

The implementation of {\mathrm{e}}^{-\mathrm{i}{H}_z\mathrm{\Delta }t}was discussed above, and the evolution under the transverse term can be realized as:

\begin{array}{}\text{(22)}& \begin{array}{ccc}{\mathrm{e}}^{-\mathrm{i}{H}_{xy}\mathrm{\Delta }t}=\hfill & & \\ R_y^1\left({\displaystyle \frac{\pi }{2}}\right)R_y^2\left({\displaystyle \frac{\pi }{2}}\right){\mathrm{e}}^{-\mathrm{i}{J}_x{I}_z^1{I}_z^2\mathrm{\Delta }t}{R}_{-y}^1\left({\displaystyle \frac{\pi }{2}}\right)R_{-y}^2\left({\displaystyle \frac{\pi }{2}}\right)\hfill & & \\ \cdot R_{-x}^1\left({\displaystyle \frac{\pi }{2}}\right)R_{-x}^2\left({\displaystyle \frac{\pi }{2}}\right){\mathrm{e}}^{-\mathrm{i}{J}_y{I}_z^1{I}_z^2\mathrm{\Delta }t}{R}_x^1\left({\displaystyle \frac{\pi }{2}}\right)R_x^2\left({\displaystyle \frac{\pi }{2}}\right)\hfill & & \end{array}\end{array}

This decomposition is exact, since the two terms commute, [I_x^1{I}_x^2,I_y^1{I}_y^2]=0. The overall evolution under H_{XYZ}=H_z+H_{xy}, however, is approximate, since the two terms do not commute. Using the Trotter expansion, we can approximate the Heisenberg evolution operator by

\begin{array}{}\text{(23)}& U_{ap}={({\mathrm{e}}^{-\mathrm{i}{H}_z\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_{xy}\mathrm{\Delta }t}{\mathrm{e}}^{-\mathrm{i}{H}_z\frac{\mathrm{\Delta }t}{2}})}^M\end{array}

This approximates the true evolution up to second order in \mathrm{\Delta }t:

\begin{array}{}\text{(24)}& U_{XYZ}(\tau )=U_{ap}+O(\mathrm{\Delta }{t}^3)\end{array}

The segment duration \mathrm{\Delta }t determines the quality of the approximation; it has to be kept short in order to minimize the higher-order terms. If the required total evolution time τ is longer than a reasonable choice of \mathrm{\Delta }t allows, \tau >\mathrm{\Delta }t, we repeat the segment M times, such that \tau =M\mathrm{\Delta }t. Figure 5 shows the resulting sequence, using the same representation as before.

As a check of the effect of the duration of the segment, Figure 6 shows the fidelity of the resulting evolution as a function of the duration \mathrm{\Delta }t of each segment. For this analysis, we calculated the fidelity of the evolution operator U_{XYZ}(\tau ) as a function of the number of segments, M=1,\mathrm{\dots },30 while keeping the total duration constant at \tau =100. Since this fidelity depends on the parameters in the Hamiltonian, we averaged the fidelity over 200 different, randomly chosen parameter sets ({\Omega }_z^1,{\Omega }_z^2,J_x,J_y,J_z) with coupling constants in the range [0…1]. For each instance, the fidelity was calculated as:

\begin{array}{}\text{(24)}& F={|Tr(U_{XYZ}{(\tau )}^{†}{U}_{ap})/4|}^2\end{array}

5 Entanglement and quantum phase transitions

5.1 Basics

Nonlocality appears to be one of the fundamental but still baffling properties of quantum mechanics. Entangled states are quantum states that cannot be understood by the properties of their constituents alone. Among the best known examples is the singlet state of two identical particles, such as two electrons.

For a long time, entanglement was considered as a curiosity only of interest to those who were interested in the fundamentals of quantum mechanics. However, this has changed recently, with the advent of quantum technologies, where entanglement is now often regarded as a resource rather than a curiosity. This is most clearly the case in quantum communication, where entanglement can even be considered as a commodity, since it allows, e.g., to securely exchange data, eliminating many possibilities for eavesdropping.

This increased and more widespread interest has also initiated research into better characterization of entanglement, including the study of different forms of entanglement, and its quantification, leading to questions such as:

• How strongly are two given subsystems entangled?

• Which subsystems of a group are entangled with what other subsystems?

• What happens if one subsystem is removed from an entangled system?

Answers to these (and many similar) questions have an impact on future tools of communication and secrecy, and possibly also to a number of other “quantum technologies”, such as quantum metrology.

While much of this research is still in its infancy, small systems, consisting of 2–3 subsystems can now be classified completely in terms of their entanglement. Nuclear spin systems have been used very successfully to study the various entangled states in small systems, since they can be quite readily prepared using the techniques of pulsed NMR. We start the discussion by reviewing some measures of entanglement for bipartite and tripartite systems.

5.2 Bipartite systems

A bipartite system is called “entangled” if it cannot be written as a product state of the two subsystems. For pure states of the bipartite system, this can easily be checked by tracing over one of the two subsystems. If the system is a product state, the resulting density operator for the subsystem is a pure state,

\begin{array}{}\text{(24)}& {\rho }_A={\mathrm{Tr}}_B({\rho }_{AB})={({\rho }_A)}^2\end{array}

The most common measure of entanglement is the concurrence [87], which is related to the “entanglement of formation”. Let A and B be a pair of qubits, whose state is described by the density matrix {\rho }_{AB} of the pair, which may be pure or mixed. We define the spin-flipped density matrix as:

\begin{array}{}\text{(25)}& {\tilde{\rho }}_{AB}=(2I_y^1\otimes 2I_y^2){\rho }_{AB}^{\ast }(2I_y^1\otimes 2I_y^2)\end{array}

where ^∗ denotes the complex conjugation. Let {\lambda }_1,{\lambda }_2,{\lambda }_3 and {\lambda }_4 be the square roots of the eigenvalues of {\rho }_{AB}{\tilde{\rho }}_{AB} in decreasing order. Then the concurrence C for {\rho }_{AB} is

\begin{array}{}\text{(26)}& C({\rho }_{AB})=\max \{{\lambda }_1-{\lambda }_2-{\lambda }_3-{\lambda }_4,0\}\end{array}

A simpler expression is possible if we restrict ourselves to pure two-qubit states {|\psi ⟩}_{AB}. Then, the concurrence can be calculated as:

\begin{array}{}\text{(27)}& C({|\psi ⟩}_{AB})=|⟨{\psi }_{AB}|2I_y^1\otimes 2I_y^2|{\psi }_{AB}^{\ast }⟩|\end{array}

Its value ranges from 0 (no entanglement) to 1 (maximum entanglement). An example of maximally entangled states are the Bell states:

\begin{array}{}\text{(28)}& \begin{array}{c}|{\Phi }_{\pm }⟩={\displaystyle \frac{1}{\sqrt{2}}}(|00⟩\pm |11⟩)\hfill \\ |{\Psi }_{\pm }⟩={\displaystyle \frac{1}{\sqrt{2}}}(|01⟩\pm |10⟩)\hfill \end{array}\end{array}

“Maximally entangled” means that if we consider only one of the two subsystems (e.g. A) by tracing over the other subsystem, the resulting state is maximally mixed, i.e., its density operator is equal to the unity operator:

\begin{array}{}\text{(28)}& {\rho }_A={\mathrm{Tr}}_B({\rho }_{AB})={\displaystyle \frac{1}{2}}\mathbf{1}\end{array}

5.3 Tripartite systems

While entanglement in bipartite systems can be quantified by a single measure (the concurrence), larger systems can support different types of entanglement. In tripartite systems, we can have either pairwise entanglement or essential three-way entanglement [88], which cannot be reduced to pairwise entanglement. To keep the situation (relatively) simple, we discuss here only tripartite systems ABC that are in a pure state. For these systems, we can calculate the pairwise entanglement between qubits A and B by first tracing over the third qubit (C) and then using Eq. (26) to calculate the bipartite concurrence C_{AB}, and similarly for C_{AC} and C_{BC}.

In addition to this pairwise entanglement, tripartite systems also support an essential three-way entanglement [88]. This is found, e.g., in the GHZ states,

\begin{array}{}\text{(28)}& |{\mathrm{GHZ}}_{\pm }⟩={\displaystyle \frac{1}{\sqrt{2}}}(|000⟩\pm |111⟩)\end{array}

Since tracing over any of the three qubits leaves the remaining system in a maximally mixed state {\displaystyle \frac{1}{4}}\mathbf{1}, these states

• are maximally entangled;

• contain no pairwise entanglement, according to Eq. (26).

For pure states of a 3-qubit system ABC, it is possible to quantitate the different types of entanglement [88]. For the essential three-way entanglement, we use the 3-tangle

\begin{array}{}\text{(29)}& {\tau }_{ABC}=C_{A(BC)}^2-C_{AB}^2-C_{AC}^2\end{array}

The first term in Eq. (29) quantifies the entanglement of qubit A with the remaining subsystem BC. It is called the “generalized I-concurrence” [89, 90] and is defined as:

\begin{array}{}\text{(30)}& C_{A(BC)}=\sqrt{2\left[1-\mathrm{Tr}\left({\rho }_A^2\right)\right]}\end{array}

Here, the density operator {\rho }_A of the A subsystem is obtained by tracing over the degrees of freedom of B and C. If A is not entangled with the remaining subsystem, i.e., if the system is in a product state, {\rho }_A represents a pure state, \mathrm{Tr}\left({\rho }_A^2\right)=1, and C_{A(BC)}=0. For a maximally entangled state, {\rho }_A={\displaystyle \frac{1}{2}}\mathbf{1}, \mathrm{Tr}\left({\rho }_A^2\right)=1/2, and C_{A(BC)}=1, as expected. The remaining 2 terms in Eq. (29) quantify the pairwise entanglement of A with B and A with C. The 3-tangle can thus be understood as the difference between the entanglement of one of the qubits with the remaining system and the total pairwise entanglement of the same qubit with each of the other qubits.

As noted above, the GHZ states are examples (actually: the only examples in a tripartite system) for maximum essential three-way entanglement, {\tau }_{ABC}=1, while all measures of pairwise entanglement vanish, C_{\alpha \beta }=0,(\alpha ,\beta =A,B or C). In contrast to this, the W states,

\begin{array}{}\text{(31)}& \begin{array}{c}|{\mathrm{W}}_{001}⟩={\displaystyle \frac{1}{\sqrt{3}}}\left(|001⟩+|010⟩+|100⟩\right)\\ |{\mathrm{W}}_{110}⟩={\displaystyle \frac{1}{\sqrt{3}}}\left(|110⟩+|101⟩+|011⟩\right)\end{array}\end{array}

have no essential three-way entanglement {\tau }_{ABC}=0, but are strongly pairwise entangled, C_{AB}^2+C_{AC}^2+C_{BC}^2=4/3.

GHZ states and W states are the only two inequivalent kinds of genuine tripartite entanglement in a 3-qubit system [91]. Unfortunately, there is currently no known analytical expression for the 3-tangle of a mixed state of three qubits, nor is it known whether expressions analogous to Eq. (29) can be found for systems consisting of more than three qubits.

Another important method for characterizing entanglement uses specifically designed observables to distinguish between entangled and separable states in a quantum system. These observables are called entanglement witnesses [92–95], which can always be used to detect various forms of multipartite entanglement, provided some a priori knowledge is available about the states under investigation. In our context, the two witness operators

\begin{array}{}\text{(32)}& {\mathcal{W}}_{\mathrm{GHZ}}={\displaystyle \frac{3}{4}}\mathbf{1}-|{\mathrm{GHZ}}_{-}⟩⟨{\mathrm{GHZ}}_{-}|\end{array}

and

\begin{array}{}\text{(33)}& {\mathcal{W}}_{\mathrm{W}}={\displaystyle \frac{2}{3}}\mathbf{1}-|{\mathrm{W}}_{001}⟩⟨{\mathrm{W}}_{001}|\end{array}

are useful for distinguishing the two types of entanglement. If their expectation value is negative, they indicate genuine tripartite entanglement; further detects GHZ-type entanglement. These witnesses do not detect all types of GHZ- or W-entanglement, e.g., they do not detect the {\mathrm{W}}_{110} state. However, this state can be converted into the {\mathrm{W}}_{001} state by local unitary transformations; in the context of NMR, this means that such a state could be transformed by rf pulses and then detected by the same entanglement witness.

5.4 Quantum phase transitions

Phase transitions are a common feature in classical physics. Familiar examples are boiling and freezing of water. In these classical phase transitions, several properties of the system change discontinuously as a control parameter (e.g., temperature or pressure) changes by a small amount.

In quantum mechanical systems, we do not consider phase transitions of ensembles, but instead abrupt changes of the ground state of a Hamiltonian when a parameter of the Hamiltonian changes continuously - so-called quantum phase transitions (QPTs) [96]. Well-known examples are the superconductor-insulator transition and the paramagnetic-antiferromagnetic transition in quantum magnets. While classical phase transitions often occur as a consequence of temperature change, quantum phase transitions occur at T=0.

Besides its own remarkable properties, the investigation of such phase transitions is often tightly associated with the design, discovery and growth of novel materials [97]. Experimental observations of QPTs were reported only recently, e.g., in magnetic systems [98], heavy-fermion metals [99], common metals [100] and Bose–Einstein condensation [24]. A collection of reviews [97, 101–105] on QPTs reported the current status and recent developments in this field.

At the critical point where QPT occurs, the ground state of the system undergoes a qualitative change in some of its properties [96]. In this context, we will emphasize the change in entanglement, which is associated with various types of QPTs. Examples that were discussed in detail include the Heisenberg model [106–109] and the Hubbard model [110]. It is believed that the ground-state entanglement also plays a crucial role in other QPTs, like the change of conductivity in the Mott insulator-superfluid transition [111] and the quantum Hall effect [112]. Many of the relevant features, like the transition from a simple product state to a strongly entangled state, occur over a wide range of parameters and persist for infinite systems as well as for systems with as few as two spins [108, 113–115].

In most systems studied, attention was focused on two-body interactions, which are most readily accessible experimentally. On the other hand, spin systems with three-body interactions have been shown to exhibit exotic ground-state properties [116–127], which are unknown in systems with only two-body interactions. In cluster Hamiltonians, the presence of three-spin interactions results in unique entanglement properties and new critical phenomena that cannot be detected by two-point correlation functions, a traditional approach to measure criticality of ground states [122].

6 A simple example

6.1 Motivation

The goal of a quantum simulation is to observe the behavior of a system that does not exist or cannot be observed directly, or to examine properties that are not directly measurable in the real system. In our example, we use a system of nuclear spins to investigate entanglement and quantum phase transitions. The systems that we discussed in the previous section are not easy to find in nature and even if they exist, it is usually not straightforward to change the relevant parameters nor to observe their behavior as a function of these parameters. We demonstrate these issues here, first using the simplest possible example. In the following section, we move to a slightly larger system, where we can discuss the issues that arise when the system becomes more complex.

The simulation experiments that we discuss here consist of three steps:

(1) Initialize the quantum simulator into the ground state of the initial Hamiltonian H(0).

(2) Generate an effective Hamiltonian that corresponds to the target Hamiltonian and drive it through the relevant trajectory, including the quantum critical points, by adiabatically changing a suitable control parameter.

(3) Detect the QPTs using a suitable indicator.

6.2 Target system Hamiltonian

The simplest system that exhibits a QPT of ground-state entanglement consists of two spins in an external magetic field, coupled by an Ising-type interaction. The corresponding Hamiltonian is

\begin{array}{}\text{(34)}& H_{2\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}={\Omega }_z(I_z^1+I_z^2)+2J_z{I}_z^1{I}_z^2,J_z>0\end{array}

where {\Omega }_z is the magnetic field strength (in frequency units), and J_z is the spin-spin coupling constant.

A symmetry-adapted basis that is an eigenbasis for this Hamiltonian consists of the states \{|\uparrow \uparrow ⟩,|{\Psi }^{+}⟩,|\downarrow \downarrow ⟩,|{\Psi }^{-}⟩\}, with |{\Psi }^{\pm }⟩={\displaystyle \frac{1}{\sqrt{2}}}(|\uparrow \downarrow ⟩\pm |\downarrow \uparrow ⟩) and |\downarrow ⟩ and |\uparrow ⟩ the spin-down (m=-1/2) and spin up (m=+1/2) states. While the states |\uparrow \uparrow ⟩ and |\downarrow \downarrow ⟩are product states, the states |{\Psi }^{\pm }⟩ are maximally entangled states. |{\Psi }^{-}⟩ is the singlet state, the others are triplet states.

Since the ground state of this system is always one of the triplet states, and transitions to the singlet state are symmetry-forbidden, we can reduce our system of interest to the triplet states. The ground state |{\psi }_1⟩ depends then only on the dimensionless field strength g_z={\Omega }_z/J_z:

The ground state of the system is therefore a maximally entangled state if the field is small, , but a product state in the large field region, |g_z|>1. g_z=\pm 1 are quantum critical points, where the ground state changes from the ferromagnetically ordered high field states to the entangled, antiferromagnetic low-field state. Figure 7 summarizes this behavior.

To observe the transitions between these different ground states, we will drive the system adiabatically through the critical points. Such an adiabatic passage is only possible if the ground state is never degenerate. We therefore must avoid the degeneracies of the Hamiltonian (34). This is possible by adding a small transverse field to the system:

\begin{array}{}\text{(36)}& H_{\mathrm{Ising}}={\Omega }_z(I_z^1+I_z^2)+2J_z{I}_z^1{I}_z^2+{\Omega }_x(I_x^1+I_x^2)\end{array}

The transverse field is always kept small, |{\Omega }_x|\ll \max \{|{\Omega }_z|,|J_z|\}. As shown in Fig. 7, this modified Hamiltonian is nondegenerate and the QPTs now have a finite width of order {\Omega }_x/J_z.

6.3 Effective Hamiltonian

To generate the effective Hamiltonian, we used a pulse sequence whose parameters could be varied in such a way that the corresponding parameter in the effective Hamiltonian went through the desired trajectory adiabatically, i.e., such that the system always remained in the instantaneous ground state. We used the approximate expansion technique to create the target Hamiltonian (see Section 4.2):

\begin{array}{}\text{(36)}& \begin{array}{cc}{U}_{\mathrm{Ising}}(\mathrm{\Delta }t)\hfill & ={\mathrm{e}}^{-\mathrm{i}{H}_{\mathrm{Ising}}\mathrm{\Delta }t}\hfill \\ & \approx {\mathrm{e}}^{-\mathrm{i}{H}_x^{\mathrm{Ising}}\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_z^{\mathrm{Ising}}\mathrm{\Delta }t}{\mathrm{e}}^{-\mathrm{i}{H}_x^{\mathrm{Ising}}\frac{\mathrm{\Delta }t}{2}}\hfill \end{array}\end{array}

with

\begin{array}{}\text{(36)}& \begin{array}{ccc}{H}_x^{\mathrm{Ising}}={\Omega }_x(I_x^1+I_x^2)\hfill & & \\ H_z^{\mathrm{Ising}}={\Omega }_z(I_z^1+I_z^2)+2J_z{I}_z^1{I}_z^2\hfill & & \end{array}\end{array}

The approximation is valid as long as the evolution period \mathrm{\Delta }t is kept sufficiently short. The decomposed operator can be realized by incorporating the free evolutions under the natural NMR Hamiltonian H_{\mathrm{NMR}} with the RF pulses.

For the experimental implementation of the adiabatic passage, we discretized the time-dependent Hamiltonian H(t) into M+1 segments with peacewise constant Hamiltonians

\begin{array}{}\text{(36)}& H(m)=H\left[C\left({\displaystyle \frac{m}{M}}T\right)\right],m=0,1,\mathrm{\dots },M\end{array}

where T is the total duration of the adiabatic passage and C is the control parameter in H(t). The adiabatic condition is achieved when both T,M\to \mathrm{\infty } and the duration of each step \mathrm{\Delta }t\to 0. Figure 8 shows the sequence of RF pulses that generates the adiabatic evolution.

The experimental implementation generates an effective Hamiltonian that is constant for a time \mathrm{\Delta }t. This stepwise approximation is faithful if the duration of each time step is chosen such that (1) the time is short enough that the average Hamiltonian approximation holds and (2) the transfer remains adiabatic. While this calls for many short steps, there is also a lower limit for the duration of each step, which is dictated by experimental aspects: pulse errors and switching transients tend to generate errors that increase with the number of cycles.

We used a numerical optimisation procedure to determine the optimal sequence of Hamiltonians, taking the full level structure into account. Choosing a hyperbolic sine as the functional form, we optimised its parameters. Furthermore, to determine the optimal number of steps, we used the same numerical simulation, keeping the functional dependence {\Omega }_z vs. t constant, but increasing the number of steps M.

The results are summarized in Fig. 9, which plots the lowest fidelity encountered during each scan against the number of steps taken in the simulation. The fidelity is calculated as the overlap of the state with the ground state at the relevant position. The simulation shows also the effect of decoherence, which reduces the achievable fidelity if the total duration of the scan becomes comparable to the decoherence time. The number of steps for the discretized scans in the experiments was chosen near the top of the curves. The model that we used to take the effect of decoherence into account is similar to that of Vandersypen et al. [128].

6.4 Experimental results

The experiments were carried out at room temperature using a 500 MHz (11.7 Tesla) NMR spectrometer from Bruker Biospin. For the two-qubit quantum simulations, we chose the ¹³C and ¹H spins of ¹³C-labeled chloroform. The relatively large spin-spin coupling constant of J_{12}/\pi =214.94 Hz makes this molecule well suited for this simulation. The pseudo-pure initial state {\rho }_{pp}(\uparrow \uparrow ) was generated with the pulse sequence (9) described in Section 2. The fidelity of this state preparation was checked by quantum state tomography and found to be better than 0.99.

The adiabatic scan through the QPT was performed by shifting the RF frequencies of both channels after each period. In the following, we present the results as a function of the dimensionless magnetic field strength g_z={\Omega }_z/J_z, i.e., normalized to the strength of the coupling. Using the sweep g_z(t) shown in Fig. 9, the offset was changed from g_z=-3 to g_z=+3 in 60 steps. The evolution of the system during the scan was checked by performing a complete quantum state tomography after every second step during the experiment.

We calculated the concurrence C from the tomographically reconstructed density matrix according to Eq. (26). Figure 10 shows the measured concurrence C_{\exp }=C({\rho }_{\exp }) as individual points and compares them with the theoretical values C_{\mathrm{th}}. Both data sets clearly show the expected QPTs near the critical points g_z=\pm 1. The entangled ground state for is characterised by a concurrence close to 1, while the high field states are only weakly entangled (the entanglement vanishes for g_x=0).

The experimentally determined concurrence C remains below \sim 0.75, significantly less than the theoretical values. To verify that this deviation is due to decoherence, we simulated the experiment on a classical computer, taking into account the details of the pulse sequence as well as the effect of decoherence. We obtained good agreement between theoretical and experimental data if we assumed a total decoherence time of 130 ms, which is comparable to the experimental scan time of \sim 110 ms. Figure 10 shows the simulated values of the concurrence as triangles; their evolution during the scans is quite similar to that of the experimental data points.

To assess the quality of the adiabatic scan, we also determined the fidelities F_{\exp }=F({\rho }_{\exp })=\mathrm{Tr}({\rho }_{\mathrm{id}}{\rho }_{\exp }) from the tomographically reconstructed density operators, where {\rho }_{\mathrm{id}} is the ideal ground state at the measured time point. The fidelities, which are shown at the top of Fig. 10, deviate from unity when the system passes through the critical points and show some overall decrease due to decoherence. It can also been seen that the biggest deviation of the fidelity occurs at the critical points, i.e., the sensitivity of the fidelity can also give some information on critical points. Again, the simulated fidelities agree remarkably well with the experimental values.

As a second-order parameter, we also determined the two-spin correlation [96] ⟨4I_z^1{I}_z^2⟩=4\mathrm{T}\mathrm{r}({\rho }_{\exp }{I}_z^1{I}_z^2). Figure 10 shows the results. As expected, the system is ferromagnetically ordered (4⟨I_z^1{I}_z^2⟩=+1) at high fields, but turns to an antiferromagnetic state (4⟨I_z^1{I}_z^2⟩=-1) at low fields between the two quantum critical points. Comparing these two plots in Fig. 10, we find that the behavior of the concurrence is similar to the two-spin correlation.

The quantum state tomography used in these experiments provides very detailed information on the quantum system and contains all the possible information about its state. However, the number of measurements required for this analysis, as well as the date processing increases exponentially with the size of the system. It is therefore important to look into alternative methods for characterizing quantum systems and quantum phase transitions. As discussed elsewhere [129], one possible technique is to couple the system to one additional qubit. If the coupling is chosen well, a single measurement on this “ancilla qubit” is sufficient to detect the quantum phase transition, independent of the size of the system.

7 Three-qubit system

7.1 Three-spin chain

If we extend our system to three qubits, we can study several qualitatively new effects. In particular, we are interested in the different types of entanglement that we find in this system: Some states correspond to pairwise entanglement, while other states correspond to genuine tripartite entanglement. We are again interested in the quantum mechanical ground state of the system and investigate to which class of states it belongs, depending on the relevant parameters of the Hamiltonian. We consider a triangle of three spins 1/2 in a uniform magnetic field, interacting by both Ising-type two-body and three-body couplings:

\begin{array}{}\text{(37)}& \begin{array}{cc}{H}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=\hfill & {\Omega }_z(I_z^1+I_z^2+I_z^3)+{\Omega }_x(I_x^1+I_x^2+I_x^3)\hfill \\ & +2J_2(I_z^1{I}_z^2+I_z^2{I}_z^3+I_z^1{I}_z^3)+4J_3{I}_z^1{I}_z^2{I}_z^3\hfill \end{array}\end{array}

where {\Omega }_z and {\Omega }_x characterize the strengths of the longitudinal and transverse magnetic fields, and J₂, J₃ represent the two-body and three-body coupling constants. The transverse field is small, |{\Omega }_x|\ll 1.

Using perturbation theory in a symmetry-adapted basis \{|000⟩,|{\mathrm{W}}_{001}⟩,|{\mathrm{W}}_{110}⟩,|111⟩\}, we find the ground state of the Hamiltonian (37) as a function of the three control parameters {\Omega }_z,J_2, and J₃. Table 1 summarizes the different phases for the half-space {\Omega }_z⩾0. Similar expressions are obtained for . The competition between the three kinds of interactions (one-, two- and three-body) results in distinct types of ground states (product-, W- and GHZ-states) in different parameter ranges.

In the two-qubit system that we discussed in the preceding subsection, only one type of entanglement exists. In the three-spin systems, we obtain a qualitatively different behavior: In addition to the familiar bipartite entanglement, the system can also exhibit genuine tripartite entanglement. In particular, we will consider two classes of states that exhibit two types of entanglement—W states and GHZ states. Compared to a two-qubit system, a three-qubit system therefore has a much richer phase diagram and requires more specific measures for characterizing the different possible QPTs.

7.2 Target Hamiltonians and phases

To observe the different types of phase transitions in the Hamiltonian (16), we chose two different Hamiltonians:

\begin{array}{}\text{(38)}& \begin{array}{cc}{H}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}=\hfill & {\Omega }_z(I_z^1+I_z^2+I_z^3)+{\Omega }_x(I_x^1+I_x^2+I_x^3)\hfill \\ & +2J_2(I_z^1{I}_z^2+I_z^2{I}_z^3+I_z^1{I}_z^3)\hfill \end{array}\end{array}

with vanishing three-spin interaction but nonvanishing two-spin interactions and

\begin{array}{}\text{(39)}& H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}={\Omega }_x(I_x^1+I_x^2+I_x^3)+4J_3{I}_z^1{I}_z^2{I}_z^3\end{array}

with only three-spin interaction.

From Table 1 and Fig. 11, we find that the ground state of the first system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)} has a transition from a product state to a W entangled state at J_2=\pm {\Omega }_z/2, while the ground state of the second system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)} has a transition from a product state to a GHZ entangled state at J_3={\Omega }_x.

7.3 Physical system and effective Hamiltonians

Our three-qubit NMR quantum simulator is based on the molecule Diethyl-fluoromalonate, dissolved in acetone-d₆. Its structure is shown in Figure 12 and we use the ¹³C, ¹H, and ¹⁹F nuclear spins in the dashed oval as qubits. Figure 12 also shows the relevant parameters of this system.

The target Hamiltonians were created as an average Hamiltonian by concatenating evolutions with short periods:

\begin{array}{}\text{(39)}& \begin{array}{cc}{U}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}(\mathrm{\Delta }t)\hfill & ={\mathrm{e}}^{-\mathrm{i}{H}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}\mathrm{\Delta }t}\hfill \\ & \approx {\mathrm{e}}^{-\mathrm{i}{H}_x^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_z^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\mathrm{\Delta }t}\hfill \\ & \cdot {\mathrm{e}}^{-\mathrm{i}{H}_{zz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\mathrm{\Delta }t}{\mathrm{e}}^{-\mathrm{i}{H}_x^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\frac{\mathrm{\Delta }t}{2}}\hfill \end{array}\end{array}

and

\begin{array}{}\text{(39)}& \begin{array}{cc}{U}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}(\mathrm{\Delta }t)\hfill & ={\mathrm{e}}^{-\mathrm{i}{H}_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}\mathrm{\Delta }t}\hfill \\ & \approx {\mathrm{e}}^{-\mathrm{i}{H}_x^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\frac{\mathrm{\Delta }t}{2}}{\mathrm{e}}^{-\mathrm{i}{H}_{zzz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\mathrm{\Delta }t}{\mathrm{e}}^{-\mathrm{i}{H}_x^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\frac{\mathrm{\Delta }t}{2}}\hfill \end{array}\end{array}

Here,

\begin{array}{}\text{(39)}& \begin{array}{cc}{H}_x^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}={\Omega }_x(I_x^1+I_x^2+I_x^3)\hfill & \\ H_z^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}={\Omega }_z(I_z^1+I_z^2+I_z^3)\hfill & \\ H_{zz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=2J_2(I_z^1{I}_z^2+I_z^2{I}_z^3+I_z^1{I}_z^3)\hfill & \\ H_{zzz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=4J_3{I}_z^1{I}_z^2{I}_z^3\hfill & \end{array}\end{array}

All the pairwise couplings between the spins in the molecule have nonzero strength. However, they are not identical and therefore cannot be used directly. We used a sutiable refocusing scheme [39,40,130] to create the Hamiltonian H_{zz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}} with equal two-spin interactions. Figure 13 (a) shows the pulse sequence used in the experiment. For all three spins, nonzero offsets FQn (n=1,2,3) were used; their values are given in the figure caption.

Three-spin interactions do not occur naturally in NMR, but they can be simulated by a combination of two-spin interactions and RF pulses [27]:

\begin{array}{}\text{(40)}& \begin{array}{cc}{\mathrm{e}}^{-\mathrm{i}{H}_{zzz}^{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\mathrm{\Delta }t}=\hfill & R_x^1\left({\displaystyle \frac{\pi }{2}}\right){\mathrm{e}}^{-\mathrm{i}{I}_z^1{I}_z^3\pi /4}{R}_y^1\left({\displaystyle \frac{\pi }{2}}\right){\mathrm{e}}^{-\mathrm{i2}{J}_3{I}_z^1{I}_z^2\mathrm{\Delta }t}\hfill \\ \hfill & \cdot R_{-y}^1\left({\displaystyle \frac{\pi }{2}}\right){\mathrm{e}}^{\mathrm{i}{I}_z^1{I}_z^3\pi /4}{R}_{-x}^1\left({\displaystyle \frac{\pi }{2}}\right)\hfill \end{array}\end{array}

To implement this, we also use the relation

\begin{array}{}\text{(40)}& {\mathrm{e}}^{\mathrm{i}{I}_z^i{I}_z^j\pi /4}=I_x^i{\mathrm{e}}^{-\mathrm{i}{I}_z^i{I}_z^j\pi /4}{I}_x^i\end{array}

The pulse sequence that generates the Hamiltonian H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)} is shown in Fig. 13 (b).

As in the two-spin system, we used a hyperbolic sine for the variation of the control parameters J₂ and J₃. To optimize the experimental parameters, we first performed a numerical simulation. The results are shown in Fig. 14. For the experiments, we chose the number of steps such that the fidelities were maximized, to M=8. The first adiabatic trajectory that we implemented corresponded to J_2=0\mathrm{\dots }2 in Hamiltonian H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}, with the duration of each segment set to \mathrm{\Delta }t=\pi /2. The longitudinal field in H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)} was set to {\Omega }_z=-2, and the transverse field to {\Omega }_x=0.09. The second trajectory implemented H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}, varying the 3-body interaction from J_3=0\mathrm{\dots }2. The other parameters were \mathrm{\Delta }t=\pi /2 and {\Omega }_x=0.12.

7.4 Experimental results

At the start of the simulation, the system was prepared in a pseudo-pure initial state {\rho }_{pp}(\uparrow \uparrow \uparrow ) with an experimental fidelity of 0.99 by using the pulse sequence (10) (see Section 2.1). The density matrix reconstructed from the experimental results is shown in Fig. 15 (a).

We first checked the final states generated at the end of the adiabatic scan by performing complete quantum state tomography [131]. The experimentally reconstructed density matrices are shown in Fig. 15, together with the initial state. They agree quite closely with the theoretical ones. To quantify the agreement between the theoretical and experimental states, we used two different measures to compare them: First, we considered the attenuated correlation [132]:

\begin{array}{}\text{(41)}& c({\rho }_{\exp })={\displaystyle \frac{\text{Tr}({\rho }_{\mathrm{th}}{\rho }_{\exp })}{\text{Tr}({\rho }_{\mathrm{th}}{\rho }_{\mathrm{th}})}}\end{array}

where {\rho }_{\mathrm{th}} is a theoretical output state, obtained by applying the ideal transformations to the experimental input state {\rho }_{000}^{\exp }, and {\rho }_{\exp } is the experimental output state obtained by density operator tomography. Eq. (41) implies that -1⩽c⩽1. The experimental values of the correlation for the two cases considered here were c({\rho }_{\exp }^{\mathrm{W}})=0.60 and c({\rho }_{\exp }^{\mathrm{GHZ}})=0.73. As the second measure, we used the definition of the state fidelity (42):

\begin{array}{}\text{(42)}& F({\rho }_{\mathrm{th}},{\rho }_{\exp })={\displaystyle \frac{\mathrm{Tr}({\rho }_{\mathrm{th}}{\rho }_{\exp })}{\sqrt{\mathrm{Tr}({\rho }_{\mathrm{th}}^2)\mathrm{Tr}({\rho }_{\exp }^2)}}}\end{array}

which yielded F({\rho }_{\exp }^{\mathrm{W}})=0.90 and F({\rho }_{\exp }^{\mathrm{GHZ}})=0.92.

In principle, we could use density operator tomography to monitor the variation of the system state during the adiabatic passage. However, this would involve a very large number of measurements and experimental data to be processed. We therefore switched to two simpler indicators of the expected critical behavior, correlation functions and entanglement witnesses. One is the spin-spin correlation including two-spin and three-spin correlation functions. We define an average two-spin correlation as:

\begin{array}{}\text{(42)}& C_2={\displaystyle \frac{4}{3}}{\displaystyle \sum _{i\ne j}}⟨I_x^i{I}_x^j⟩\end{array}

and a three-spin correlation

\begin{array}{}\text{(42)}& C_3=8⟨I_x^1{I}_x^2{I}_x^3⟩\end{array}

The other one is the entanglement witnesses [see Eqs. (32) and (33)].

Figure 16 shows the result of these measurements. The two-spin correlation function C₂, which is frequently used as a measure of the critical behavior, clearly shows the expected transition for the first system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}, but nothing for the second system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}. Inversely, the higher-order correlation function C₃ reveals a critical behavior for H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}, but nothing for H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)}. Therefore, two types of QPTs induced by varying different control parameters are detected by simultaneously employing the multi-spin correlation functions.

The experimentally determined correlations are generally smaller than the theoretical values. This deviation can be attributed to decoherence. To verify this, we also simulated the experiments on a classical computer, using the decoherence model of Ref. [128]. We obtained good agreement between theoretical and experimental data if we assumed an average decoherence time of 150 ms for H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)} and 600ms for H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}. In system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}, the symmetric pulse sequence eliminated the inhomogeneity of the static magnetic field to some extent. The run time of the experiments was about 146 ms for the system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)} and 62 ms for the system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)}. The simulated values are also shown in Fig. 16 (a) and (b), as circles and triangles.

Figure 16 also shows the result of monitoring the QPTs using the entanglement witnesses [parts (c) and (d)]. The data points were obtained by first fitting the experimental signals by an effective decoherence time T_2^{\mathrm{eff}}, which had been estimated by the decay of the signal with the duration of the experiment. These data also allow us to verify the generation of a tripartite entangled state in the system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(1)} as the strength J₂ of the two-body interaction increases, i.e., ⟨{\mathcal{W}}_{\mathrm{GHZ}}⟩>0 and after the critical point J_2=1. Conversely, we observe the generation of a GHZ state in the system H_{3\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}^{(2)} if the strength J₃ of the three-body interaction increases, in agreement with and ⟨{\mathcal{W}}_{\mathrm{W}}⟩>0 when J_3\gg 0.

8 Conclusion

The main goal of this review was an introduction into the basics of quantum information processing and an illustration of its application to the simulation of quantum systems, using nuclear spins as qubits. The examples that we used consisted of small systems in which we adiabatically varied the Hamiltonian in such a way that the systems went through quantum critical points. This allowed us to discuss some important features of entanglement and quantum phase transitions.

Many physical systems or models show a QPT when a single parameter is varied through a critical point, where the ground state and the first excited state become degenerate or go through an avoided crossing. In this paper, we investigated some of these systems by NMR QIP, using adiabatic evolutions under an effective Hamiltonian. A related simulation of a transition from a quantum para- to a quantum (anti-) ferromagnet in a Heisenberg transverse Ising Hamiltonian was also demonstrated with two trapped ions [133].

While the small (toy-) systems that we used here are ideal for demonstrating the basic issues, it is of course tempting to consider larger systems. While this is beyond the scope of this article, we would like to mention some of the issues that will arise in this context.

Generalization. — If we consider systems consisting of a larger number of subsystems, the complexity increases. In particular, larger systems support even more types of qualitatively different entanglement, e.g., nine ways for four qubits [134]. Clearly, the increasing complexity, as well as the increasing size of the relevant Hilbert space, will make it even more difficult to completely characterize these systems. Most of the previous works used the nearest-neighbor pairwise entanglement [106, 108, 115, 135–139] or entanglement entropy [140–146]. However, it has been shown that these measures are not always suitable and should be considered incomplete [147]. Other entanglement measures, such as entanglement witnesses, a global entanglement [109, 148, 149] and average local entanglement \overline{C} [150] will lead to more useful tools for detecting possible phase transitions in multi-spin systems.

Efficiency. — Simulating a QPT on a classical computer is known to be a hard problem. The experimental examples of quantum simulations presented here demonstrate that this simulation can in principle be completed on a quantum computer programmed to pass near a quantum phase transition using adiabatic evolution. The adiabatic simulation of a QPT is thus equivalent to a very specific quantum adiabatic computation [151–154]. Therefore, in analogy to AQC, the performance of the simulation is determined by the minimum energy gap g_{\min } between the lowest two energy eigenvalues of H(t), since the total run time T_{\mathrm{global}}\propto g_{\min }^{-2} in a global adiabatic evolution scheme and T_{\mathrm{local}}\propto g_{\min }^{-1} in a local adiabatic evolution scheme [75, 78] (used in our simulation scheme). If g_{\min } along the adiabatic evolution is only polynomially small in the number of qubits, the simulation is efficient. Schützhold and Schaller [152] suggested that in an AQC, transitions of second or higher order are advantageous in comparison to those of first order. Transitions of second or higher order are thus easier to be simulated in comparison to those of first order.

In our simulation scheme, we introduced a small perturbation into the simulated Hamiltonian to eliminate the degeneracy between the ground state and the first excited state. This allows one to adiabatically go through the critical points. Different choices of time-dependent perturbations, such as Gauss functions, can further improve the energy gap g_{\min } while the phases on the two sides of the critical points remain the same. The implementation of such a scheme requires some prior knowledge of the spectrum. The scaling of the energy gap with the system size is very important for this type of quantum simulations and must be further investigated.

Simulations of quantum mechanical systems may well be the first applications where quantum computers outcompete classical computers, since the requirements on the number of qubits are relatively modest and known quantum algorithms are exponentially faster than classical algorithms. This expectation is also reflected in the European QIPC roadmap [155].

A number of possible approaches to the implementation of quantum information processing have been suggested or demonstrated. Out of these, liquid-state NMR is clearly the most advanced technique. Fully operational systems with up to 12 quantum bits have been demonstrated and used to confirm the basic principles of quantum computing. Even larger system, with up to several thousand qubits have been created using solid-state NMR techniques [156–159]. These systems were designed to study the decoherence of large quantum systems, but they do not allow addressing of individual qubits and therefore cannot be used for information processing.

For the implementation of large-scale general-purpose quantum computers, different types of qubits may well turn out to be more promising. Independent of the physical basis of these future systems, however, it is more than likely that the concepts and techniques for coherent control of coupled nuclear spins will be transferred to other fields of physics. This process has already started in several systems, such as trapped ions [160], excitons in quantum dots [161] and Cooper pair boxes [162].

As Vandersypen and Chuang mentioned “NMR has been the workhorse for the experimental implementation of quantum protocols” [37]. We expect that NMR will remain a testbed for this development for the foreseeable future, especially for the first short term application of quantum computers – quantum simulators.

References

[1] G. E. Moore, Electronics, 1965, 38: 114

[2] P. S. Peercy, Nature, 2000: 406: 1023

[3] L. B. Kish, Phys. Lett. A, 2002, 305: 144

[4] R. Landauer, Phys. Today, 1991, May: 23

[5] R. P. Feynman, International Journal of Theoretical Physics, 1982, 21: 467

[6] P. Benioff, J. Stat. Phys., 1982, 29: 515

[7] E. Bernstein and U. Vazirani, Quantum complexity theory, in: Proc. 25th ACM Symp. Theory Comp., 1993: 11

[8] D. Coppersmith, arXiv: quant-ph/0201067, 1994

[9] Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, Piscataway, NJ: IEEE Press, 1994

[10] S. Lloyd, Science, 1996, 273: 1073

[11] R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Phys. Rev. A, 2002, 65: 042323

[12] D. S. Abrams and S. Lloyd, Phys. Rev. Lett., 1997, 79: 2586

[13] D. S. Abrams and S. Lloyd, Phys. Rev. Lett., 1999, 83: 5162

[14] C. Zalka, Proc. R. Soc. Lond. A, 1998, 454: 313

[15] S. Wiesner, arXiv: quant-ph/9603028, 1996

[16] B. M. Boghosian and W. Taylor, arXiv: quantph/ 9701016v2, 1997

[17] L. A.Wu, M. S. Byrd, and D. A. Lidar, Phys. Rev. Lett., 2002, 89: 057904

[18] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Phys. Rev. A, 2001, 64: 022319

[19] H. Wang, S. Kais, A. Aspuru-Guzik, and M. R. Hoff- mann, Phys. Chem. Chem. Phys., 2008, 10: 5388

[20] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head- Gordon, Science, 2005, 309: 1704

[21] D. A. Lidar and H. Wang, Phys. Rev. E, 1999, 59: 2429

[22] A. Y. Smirnov, S. Savel’ev, L. G. Mourokh, and F. Nori, Europhys. Lett., 2007, 80: 67008

[23] I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-Guzik, Proc. Nat. Acad. Sci. USA, 2008, 105: 18681

[24] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature, 2002, 415: 39

[25] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett., 1998, 81: 3108

[26] J. J. García-Ripoll, E. Solano, and M. A. Martin-Delgado, Phys. Rev. B, 2008, 77: 024522

[27] C. H. Tseng, S. Somaroo, Y. Sharf, E. Knill, R. Laflamme, T. F. Havel, and D. G. Cory, Phys. Rev. A, 1999, 61: 012302

[28] S. Somaroo, C. H. Tseng, T. F. Havel, R. Laflamme, and D. G. Cory, Phys. Rev. Lett., 1999, 82: 5381

[29] C. H. Tseng, S. Somaroo, Y. Sharf, E. Knill, R. Laflamme, T. F. Havel, and D. G. Cory, Phys. Rev. A, 2000, 62: 032309

[30] A. K. Khitrin and B. M. Fung, Phys. Rev. A, 2001, 64: 032306

[31] C. Negrevergne, R. Somma, G. Ortiz, E. Knill, and R. Laflamme, Phys. Rev. A, 2005, 71: 032344

[32] U. Haeberlen and J. S. Waugh, Phys. Rev., 1968, 175: 453

[33] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge: Cambridge University Press, 2001

[34] J. Stolze and D. Suter, Quantum Computing: A Short Course from Theory to Experiment, 2nd Ed., Berlin: Wiley-VCH, 2008

[35] D. Deutsch, Proc. R. Soc. Lond. A, 1989, 425: 1934

[36] D. Deutsch, Proc. R. Soc. Lond. A, 1985, 400: 1934

[37] L. M. K. Vandersypen and I. L. Chuang, Rev. Mod. Phys., 2004, 76: 1037

[38] R. Laflamme, E. Knill, D. Cory, E. Fortunato, T. Havel, C. Miquel, R. Martinez, C. Negrevergne, G. Ortiz, M. Pravia, et al., arXiv: quant-ph/0207172v1, 2002

[39] J. A. Jones and E. Knill, J. Magn. Res., 1999, 141: 322

[40] N. Linden, B. Herve, R. J. Carbajo, and R. Freeman, Chem. Phys. Lett., 1999, 305: 28

[41] N. Sinha, T. S. Mahesh, K. V. Ramanathan, and A. Kumar, J. Chem. Phys., 2001, 114: 4415

[42] A. K. Khitrin and B. M. Fung, J. Chem. Phys., 2000, 112: 6963

[43] J. Du, M. Shi, J. Wu, X. Zhou, and R. Han, Phys. Rev. A, 2001, 63: 042302

[44] K. Dorai, Arvind, and A. Kumar, Phys. Rev. A, 2000, 61: 042306

[45] T. S. Mahesh, N. Sinha, K. V. Ramanathan, and A. Kumar, Phys. Rev. A, 2002, 65: 022312

[46] J. Du, J. Wu, M. Shi, L. Han, X. Zhou, B. Ye, H. Weng, and R. Han, Chin. Phys. Lett., 2000, 17: 64

[47] K. V. R. M. Murali, N. Sinha, T. S. Mahesh, M. H. Levitt, K. V. Ramanathan, and A. Kumar, Phys. Rev. A, 2002, 66: 022313

[48] T. S. Mahesh, K. Dorai, Arvind, and A. Kumar, J. Magn. Res., 2001, 148: 95

[49] N. Linden, H. Barjat, and R. Freeman, Chem. Phys. Lett., 1998, 296: 61

[50] D. G. Cory, M. D. Price, and T. F. Havel, Physica D, 1998, 120: 82. In: Proceedings of the Fourth Workshop on Physics and Consumption

[51] K. Dorai, Arvind, and A. Kumar, Phys. Rev. A, 2001, 63: 034101

[52] R. Das, T. S. Mahesh, and A. Kumar, J. Magn. Res., 2002, 159: 46, ISSN 1090-7807

[53] R. Das, T. S. Mahesh, and A. Kumar, Chem. Phys. Lett., 2003, 369: 8, ISSN 0009-2614

[54] D. G. Cory, A. F. Fahmy, and T. F. Havel, Proc. Nat. Acad. Sci. USA, 1997, 94: 1634

[55] N. A. Gershenfeld and I. L. Chuang, Science, 1997, 275: 350

[56] E. Knill, I. Chuang, and R. Laflamme, Phys. Rev. A, 1998, 57: 3348

[57] L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, R. Cleve, and I. L. Chuang, Phys. Rev. Lett., 2000, 85: 5452

[58] X. Peng, X. Zhu, X. Fang, M. Feng, M. Liu, and K. Gao, Phys. Rev. A, 2002, 65: 042315

[59] X. Peng, X. Zhu, X. Fang, M. Feng, K. Gao, X. Yang, and M. Liu, Chem. Phys. Lett., 2001, 340: 509

[60] Y. Sharf, T. F. Havel, and D. G. Cory, Phys. Rev. A, 2000, 62: 052314

[61] U. Sakaguchi, H. Ozawa, and T. Fukumi, Phys. Rev. A, 2000, 61: 042313

[62] Z. L. Mádi, R. Brüschweiler, and R. R. Ernst, J. Chem. Phys., 1998, 109: 10603

[63] X. Peng, X. Zhu, X. Fang, M. Feng, M. Liu, and K. Gao, J. Chem. Phys., 2004, 120: 3579

[64] E. Knill, R. Laflamme, R. Martinez, and C.-H. Tseng, Nature, 2000, 404: 368

[65] B. M. Fung, Phys. Rev. A, 2001, 63: 022304

[66] X. Peng, X. Zhu, X. Fang, M. Feng, X. Yang, M. Liu, and K. Gao, arXiv: quant-ph/0202010, 2002

[67] W. S. Warren, Science, 1997, 277: 1688

[68] D. Suter and T. S. Mahesh, J. Chem. Phys. 2008, 128: 052206

[69] G. L. Long, H. Y. Yan, and Y. Sun, J. Opt. B, 2001, 3: 376

[70] E. M. Fortunato, M. A. Pravia, N. Boulant, G. Teklemariam, T. F. Havel, and D. G. Cory, J. Chem. Phys., 2002, 116: 7599

[71] R. Das, T. S. Mahesh, and A. Kumar, Phys. Rev. A, 2003, 67: 062304

[72] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science, 2001, 292: 472

[73] A. Mizel, D. A. Lidar, and M. Mitchell, Phys. Rev. Lett., 2007, 99: 070502

[74] M. H. S. Amin, Phys. Rev. Lett., 2008, 100: 130503

[75] J. Roland and N. J. Cerf, Phys. Rev. A, 2002, 65: 042308

[76] M. Steffen, W. van Dam, T. Hogg, G. Breyta, and I. Chuang, Phys. Rev. Lett., 2003, 90: 067903

[77] X. Peng, Z. Liao, N. Xu, G. Qin, X. Zhou, D. Suter, and J. Du, Phys. Rev. Lett., 2008, 101: 145501

[78] A. Mitra, A. Ghosh, R. Das, A. Patel, and A. Kumar, J. Magn. Res., 2005, 177: 285

[79] J. Roland and N. J. Cerf, Phys. Rev. A, 2005, 71: 032330

[80] A. M. Childs, E. Farhi, and J. Preskill, Phys. Rev. A, 2001, 65: 012322

[81] M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness, San Francisco: Freeman, 1979

[82] M. Žnidarič and M. Horvat, Phys. Rev. A, 2006, 73: 022329

[83] T. Hogg, Phys. Rev. A, 2003, 67: 022314

[84] M. Žnidarič, Phys. Rev. A, 2005, 71: 062305

[85] S. Blanes, F. Casas, J. Oteo, and J. Ros, Physics Reports, 2009, 470: 151

[86] M. Suzuki, Quantum Monte Carlo Methods in Condensed-Matter Physics, Singapore: World Scientific, 1993

[87] W. K. Wootters, Phys. Rev. Lett., 1998, 80: 2245

[88] V. Coffman, J. Kundu, and W. K. Wootters, Phys. Rev. A, 2000, 61: 052306

[89] P. Rungta and C. M. Caves, Phys. Rev. A, 2003, 67: 012307

[90] P. Rungta, V. Bužek, C. M. Caves, M. Hillery, and G. J. Milburn, Phys. Rev. A, 2001, 64: 042315

[91] W. Dür, G. Vidal, and J. I. Cirac, Phys. Rev. A, 2000, 62: 062314

[92] B. M. Terhal, Phys. Lett. A, 2000, 271: 319

[93] M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Lett. A, 1996, 223: 1

[94] M. Lewenstein, B. Kraus, J. I. Cirac, and P. Horodecki, Phys. Rev. A, 2000, 62: 052310

[95] A. Sanpera, D. Bruß, and M. Lewenstein, Phys. Rev. A, 2001, 63: 050301

[96] S. Sachdev, Quantum Phase Transition, Cambridge: Cambrige: University Press, 1999

[97] P. C. Canfield, Nature Phys., 2008, 4: 167

[98] H. M. Rønnow, R. Parthasarathy, J. Jensen, G. Aeppli, T. F. Rosenbaum, and D. F. McMorrow, Science, 2005, 308: 389

[99] J. Custers, P. Gegenwart, H. Wilhelm, K. Neumaier, Y. Tokiwa, O. Trovarelli, C. Geibel, F. Steglich, C. Pépin, and P. Coleman, Nature, 2003, 424: 524

[100] A. Yeh, Y. A. Soh, J. Brooke, G. Aeppli, T. F. Rosenbaum, and S. M. Hayden, Nature, 2002, 419: 459

[101] T. Giamarchi, C. Ruegg, and O. Tchernyshyov, Nature Phys., 2008, 4: 198

[102] P. Gegenwart, Q. Si, and F. Steglich, Nature Phys., 2008, 4: 186

[103] S. Sachdev, Nature Phys., 2008, 4: 173

[104] D. M. Broun, Nature Phys., 2008, 4: 170

[105] Editorial, Nature Phys., 2008, 4: 157

[106] A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature, 2002, 416: 608

[107] T. J. Osborne and M. A. Nielsen, Phys. Rev. A, 2002, 66: 032110

[108] M. C. Arnesen, S. Bose, and V. Vedral, Phys. Rev. Lett., 2001, 87: 017901

[109] R. Somma, G. Ortiz, H. Barnum, E. Knill, and L. Viola, Phys. Rev. A, 2004, 70: 042311

[110] S. J. Gu, S. S. Deng, Y. Q. Li, and H. Q. Lin, Phys. Rev. Lett., 2004, 93: 086402

[111] F. Gebbhard, The Mott Metal-Insulator Transition: Models and Methods, New York: Springer-Verlag, 1997

[112] R. B. Laughlin, Phys. Rev. Lett., 1983, 50: 1395

[113] L. Zhou, H. S. Song, Y. Q. Guo, and C. Li, Phys. Rev. A, 2003, 68: 024301

[114] X. Wang, Phys. Rev. A, 2002, 66: 034302

[115] G. Lagmago Kamta and A. F. Starace, Phys. Rev. Lett., 2002, 88: 107901

[116] R. J. Baxter and F. Y. Wu, Phys. Rev. Lett., 1973, 31: 1294

[117] F. Igloi, J. Phys. A: Math. Gen., 1987, 20: 5319

[118] P. Lou, W. C. Wu, and M. C. Chang, Phys. Rev. B, 2004, 70: 064405

[119] P. Suranyi, Phys. Rev. Lett., 1976, 37: 725

[120] C. D Cruz and J. K. Pachos, Phys. Rev. A, 2005, 72: 043608

[121] D. I. Tsomokos, J. J. García-Ripoll, N. R. Cooper, and J. K. Pachos, Phys. Rev. A, 2008, 77: 012106

[122] J. K. Pachos and E. Rico, Phys. Rev. A, 2004, 70: 053620

[123] J. K. Pachos and M. B. Plenio, Phys. Rev. Lett., 2004, 93: 056402

[124] H. P. Buchler, A. Micheli, and P. Zoller, Nature Phys., 2007, 3: 726

[125] J. C. Anglès d’Auriac and F. Iglói, Phys. Rev. E, 1998, 58: 241

[126] K. A. Penson, R. Jullien, and P. Pfeuty, Phys. Rev. B, 1982, 26: 6334

[127] K. A. Penson, J. M. Debierre, and L. Turban, Phys. Rev. B, 1988, 37: 7884

[128] L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, M. H. Sherwood, and I. L. Chuang, Nature, 2001, 414: 883

[129] J. Zhang, et al., Phys. Rev. Lett., 2008, 100

[130] N. Linden, E. Kupce, and R. Freeman, Chem. Phys. Lett., 1999, 311: 321

[131] I. L. Chuang, N. Gershenfeld., M. G. Kubinec, and D. W. Leung, Proc. R. Soc. Lond. A, 1998, 454: 447

[132] G. Teklemariam, E. M. Fortunato, M. A. Pravia, T. F. Havel, and D. G. Cory, Phys. Rev. Lett., 2001, 86: 5845

[133] A. Friedenauer, H. Schmitz, J. T. Glueckert, D. Porras, and T. Schaetz, Nature Phys., 2008, 4: 757

[134] F. Verstraete, J. Dehaene, B. De Moor, and H. Verschelde, Phys. Rev. A, 2002, 65: 052112

[135] O. Osenda, Z. Huang, and S. Kais, Phys. Rev. A, 2003, 67: 062321

[136] Z. Y. Sun, K. L. Yao, W. Yao, D. H. Zhang, and Z. L. Liu, Phys. Rev. B, 2008, 77: 014416

[137] V. Subrahmanyam, Phys. Rev. A, 2004, 69: 022311

[138] F. C. Alcaraz, A. Saguia, and M. S. Sarandy, Phys. Rev. A, 2004, 70: 032333

[139] X. Wang, Phys. Rev. A, 2001, 64: 012313

[140] J. Zhao, I. Peschel, and X. Wang, Phys. Rev. B, 2006, 73: 024405

[141] A. Kopp and K. L. Hur, Phys. Rev. Lett., 2007, 98: 220401

[142] X. Jia, A. R. Subramaniam, I. A. Gruzberg, and S. Chakravarty, Phys. Rev. B, 2008, 77: 014208

[143] L. Cincio, J. Dziarmaga, M. M. Rams, and W. H. Zurek, Phys. Rev. A, 2007, 75: 052321

[144] C. Wellard and R. Orùs, Phys. Rev. A, 2004, 70: 062318

[145] J. I. Latorre, E. Rico, and G. Vidal, Quant. Inf. Comput., 2004, 4: 48

[146] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett., 2003, 90: 227902

[147] M. F. Yang, Phys. Rev. A, 2005, 71: 030302

[148] T. R. de Oliveira, G. Rigolin, and M. C. de Oliveira, Phys. Rev. A, 2006, 73: 010305(R)

[149] D. A. Meyer and N. R. Wallach, Journal of Mathematical Physics, 2002, 43: 4273

[150] X. Peng, X. Zhu, D. Suter, J. Du, M. Liu, and K. Gao, Phys. Rev. A, 2005, 72: 052109

[151] G. Schaller, Phys. Rev. A, 2008, 78: 032328

[152] R. Schützhold and G. Schaller, Phys. Rev. A, 2006, 74: 060304

[153] J. I. Latorre and R. Orùs, Phys. Rev. A, 2004, 69: 062302

[154] T. Caneva, R. Fazio, and G. E. Santoro, arXiv: 0706. 1832v1, 2007

[155] Quantum information processing and communication: strategic report on current status, visions and goals for research in Europe

[156] H. G. Krojanski and D. Suter, Phys. Rev. Lett., 2004, 93: 090501

[157] H. G. Krojanski and D. Suter, Phys. Rev. Lett., 2006, 97: 150503

[158] H. G. Krojanski and D. Suter, Phys. Rev. A, 2006, 74: 062319

[159] M. Lovric, H. G. Krojanski, and D. Suter, Phys. Rev. A, 2007, 75: 042305

[160] R. Gulde, M. Riebe, G. Lancaster, C. Becher, J. Eschner, H. H., F. Schmidt-Kaler, I. Chuang, and R. Blatt, Nature (London), 2003, 421: 48

[161] P. Chen, C. Piermarocchi, and L. J. Sham, Phys. Rev. Lett., 2001, 87: 067401

[162] E. Collin, G. Ithier, A. Aassime, P. Joyez, D. Vion, and D. Esteve, Phys. Rev. Lett., 2004, 93: 157005

[163] Obviously, the populations are not really negative; since we omitted terms proportional to the unity operator, the populations quoted here are just the difference from the mean populations p_av = n^–n.

Rights & permissions

Higher Education Press and Springer-Verlag 2009