Quantum process tomography and Linblad estimation of a solid-state qubit

We present an example of quantum process tomography (QPT) performed on a single solid-state qubit. The qubit used is two energy levels of the triplet state in the nitrogen vacancy defect in diamond. QPT is applied to a qubit which has been allowed to decohere for three different time periods. In each case, the process is found in terms of the χ matrix representation and the affine map representation. The discrepancy between experimentally estimated process and the closest physically valid process is noted. The results of QPT performed after three different decoherence times are used to find the error generators, or Lindblad operators, for the system, using the technique introduced by Boulant et al (2003 Phys. Rev. A 67 042322).


Introduction
Quantum process tomography (QPT) is a method of experimentally determining the unknown dynamics of a quantum system. This is crucial as a method of checking the functionality of wouldbe quantum information processing devices. Using knowledge obtained from QPT applied to such a device, one can ideally locate and possibly rectify any source of errors or decoherence. The standard QPT technique, which involves the use of multiple test states, was first developed by Nielsen and Chuang [1] and Poyatos et al [2]. A different technique, exploiting entangled states, was subsequently proposed by Leung [3] and D'Ariano and Lo Presti [4]. More recently, methods have been developed to determine a process from a tomographically incomplete set of measurements [5], as well as techniques to ascertain the master equation describing time evolution of the system [6,7]. QPT has been performed in liquid NMR implementations [8]- [10], numerous optical systems [11]- [16], atoms in an optical lattice [17] and in bulk solid state NMR [18]. This paper contains an example of QPT performed on a single solid-state qubit.

The nitrogen vacancy defect in diamond
The nitrogen vacancy defect (or N-V centre) is a naturally occurring defect in diamond with nitrogen impurities. It can be manufactured in type IB synthetic diamond by irradiation and subsequent annealing at temperatures above 550 • C; radiation damage causes vacancies in the diamond lattice and annealing leads to migration of vacancies towards the nitrogen atoms (see figure 1(a)). The N-V centre can also be produced in type IIA diamonds by N + ion implantation [19,20].
The energy level scheme of the N-V centre is shown in figure 1(b). It consists of a triplet ground state 3 A and a triplet excited state 3 E with a metastable singlet state 1 A. At zero magnetic field, the ground triplet state energy levels are split into degenerate m s = ±1 sublevels and the m s = 0 sublevel. This is due to the magnetic dipole-dipole interaction between the two unpaired electrons. This dipole-dipole interaction can be described by a Hamiltonian term H D = S · D · S where D is the zero-field splitting (or fine structure) tensor and S = (S x , S y , S z ). Using standard techniques [22], this dipolar term can be rewritten as where D and E are zero-field splitting parameters which are dependent on the symmetry of the molecule. The axial symmetry of the N-V centre means that E = 0 MHz, while D for this system has previously been characterized [23,24] as D = 2880 MHz. The end result is that, at zero magnetic field, the m s = ±1 sublevels are degenerate and at an energy 2880 MHz greater than that of the m s = 0 sublevel. With an applied magnetic field aligned along the z-axis of the molecular system (i.e., B = (B z , 0, 0)) the total Hamiltonian, including zero-field splitting term, takes the simple form where the gyromagnetic ratio, gβ, is 2.8025 MHz Gauss −1 . This Hamiltonian leads to a |0 ↔ |1 transition of energy where B z is the applied magnetic field (figure 2(a)), measured in Gauss. The fluorescence intensity corresponding to the m s = 0 transition between the ground ( 3 A) and excited ( 3 E) triplet states is strongest. Optically detected magnetic resonance (ODMR) is used to perform state readout [25]- [27], i.e., the population of the m s = 0 sublevel is identified by the amount of measured fluorescence (figure 2(b)). Manipulation of the qubit state is via standard electron spin resonance (ESR) techniques, using microwave pulses resonant with the |0 ↔ |1 transition (3).
The spin longitudinal relaxation time, T 1 , is of the order of milliseconds at room temperature [28] (relaxation time of the order of seconds is expected at low temperature). Transverse relaxation times (or decoherence times), T 2 , of up to 60 µs have been reported [29], for samples with low-nitrogen concentration.
In the N-V centre, initialization into the m s = 0 state is achieved by optical pumping. Optically induced spin polarization is thought to be related to spin-selective intersystem crossing from the photoexcited 3 E triplet state to the metastable 1 A singlet state (figure 1(b)) i.e., there is an irreversible transition between 3 E → 1 A. The spin polarization achieved by this pumping corresponds to at least 70% population in the m s = 0 sublevel at room temperature. Instead of a pure state, corresponding to 100% polarization, we have what is known as a pseudopure state which takes the general form where n is the number of qubits and 0 < α < 1 is related to the amount of polarization.
In the case of the N-V centre, the 70% population of the m s = 0 sublevel, after optical pumping, corresponds to Tr(ρ pseudo |0 0|) = 0.7 and so It is hoped that, in future experiments, the use of projective readouts, which are possible at low temperatures, will lead to increased polarizations.
A controlled two-qubit quantum gate in which the vacancy electron spin is hyperfine coupled to a nearby C 13 nucleus, has recently been performed using this system [30].
While each of the reconstructed E(ρ j ) will be physically valid density matrices (by appropriate use of quantum state tomography) it is possible that the process E describes is unphysical e.g., the E(ρ j ) are mutually inconsistent or E is not completely positive. Experimental noise and a finite number of measurements to determine expectation values can result in a reconstructed E which is not physically valid.
If it is possible to write a process in the Kraus representation then we can be sure it represents a physically valid process (one that is described by a completely positive map). We can specify an unknown process E by experimentally determining the Kraus operators [31] (6) or, using a fixed (complete) basis of operators {A i }, χ mn is a matrix of coefficients which completely describes the process E and is positive Hermitian by construction. The trace preserving constraint i E † i E i = I becomes m,n χ mn A † n A m = I. • From the measurement results {E(ρ 1 ) . . . E(ρ d 2 )}, the matrix λ jk (which could be viewed as a superoperator) can be determined, given the relation E(ρ j ) = k λ jk ρ k . • In order to determine χ from the matrix λ, one operates on λ with the pseudoinverse of β, where β is derived theoretically from the relationÂ m ρ jÂ † n = k β mn jk ρ k . • Using this last relation and (7), we can see λ ij = mn β mn ij χ mn and so inverting β gives us χ, as required.
To derive {E i } from χ, we first diagonalize χ with a unitary U † Note that this procedure only works if d i 0 which follows from the positivity (semidefinite) of χ. It is this property which we use to check the physicality of a process E; if the χ matrix reconstructed from experimental data has negative eigenvalue(s) then this indicates that noise and/or finitely sampled expectation values has caused the output data to infer an unphysical process. To overcome this problem, a physical matrixχ is found which is as close as possible to the original χ in some sense. Specifically, we used a technique analogous to MLE for state estimation [14,21] by minimizing a deviation function, (t), incorporating a general parameterization for a positiveχ (which, by extension, enforces complete positivity of the process E): where λ is a Lagrange multiplier (which ensures a trace-preserving process) and where T(t) is a d 2 × d 2 complex, lower triangular matrix with d 4 real parameters, t(i). In the case of a single qubit we have: The algorithm used to minimize (t) in (10) was the Nelder-Mead simplex algorithm as implemented in Matlab ® .

Avoiding local minima
As is often the case with numerical minimization problems, (10) typically contains numerous local minima. Our preliminary investigations suggested that, given a random initial point in parameter space, the algorithm would often fail to find the global minimum. As such, a good starting point, t (i), in (12) was deemed necessary if the results were to be meaningful. In order to obtain t (i), we used a technique based on principal component analysis [35]. If the experimentally determined matrix χ is not positive semidefinite, one can 'filter' it by setting any (presumably quite small) negative eigenvalues to zero. Specifically, if we decompose χ as where D is a diagonal matrix containing the eigenvalues of χ, then we can construct a similar but positive matrix χ via where D is identical to D except that any negative eigenvalues have been set equal to zero. From χ we can extract good initial parameters t (i) by performing a Cholesky decomposition (e.g., using the in-built 'Chol' function in Matlab ® ).

Process visualization
Any one qubit state can be parameterized by a so-called Bloch vector r. Equivalently, one can explicitly include the unit coefficient of the identity basis component and parameterize the qubit state by a 4-vector (1, r x , r y , r z ) T : where r 2 x + r 2 y + r 2 z 1 [5]. In this basis any linear trace-preserving evolution E takes the affine form and the process applied to an arbitrary input state is an affine map on r, (although it remains a linear map on ρ): or, more explicitly, Identifying the surface of the unit Bloch sphere as the set of all possible input states, we can visualize a process E by its action on the sphere. The 3 × 3 matrix E is responsible for deformation and rotation of the Bloch sphere, while t denotes displacement from r = (0, 0, 0).
We can gain crude but immediate insight into the nature of a single qubit process by using this visualization technique. Bloch vectors with | r| > 1 do not correspond to any valid density matrix. Hence, when a process is depicted this way, any protrusions of the output ellipsoid outside the unit Bloch sphere mean that the process is not trace preserving.
Complete positivity also places constraints on physically allowable output ellipsoids. For a more detailed analysis of this problem we refer to [36]- [38]. The crucial point to note is that 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT seemingly innocuous ellipsoids (i.e., contained within the unit Bloch sphere and not overly deformed) can represent a process which is not physically valid. The variety of allowable ellipsoids is particularly reduced when the translation vector t in E A is nonzero.

Experimental method
The QPT experiment was performed at room temperature using diamond nanocrystals obtained from type Ib synthetic diamond. Diamond nanocrystals were spin coated on a glass substrate, and single nanocrystals were observed with a homebuilt sample-scanning confocal microscope. In order to ensure the presence of a single defect in the laser focus, the second-order coherence was measured using Hanbury-Brown and Twiss interferometer and the contrast of the antibunching depth was determined. In order to perform many repetitions (to obtain good expectation values) in a reasonable amount of time, a sample with a relatively short coherence time (2 µs) was chosen. Microwaves were coupled to the sample by a ESR microresonator connected to a 40 W travelling wave tube amplifier. A magnetic field ( B = (0, 0, B z ), B z = 200 Gauss) was applied to the defect in order to remove degeneracy of the |±1 energy levels.
Rabi oscillations of fluorescence intensity were used to obtain expectation values. A reference nutation was initially taken in order to normalize the measurement nutations (i.e. set 0, 1 levels for expectation values for the pseudopure initial state). In addition this reference oscillation was used to derive the microwave pulse time required to perform π 2 and π rotations of the Bloch vector. A complete basis of four input states was then prepared using microwave pulses resonant with qubit transitions. The ρ 0 state was obtained directly by optical pumping and three remaining input states ρ 1 , ρ 2 and ρ 3 were obtained by application of suitable π or π 2 pulses. Each of these states was left to decohere for a series of time intervals τ. As a last step, measurements of the diagonal and off-diagonal elements of the density matrix were performed.
Estimates of the density matrix elements were extracted from experimental data (Rabi oscillations) using the maximum entropy (MaxEnt) technique [39]. This method returns a physically valid density matrix which satisfies, as closely as possible, the expectation values of measured observables. In cases where only incomplete knowledge of the output state is known, an additional constraint is used; the reconstructed state must also have the maximum allowable von Neumann entropy.
• Using σ x , σ y and σ z the final (output) states E(ρ i ), after the decoherence, were reconstructed via the MaxEnt technique. • Using these output states, the process for each time period was reconstructed using the technique described above i.e. minimization of (10) The processes, both from raw data and reconstructed (i.e. completely positive), are depicted in figures 3-5.

Quantification of unphysicality
where | | is a projector on to a maximally entangled state in H d 2 i.e., where {|j } is some orthonormal basis set. When the process E is physical, one then obtains a physically valid ρ E which can be compared to the ideal process (converted into ρ id ) using distance measures on quantum states [33]. The trace distance between density matrices ρ id and ρ E is where Similarly one can define the Fidelity: Using this definition, we define the Bures Metric and the C Metric An unphysical process, however, can lead to an unphysical ρ E , possibly resulting in a process fidelity which is greater than one. The application of the preceding fidelity-based distance measures can, therefore, produce nonsensical results. In such cases it is necessary to use other techniques in order to estimate the disparity between χ andχ for example. If one defines X = χ −χ, then possible measures are the matrix p-norms (p = 1, 2, ∞) of X and the Frobenius norm of X ( X Fro ) as well as the trace distance (D pro ). As stated above, these quantities gives a measure of how wellχ describes the experimental results. The difference matrix X varies depending on the basis of operators chosen for χ. We propose that if a standard basis is chosen then the norms of X, defined above, provide a method of 'benchmarking' the quality of all such QPT experiments, regardless of implementation. Table 1 describes the discrepancies between experimental and reconstructed processes. Here, we have chosen the normal basis {A k }(A di+j+1 = |i j|) with which to construct χ andχ. This  (27) and proportional to the characteristic state ρ E for the process:

Markovian process tomography
Standard QPT utilizes a 'black box' approach to studying dynamics. It predicts the resultant output states given arbitrary initial states, but fails to describe the path taken through state space, over time, from initial to final state. If we are prepared to make certain assumptions about the relationship between system and environment, the Born and Markov approximations, then we can construct a Markovian quantum master equation which describes the time evolution of the system over the time studied. In the area of quantum information processing, understanding open system dynamics is important in studying quantum noise processes, designing error correcting codes [40] and locating decoherence free subspaces [41]. We discuss a method of experimentally reconstructing a master equation by using estimates of the state of the system at a number of different timepoints. This particular technique was developed and implemented by Boulant et al [6], in the context of liquid NMR quantum computation. The method is different to standard QPT techniques insofar as it assumes Markovian dynamics and, consequentially, the reconstructed equation has separate terms which describe the unitary and non-unitary aspects of the evolution. We then proceed to apply this technique to the N-V centre, in an effort to better identify the decoherence processes it suffers. By invoking the Born and Markov approximations one may derive the so-called Lindblad master equation [44], commonly expressed in a number of different but equivalent ways, including: The Lindblad master equation (30) is the diagonal form of the GKS master equation derived by Gorini et al [42]: where a αβ is a (d 2 − 1) × (d 2 − 1) Hermitean matrix and {F α } is a linear basis of traceless operators on density matrices. If the matrix of coefficients a αβ (sometimes called the GKS matrix [43]) is positive, then the process described by (32) is completely positive. Diagonalizing a αβ in (32) leads to the Linblad form (30). We can clearly see one advantage gained by making the Born-Markov approximations-the resulting master equation can be separated into components which cause unitary and non-unitary evolution of the system: Non-unitary evolution: 1 2 The operators L k are Lindblad operators and they describe the decoherence processes.
We now describe a technique for estimating the generatorL of a Markovian process. We transform to a Liouville space basis, as in [6], where density matrices become column vectors (denoted |ρ ), and dynamical maps become d 2 × d 2 matrix superoperators (denoted byˆ): The matrix superpropagatorP(t) is then defined by exponentiation: Henceforth, we adopt a similar notation to [6]: becomes ∂|ρ ∂t We callĤ the Hamiltonian superoperator andR the relaxation superoperator. Exponentiating the generator gives us the superpropagator We can manipulate (40) to isolate the relaxation superoperator: An alternative method of derivingR from estimates ofP andĤ is to estimate the derivative at t = 0 of e t 2 iĤ e −(iĤ+R)t e t 2 iĤ .
Here we need to invoke the symmetric Baker-Campbell-Hausdorff formula, for arbitrary operators A and B, Making the identification we get Clearly, differentiating (46) at t = 0 gives usR, as required. In order to obtain this differential, we used a numerical differentiation technique-Richardson extrapolation [6,45].
To use this technique we require estimates of the propagator,P m =P(t m ), at time points t m = 2 m t 1 . This alone, however, is not a robust method for estimatingR. It magnifies the noise present in the estimate ofP 1 [6]. More importantly, there are no constraints on the physicality (i.e. complete positivity) of the generator. This estimate ofR which we will callR RE is instead used as the starting point for a constrained fit to the experimental data.
In order to search for a physical process which was close to the measured data we used a parameterization based on the Gorini-Kossakowski-Sudarshan master equation: Gorini et al have proven [42] thatL is the generator of a Markovian semigroup on d-dimensional Hilbert space if and only if it can be written in the form (47), where the GKS matrix a αβ is a (d 2 − 1) × (d 2 − 1) positive (semidefinite) matrix and {F α } is a linear basis of traceless operators on ρ. We will exploit the requirement that a αβ be positive in a way that is completely analogous to a technique we used in standard QPT. During process reconstruction we enforced positivity ofχ, and, by extension, the complete positivity of E(ρ). Similarly, we will enforce complete positivity of the decoherence processR by constraining the GKS matrix a αβ to be positive. For one qubit we have a parameterization in terms of nine real numbers x(i): −R(x(i)) = 1 2 There is considerable freedom in which basis {F α } to use but we chose to use generators for SU(d), rescaled to be trace orthonormal: A numerical search, using Matlab's ® fminsearch algorithm, was used to find the minimum of whereP m are the experimentally determined propagators, for evolution time t m i.e., where We label the GKS matrix which most closely fits the data, while remaining physical, as a αβ (see figures 6 and 7). To derive the Lindblad operators {L k } from a αβ we first diagonalize a αβ with a unitary U † , where d i are the eigenvalues of a αβ . The Lindblad operators can then be obtained by    figure 11. It is natural now to ask given the estimations for L i , whether one can deduce any information regarding the physical decoherence processes responsible. Unfortunately, this is not generally possible. The most general Markovian master equation for a qubit can be fully described by three Linblad operators. However, this compressed description often will not correspond to the most appropriate physical description of the decoherence. It can be argued however that in most situations, e.g. decoherence free subspace investigation, decoupling pulse generation, all the relevant information is contained in the compressed Lindblad description.
The relative contribution of each Lindblad operator, L i , to the overall decoherence process can be quantified using their Frobenius norms: In summary, the technique we used for ascertaining the relaxation superoperatorR is as follows: • Initialization via optical pumping and appropriate ESR pulses to create a complete basis of input states ρ i . • Each ρ i was left to decohere for a sequence of times t m = 2 m t 1 i.e. {20 ns, 40 ns and 80 ns}.
• State tomography was used to determine the output state ρ i (t m ) after each decoherence time t m . The results were used to construct the matrix superpropagatorsP m .
• Using Richardson extrapolation, an estimate of the generator,R RE , was extracted from the measuredP m and used as a starting point for the fitting procedure (53).
• The minimization of (53) produced a Markovian generator,ˆ R, which was physically valid and best fit the measured data at a sequence of timepoints t m . • The GKS matrix which minimized (53), a αβ , was then diagonalized in order to find the Lindblad operators {L k }.

Summary
We have presented a quantum process tomographic analysis of an individual single solid-state qubit, a nitrogen-vacancy centre in diamond. This analysis is only possible due to the enormous advances made in recent years in single-molecule spectroscopy [25]- [27], where the resultant ODMR technique provides us here with high-fidelity single-qubit readout. As experimental refinements and technological advances improve the readout and control of this system, QPT will become an even more important diagnostic tool. For example, improvement of the m s = 0 polarisation by optical pumping and increased accuracy of rotations by the use of composite pulse sequences, are both feasible short-term goals. If states can be prepared, controlled and read out with very high accuracy, then any deviations from Markovian dynamics (e.g. the disparity in expectation values, between experimental and reconstructed processes apparent in figure 11) cannot be dismissed as noise and the nature of the non-Markovian environment should be investigated. If two-qubit gates can be performed in this system then QPT can be used to verify the robustness of encoded information in decoherence free subspaces. As quantum devices develop and increase in size, the task of 'debugging' the device, or actively identifying the noise present in the device, will pose significant challenges. The paper presented here represents an initial step towards the testing of quantum devices in solid state.