# Discrete time quantum walk with nitrogen-vacancy centers

in diamond coupled to a superconducting flux qubit

###### Abstract

We propose a quantum-electrodynamics scheme for implementing the discrete-time, coined quantum walk with the walker corresponding to the phase degree of freedom for a quasi-magnon field realized in an ensemble of nitrogen-vacancy centres in diamond. The coin is realized as a superconducting flux qubit. Our scheme improves on an existing proposal for implementing quantum walks in cavity quantum electrodynamics by removing the cumbersome requirement of varying drive-pulse durations according to mean quasiparticle number. Our improvement is relevant to all indirect-coin-flip cavity quantum-electrodynamics realizations of quantum walks. Our numerical analysis shows that this scheme can realize a discrete quantum walk under realistic conditions.

###### pacs:

05.40.Fb, 03.67.Lx, 03.65.Yz^{†}

^{†}preprint: APS/123-QED

## I Introduction

Spin ensembles that are coupled to quantum electrodynamics (QED) systems and solid-state qubits are promising hybrid architectures Lloyd (2000) for quantum memories Simon et al. (2010) and for distributed quantum networking Müstecaplioğlu (2011a). Nitrogen-vacancy (NV) centers in diamond are outstanding among the solid-state spin ensembles Benjamin et al. (2009) due to their long spin relaxation and dephasing times, scalability Toyli et al. (2010) and existence of well-established techniques for their manipulation Bernien et al. (2012); Gruber et al. (1997); Kagami et al. (2011). NV center quantum memory has recently been demonstrated Fuchs et al. (2011).

Quantum memories require fast and reliable quantum state transfer, which is closely related to the quantum walk (QW) problem Kempe (2003); Venegas-Andraca (2012); Ambainis (2003). In particular, the discrete time QW (DTQW) Shikano (2012); Kitagawa (2012); Venegas-Andraca (2012); Shikano (2013) has been highlighted for universal quantum computation Lovett et al. (2010) and for quantum state transfer Kurzynski and Wojcik (2011). NV centers that are coupled to superconducting transmission line resonators are considered before for continuous-time QW Müstecaplioğlu (2011b). Here we propose a scheme for implementing DTQW with NV centers.

Our scheme is based upon the proposals for DTQW in phase space introduced for cavity QED (CQED) Sanders et al. (2003); Agarwal and Pathak (2005); Xue and Sanders (2008) and circuit QED (cirQED) Xue et al. (2008). Like some other intriguing proposals Lehman et al. (2011); Chandrashekar (2011); Goyal et al. (2012), QWs in CQED and cirQED schemes have not been demonstrated yet; despite some successful realizations in different systems Sansoni et al. (2012); Schmitz et al. (2009); Du et al. (2003); Perets et al. (2008); Broome et al. (2010); Kitagawa et al. (2012). Our procedure could lead to the first solid-state realization of the proposed DTQW in phase space.

We consider an NV center ensemble coupled to a superconducting flux qubit Marcos et al. (2010), which plays the role of the coin in DTQW. Weak excitations of NV centers can be described in terms of bosonic quasiparticles. We take into account all the crystallographic classes of the NV centers so that the single mode coupled to the qubit becomes an effective privileged mode, which is a particular superposition of all the modes corresponding to the quasiparticle excitations in four different NV center crystallographic classes. We derive an exact bilinear spin-boson model that reduces to Jaynes-Cummings model (JCM), same with that in CQED or cirQED settings, under parametric approximation.

We employ the similar but improved DTQW scheme proposed for JCM to our case, by removing the cumbersome requirement of varying drive-pulse durations. This allows for the direct realization of DTQW with regular time steps and simplifies the implementation in all indirect-coin-flip cavity QED realizations of QW including our full solid-state proposal.

This work is organized as follows. Section II presents the model Hamiltonian, describing the coupling of the NV center excitations to the flux qubit. The QW Hamiltonian is obtained by employing canonical transformations. Section III describes how the model Hamiltonian is used to evolve the initial state of the QW using the quantum master equation by taking into account the decay channels from the dephasing and relaxation of the qubit as well as the relaxation of the NV centers. The dynamics of the mean number of NV center quasiparticles (called quasimagnons), the populations of the qubit levels representing the coin operation, and the Holevo standard deviation are reported. The phase distribution and the Wigner function plots are also provided to reveal details of the QW dynamics in the quasimagnon phase space. A linear fit analysis for the Holevo measure is presented and used to distinguish a QW from the classical random walk evolution. Section IV is devoted to the summary and conclusions.

## Ii Model System

### ii.1 NV center ensemble Hamiltonian

A theoretical model of NV centers in diamond can be developed based on the six-electron description of the electron spin triplet ground state of NV Vanoort et al. (1988); Redman et al. (1991). The hyperfine ( MHz for N), quadrupole ( MHz for N), and strain ( MHz) interactions are negligible Jelezko and Wrachtrup (2006) relative to zero-field splitting GHz. Similar conclusion is applicable to C as well Smeltzer et al. (2011). The Hamiltonian of a single NV center in the absence of external electric and magnetic fields can be written as Loubser and Wyk (1978); Doherty et al. (2012)

(1) |

where is the spin-1 matrix and is the zero-field splitting matrix. The NV center Hamiltonian in the principal axes of the molecular body frame becomes diagonal so that

(2) |

where and GHz Manson et al. (2006). By taking the quantization axis as the alignment of the NV center, the Hamiltonian can be expressed in the spin basis as

(3) |

where we dropped the constant term.

There are four possible quantization axes in the NV center, which are all equivalent. We need to take them into account when we consider an ensemble of NV centers. Let us denote these axes by and introduce collective operators, which are analogs of Hubbard operators for strongly correlated electron systems Ovchinnikov and Val’kov (2004), such that

(4) |

Here, is the spin state of the -th NV center aligned along the crystallographic (molecular) axis- which is taken as the quantization axis. denote the number of NV centers aligned along the quantization axis . The Hamiltonian for the NV center ensemble becomes

(5) |

### ii.2 Flux qubit Hamiltonian

A flux qubit consists of three or more superconducting Josephson junctions, one of which is smaller than the others Orlando et al. (1999); Mooij et al. (1999). The phase quantization condition for the superconducting wave functions around the qubit loop reduces the system to a two-dimensional phase space, where the qubit can be effectively considered as a two-particle system subject to a double-well potential Forn-Díaz (2010). A two-level approximation can be made to the band structure of the potential by assuming that the external flux is in the vicinity of half of the flux quantum, where the two lowest bands are closest to each other. A small overlap of the wave functions of the two levels opens the qubit gap , which can be tuned Paauw et al. (2009); Fedorov et al. (2010); Zhu et al. (2010) by replacing the smallest Josephson junction with two parallel junctions, a DC-SQUID. On the basis of the two counter-rotating persistent currents in the qubit loop, the qubit Hamiltonian can be written as

(6) |

where is the magnetic energy bias of the qubit that contributes away from the symmetry point, when there is a net persistent current in the qubit loop.

Pauli spin operators can be conveniently transformed from the bare current basis to the dressed current basis by with . The Hamiltonian becomes

(7) |

with the qubit frequency .

### ii.3 Flux qubit-NV ensemble coupling

We assume the same conditions as in a recent experiment Zhu et al. (2011), wherein the flux qubit was coupled to the NV center ensemble located at the top of the central region of the flux qubit loop. The dimensions of the ensemble are sufficiently small for the spatial variations in the magnetic field over the ensemble to be ignored. In such a case, the magnetic fields generated by the supercurrents are expected to interact with the ensemble homogeneously. In the persistent current basis the electronic Zeeman interaction is given by

(8) |

where are the total spin-1 operators for the NV ensemble that are related to the generators of the SU(3). Nuclear Zeeman interaction is three orders of magnitude smaller and neglected. Here, denotes the coordinate in the molecular body frame; and is associated with the magnetic field in the persistent current basis.

The electronic Zeeman interaction is diagonal in the molecular body frame. The symmetry axes , shown in Fig. 1,

which correspond to four crystalographic classes , respectively, are taken to be the axes of the molecular body frames. The coupling strengths , where is the electronic Landé g factor, is the Bohr magneton and is the -component of the magnetic field acting in the molecular frame of the crystallographic class , can be adjusted by the orientation of the diamond sample. In a typical experimental setting Zhu et al. (2011), the crystal axis is parallel to the magnetic field. The crystal surface is perpendicular to the magnetic field and placed on top of the flux qubit.

If we take the direction as the crystal lab frame axis, then the Zeeman interaction coefficients become proportional to the direction cosines of the lab and body frame axes. The angle between the and axes and the angle between the and axes are the same in the diamond unit cell. The magnetic field axis bisects them and thus the magnetic field makes the same angle of with the every body frame axes. The molecular body frames are chosen such that the axes are perpendicular to the magnetic field and the axes make an angle of with the magnetic field. This makes the same for all crystallographic classes and , associated with the smaller direction cosine, becomes smaller relative to while Marcos et al. (2010); Zhu et al. (2011). These lead to the reduced interaction of the form

(9) | |||||

where H.c. stands for Hermitian conjugate. Here we replaced the collective spin operator by the analogs of the Hubbard operators introduced in Eq. (4) such that and .

We shall proceed under the assumptions of weak excitations and a large number of spins in the ensemble. Experiments can achieve coupling of NV centers directly above the flux qubit Zhu et al. (2011). This number is sufficiently large for our weak-excitation conditions. We consider about excitations in our simulations. Under these conditions, operators can be conveniently related to bosonic operators such that

(10) |

in terms of which the corresponding terms in the Hamiltonian becomes

(11) | |||||

(12) |

where we neglected and redefined in Eq. (12). The operators play the role of Schwinger bosons of subject to the constaint Li and Shen (2004); Shen and Zhang (2002); Bauer et al. (2012). Schwinger bilinear boson representation of the Hubbard operators preserve their algebra. The algebraic relations are summarized in Appendix A.

The quasiparticles of NV center excitations are analogous to quantum particles associated with the spin waves, and hence, we call them quasimagnons. We introduce quasimagnon operators as follows: , and proceed by replacing which is analogous to the parametric pump approximation in quantum optics Mollow and Glauber (1967). Parametric approximation has quite a wide range of validity, which is especially enhanced for initially coherent states D’Ariano et al. (1999). To the first order the error is as large as the of the square root of the mean quasimagnon number. Fractional error, relative to the mean number of quasimagnons, is about the size of relative fluctuation which is negligible in our case . The change from the algebra preserving Schwinger bosons to quasimagnons by the parametric aproximation is equivalent to deforming the Holstein-Primakoff map under group contraction by large total spin approximation Holstein and Primakoff (1940).

We consider initially coherent spin ensembles for our quantum walk procedure as well, which could make our procedure potentially applicable for higher number of excitations, or equivalently for smaller number of NV centers coupled to the flux qubit. A systematic study of how tolerant our procedure to smaller number of NV centers however needs further examinations without parametric approximation, which is beyond the scope of our present investigations.

The Hamiltonian of the NV centers then reduces to

(13) |

We interpret the presence of two quasimagnon operators by associating them with quasispin waves with two possible polarizations. The final Hamiltonian is then an eight-mode model, in which there are four quasimagnon modes, each with two different polarization degrees of freedom.

By employing the rotating wave approximation for , transforming to the dressed current basis at the qubit symmetry point, and using the quasimagnon picture, we can express the interaction as

(14) |

with being the collective coupling strengths for the quasimagnons and flux qubit system. Here are the Pauli spin ladder operators.

This Hamiltonian can be simplified by distinguishing sets of “coupled modes” and “uncoupled modes,” analogous to Morris-Shore transformations Morris and Shore (1983), which are combinations of quasimagnon modes such that they are respectively interacting and noninteracting with the flux qubit. We list these modes in Appendix B. Eventually, only a single effective privileged mode becomes coupled to the qubit. The annihilation operator of that mode is given by

(15) |

The interaction Hamiltonian is then simplified to the usual Jaynes-Cummings model, though the single privileged (bright) mode is actually a spesific superposition of all the quasimagnon modes,

(16) |

where the collective spin-qubit coupling strength is . We ignore the effects of inhomogenous broadening and spatial inhomogeneity here and provide a brief discussion in Appendix C.

### ii.4 Quantum walk Hamiltonian

We introduce an additional term that represents the driving field on the diamond sample

(17) |

where is the drive frequency and is the time-dependent amplitude of the drive pulse. Consistent with our spatial homogenity of the interaction assumption, the drive field is coupled with the quasimagnon mode defined in Eq. (15). We assume that the drive is of square-wave form with pulse duration and amplitude . The total Hamiltonian of the qubit ensemble becomes the Jaynes-Cummings model Hamiltonian such that

(18) |

This Hamiltonian is the same as that proposed for a QW in cirQED, although in our case, the role played by the microwave photons is assumed by the quasimagnons. Here, the qubit-NV center ensemble coupling is about MHz Zhu et al. (2011). In contrast to spontaneous decay and the dephasing of the qubit, the qubit-NV center ensemble interaction is in the strong-coupling regime. By employing a rotating frame transformation with , where we take , the Hamiltonian can be written in the form

(19) |

Here, and are introduced.

The detuning between the qubit and the NV center ensemble is , which can be adjusted by to satisfy so that the Fröchlich canonical transformation Fröhlich (1950); Schrieffer and Wolff (1966); Nakajima (1955) , with the transformation operator

(20) |

can be applied to obtain an effective Hamiltonian

(21) | |||||

where is the unit matrix; and the parameters

(22) |

are associated with the conditional rotational shift and the Hadamard coin terms in the Hamiltonian. The Hadamard coin term only operates when the pulse is on, whereas the shift operation is constantly active. The conditional shift operator rotates the phase states

(23) |

on the phase space such that . The rotation can be clockwise or counterclockwise depending on the state of the qubit . Here are the Fock states in an Hilbert space of dimension . The conditional rotation is restricted to a cycle whose radius is determined by the mean quasimagnon number and the phase space coordinates correspond to the mean quadratures and . The effect of the conditional rotation operator on a coherent state is similar; it translates the phase of the coherence parameter to .

We note that in the QW proposal for cirQED, the parameter is found to depend on by the relation . The dependence of on the drive frequency results in a QW with varying step sizes in time, and complicated pre-adjustments in pulse duration are therefore required. Here, we do not have a drive-frequency-dependent coin operation and we have a fixed step size as well as fixed pulse duration. This makes the implementation of the QW more straightforward. In contrast to our approach described here, a displaced frame transformation relative to the resonator field was previously used in the cirQED model of a bit-flip gate Blais et al. (2007), in addition to the Fröchlich transformation. This lead to vanishingly small cavity photon numbers and negligible ac Stark shifts. In our approach, we do not use a displaced frame, which allows for the presence of a larger number of quasimagnons.

## Iii Results and Discussion

The dynamics of the effective Hamiltonian is subject to decoherence processes. By denoting the qubit dephasing rate as , the qubit relaxation rate as , and the decay of the NV center excitations as , we can write the quantum master equation of the system, under the Markovian and Born-Markov approximations, which are suitable for our weakly coupled system Gardiner and Zoller (2004), as

(24) |

Decoherence processes associated with an operator are described by the Liouvillian superoperators in Lindblad form. We introduced the effective Hamiltonian to illustrate the underlying QW dynamics, but obtain quantitative results from the Hamiltonian of the system given in Eq. (19).

The effective Hamiltonian is in We do not simulate the effective Hamiltonian to avoid further approximations. The effective Hamiltonian only serves to interpret and design the dynamics of the system.

We solve the master equation using the QuTip package Johansson et al. (2012) in python software. The initial condition is considered to be

(25) |

where the quasimagnons are in a coherent state with amplitude and the qubit is in a superposition of its ground and excited states. It is assumed that excitations of the NV center ensemble is initially prepared in a coherent state. We consider a Fock space of quasimagnons with the dimension which is sufficiently large to accommodate the selected in the numerical analysis. A classical magnetic field pulse can be used to prepare the quasimagnon coherent state analogous to the typical driven oscillator case Shumovsky and Müstecaplioğlu (1995), as for large number of NV centers the atomic coherent states converges to the bosonic coherent states Mandel and Wolf (1995). The modes to be included in the initial coherent state could be selected by external bias magnetic fields. The initial coherent superposition of the flux qubit levels can be generated by a resonant magnetic pulse of area .

Most of the parameters we use in our simulations can be taken as the same values with those used in the original circuit QED and cavity QED proposals Xue et al. (2008); Xue and Sanders (2008) except and . Flux qubit has a wide range of tunability and we set its gap frequency as GHz Mooij et al. (1999). Recently developed flux qubits can have dephasing times in the order of few microseconds Steffen et al. (2010) and we set MHz. Recent experiments allow for direct coupling between the flux qubits and the NV centers with strengths of about MHz Zhu et al. (2011). Number of NV centers and the flux qubit size can be used to control the strength of the coupling, which is taken here as MHz.

As our improved procedure is general and applicable to previous circuit QED or cavity QED settings Xue et al. (2008); Xue and Sanders (2008), we shall first use the original values MHz and GHz for comparison of our improved procedure with the original ones in Refs. Xue et al. (2008); Xue and Sanders (2008). The minor difference of no consequence is that we use GHz Manson et al. (2006) for the zero-field splitting between the and ground state spin triplet levels of the NV center. This corresponds to cavity or transmission line mode frequency, which is taken as MHz in Refs. Xue et al. (2008); Xue and Sanders (2008).

Following the comparison, we present the results for the actual flux qubit and NV ensemble parameters. Relaxation rate of the flux qubit is of the same order with its dephasing rate so that we set MHz Steffen et al. (2010). Longitudinal relaxation time of the NV center electron spin is in the order of few milliseconds at room temperature Redman et al. (1991) and gets longer by reducing the temperature. The dephasing time of the spin is in the order of ms and gets longer by increasing the purity of the diamond sample isotropically Balasubramanian et al. (2009). In the case of quasimagnons, the decoherence rate would be collectively enhanced, similar to the collectively enhanced coupling strength. For and ms, we find which is larger than the circuit QED and cavity QED cases where MHz Xue et al. (2008); Xue and Sanders (2008). To compensate it we take larger drive strength GHz.

Relative to the decoherence Yoshihara et al. (2006) and dephasing Kakuyanagi et al. (2007) of the qubit, as well as to the relaxation of the NV centers, the system is in the strong-coupling regime, though still not in the ultrastrong-coupling regime, which would challenge the validity of rotating wave and Born-Markov approximations.

The repetition rate of the driving pulsed laser is taken to be , with being the number of sites on the cycle of a QW. This yields the duration of each time step of the QW as . This time covers the pulse duration of with as well as the time between pulses. Hadamard coin operation induced by the pulse can be described by the Hamiltonian

(26) |

if we take the operation time as ; and if we decouple the quasi-magnon dynamics from the qubit dynamics. For that aim, let us drop the constant term and rewrite the effective QW Hamiltonian in Eq. (21) as

(27) |

under the condition of

(28) |

This condition can be satisfied by taking the pulse frequency to be

(29) |

where is the initial number of quasi-magnons. is the amplitude of the square pulse. If the number fluctuations are sufficiently small, then the dynamics of the quasi-magnons are uncoupled from the qubit degrees of freedom, whose evolution would be a faithful representation of that of an ideal Hadamard coin. Quantitatively, if we denote the number fluctuations by , the condition of validity can be expressed as or , which becomes for our parameters. Our simulations indicate that number fluctuations remain in the order of few quasi-magnons; and thus the coin operation can be simulated. In contrast to the cirQED and CQED QW proposals, we do not change by changing at every step. Furthermore, the pulse duration, and hence, the time steps, are fixed in our approach. These two differences allow for simpler implementation as well as for direct interpretation in terms of DTQW with regular time steps.

We first examine the behavior of the number of quasimagnons, , and the qubit excited and ground level populations, which are respectively defined to be , and . Their time evolutions are shown in Fig. 2. The initial number of quasimagnons is . Fig. (a)a indicates that during the course of dynamics mean number of quasimagnons fluctuates because of the drive and dissipation. Relatively Strong fluctuations in the mean quasimagnon number occur mainly when the pulse is on, or when the Hadamard coin is tossed. When the pulse is off, a conditional phase shift operation characterized by the parameter is set into action. During this rotational QW step in quasimagnon phase space, fluctuations of the are much smaller. The dynamics of the coin, or the populations of the qubit levels, can be followed by the Fig. (b)b. We have marked the regions during which the pulse is on. The dynamics and exchanges between the upper and lower level populations reveal that the qubit evolves in close correlation with the simulation of a pure Hadamard coin. Exchanges only occur when the pulse is on, or at the end of the coin toss.

QW in the quasimagnon phase space needs to be characterized by examination of the dispersion of the angular variable . For cyclic variables the variance is conveniently determined by the Holevo standard deviation which is defined as Holevo (1984)

(30) |

where

(31) |

is the sharpness of the phase distribution. Sharpness is bounded between , for spiked distribution at the mean phase, and , for a flat distribution. It can be interpreted as the success or fidelity of the phase distribution to estimate the mean phase Berry and Wiseman (2000); Bagan et al. (2008). Here we take the mean phase as and choose the domain as . The phase distribution is determined from the reduced density matrix of the quasimagnons by using

(32) |

where the phase states in an -dimensional Fock space are defined in Eq. (23). The quasimagnon reduced density matrix is calculated by, , taking the trace over the qubit degrees of freedom of the density matrix of the system , which is found by solving Eq. (24). If the phase variance is small () then becomes equivalent to the . In contrast to localized distributions, Holevo variance is infinitely large for flat distributions with zero sharpness. In our case, this situation corresponds to large phase diffusion due to decoherence in the system. We examine the Wigner function evolution to determine the number of steps for which the Holevo measure can be used for reliable information on the phase variance.

Fig. 3 shows that increases in an almost smooth manner during the conditional phase shifts, whereas during the coin toss operation, it shows rapid oscillations. At the eighth step, more than half of the cycle is covered by the Wigner function (see Fig. 6). The Holevo measure becomes unsuitable for describing the second order phase variance beyond this point. Similar behavior is encountered in the cirQED QW proposal Xue and Sanders (2008). To check whether the phase variance is faster than that of the classical random walk, we employ a linear fit analysis for the first steps using MATLAB. The result is shown in Fig. 4.

The upper panel in Fig. 4 reveals that there is a linear relation between the number of regular time steps and the logarithm of the Holevo standard deviation. The slope is determined to be , where is the linear-regression standard deviation. This result is close enough to the ballistic case to declare it to be ballistic given that phase diffusion due to decoherence and wrapping due to periodicity of phase have occured. In the case of the classical random walk, the slope would be .

The observation of further QW signatures is to be expected, similar to earlier studies on the same Hamiltonian for similar parameters. For example, the phase distributions at the end of the first four steps are reported in Ref. Xue and Sanders (2008). Interpretation of the phase distribution is not trivial because of the effects of the pump field. It is found that as the pump intensity grows, the distribution becomes asymmetrically skewed in one direction. Moreover, some peaks become narrower. We present the phase distribution in Fig. 5. We used phase states in our numerical calculations. Phase distribution is skewed counterclockwise because of the pumping field, consistent with the observations reported in Ref. Xue and Sanders (2008). Because the walk occurs on a circle, the left and the right boundaries in Fig. 5 are the same. The initial phase distribution of the QW starts at an angular location of . The spread of the distribution is quadratically enhanced, relative to that of a classical random walk.

Wigner functions at the end of the first eight steps are shown in Fig. 6. The spread along the angular direction as the QW progress is accompanied by the increasing deformation along the radial direction due to the number fluctuations. QW splits into other circles at different mean number of quasimagnons and dephasing occurs. Holevo measure increases beyond as the distribution gets more and more flat due to significant contribution of the phase diffusion after steps. Thus we did not use the Holevo measure to characterize the QW for more than steps.

We now examine the effect of decoherence rate of presently available NV center ensembles with MHz Balasubramanian et al. (2009) which is about times larger than the photons in cirQED and CQED QW proposals.

To compensate faster loss of quasimagnons, we take times stronger drive with GHz in the simulations, which is found to be sufficient enough to stabilize the circular trajectory of quantum walkers in the phase space by maintaining as can be seen in Fig. 7. For this drive strength, the drive duration becomes times shorter, due to the relation . The spikes in the figure correspond to the times when the drive is on or the action of coin toss. The Holevo standard deviation shown in the Fig. 8. In this case phase diffusion becomes significant after four steps. When we determine the linear fit in logarithmic scale we find the similar behavior as in Fig. 4 with the slope , where is the linear-regression standard deviation for the first four steps. Such a few number of steps yields large linear-regression standard deviation error in the slope. The slope is nevertheless close enough to be declared to be ballistic, given that phase diffusion due to significantly large decoherence and wrapping due to periodicity of phase have occured.

## Iv Conclusions

We proposed a fully solid-state implementation of a DTQW of NV centers in diamond coupled to a flux qubit. The QW takes place in the phase space of the quasimagnon field by indirectly flipping the flux-qubit coin. The proposed scheme is an improved version of the one originally developed for cirQED and cavity QED systems. In our method, the time steps are of fixed duration and no alteration of the pulse is required. This improved method is applicable to general cirQED and cavity QED QW scenarios. We did not find any significant reduction in the quality of the QW, characterized by the spread of the Holevo standard deviation.

If the NV center ensemble is sufficiently large and is weakly excited, its excitations can be treated as bosonic quasiparticles, called quasimagnons. A driven Jaynes-Cummings model for the quasimagnons and the flux qubit system is established. Further canonical transformation of the Jaynes-Cummings Hamiltonian to a QW on a cycle Hamiltonian is then applied, following a procedure similar to that used for a QW in cirQED. Calculation of the Holevo measure, by taking into account the dissipation channels, and then fitting a curve to it on a logarithmic scale, reveal that the spread of the QW is quadratically faster than that of a classical random walk. Further, QW signatures are explored in the phase distribution and the Wigner function, and features that are similar to those found in the cirQED and cavity QED cases are confirmed.

###### Acknowledgements.

The authors thank William J. Munro, Yuichiro Matsuzaki, Kousuke Kakuyanagi, and Shiro Saito for helpful discussions. Y. S. acknowledges financial support from Grant-in-Aid for Young Scientists (25800181). Ö. E. M. acknowledges financial support from TÜBİTAK (Grant. No. 112T049). B. C. S. acknowledges financial support from AITF, NSERC, and CIFAR. P. X. acknowledges financial support from NSFC 11004029 and 11174052.## Appendix A Hubbard, Schwinger and Quasimagnon Operators

Here we summarize the algebraic properties and relations among the Hubbard operators, bilinear Schwinger boson forms and the quasimagnon operators. For notational simplicity we drop the NV center crystallographic class index . Let us introduce the Gell-mann isotopic spin representation of SU(3) Gell-Mann (1953) in terms of the Hubbard operators,

(33) | |||||

(34) | |||||

(35) |

with and . Each spin- subgroup obeys SU(2) algebra, with ; but they are all correlated according to the full set of commutation relations of the Lie algebra of SU(3).

SU(3) commutation relations can be obtained using the Hubbard operator algebra

(36) |

which is preserved under the bilinear Schwinger map . We associate the collective spin Hubbard operators introduced for the NV center spin ensemble with the Schwinger map subject to weak excitation condition relative to large so that . In particular the and spins become decoupled

(37) |

and contracted into Weyl-Heisenberg algebra according to

(38) |

The bosonization of the and spins in terms of the quasimagnon operators

(39) | |||||

(40) |

physically means that the particle label counting the individual (NV center) spins in the collective spin has lost its meaning for weak excitations in a large ensemble. This allows for quasimagnon description of weak excitations and spin liquid interpretation of NV center ensemble. Mathematicaly this construction is equivalent to contraction of into the semidirect product of spin algebra, which remains unaffected by the bosonization condition, and quasimagnon Weyl-Heisenberg algebra Sun et al. (2003).

## Appendix B Coupled and Uncoupled Quasimagnon Modes With a Flux Qubit

Two quasimagnon modes corresponding to two quasispin waves with different polarizations are nearly degenerate. This allows for decoupling a combination of them from the flux qubit. The annihilation operators for the coupled and uncoupled modes are and , respectively. These are also nearly degenerate, whereby subsequent decouplings are possible. New annihilation operators, given by

are defined to decouple and from the qubit.

When the procedure is applied once more, it leads to the coupling of the qubit to a mode described by the annihilation operator

(42) | |||||

with as the total number of NV centers. The uncoupled mode is given by

(43) |

## Appendix C Inhomogeneous Coupling and Broadening of Nv Centers

Our model was developed under the assumptions of spatially homogeneous coupling of the NV center ensemble to the flux qubit and absence of inhomogeneous broadening. Here we consider the more general case of spatial inhomogeneity for which the model Hamiltonian should be generalized for the position dependent zero field splitting and spatially inhomogeneous coupling .

The NV center ensemble Hamiltonian becomes

(44) |

where the Hubbard operators for each NV is defined to be

(45) |

For the single NV center the weak excitation condition can be used for introducing bilinear Schwinger map with . Using the notation and , we write the as

(46) |

Introducing and dropping the uncoupled mode , the Hamiltonian becomes a spin-boson model

(47) |

We can now identify the bosonic collective mode that is coupled to the flux qubit as

(48) |

where

(49) |

and

(50) |

This shows that quasimagnon picture can be used in the presence of spatially dependent coupling; albeit the quasimagnon mode would be subject to dephasing effect of the inhomogeneous broadening.

The coupling of the flux qubit and the quasimagnon mode leads to a qubit-polariton energy gap. This gap can protect the quasimagnon against dephasing if it is sufficiently strong. A more intriguing proposal would be to exploit a particular inhomogeneous coupling density profile to optimize linewidth reduction Kurucz et al. (2011). This method has been suggested for cavity-spin ensemble coupling Kurucz et al. (2011), where it has been shown that a Gaussian coupling density profile reduce the polariton linewidth to a limit independent of inhomogeneous linewidth. In general, it has been found that the coupling density should fall off as or faster for the manifestation of the polariton gap induced line narrowing.

Another challenge for the experimental implementation of our model in the presence of spatial inhomogeneity would be to drive the exact quasimagnon mode. For spatially homogeneous case this mode is always accessed by the drive. If the drive field spatial profile can be adjusted to match with the magnetic field profile of the flux qubit that couples to the NV center ensemble, then the quasimagnon mode in the spatial inhomogenity case can be driven as well.

## References

- Lloyd (2000) S. Lloyd, arXiv:quant-ph/0008057 (2000).
- Simon et al. (2010) C. Simon, M. Afzelius, J. Appel, A. B. de la Giroday, S. J. Dewhurst, N. Gisin, C. Y. Hu, F. Jelezko, S. Kroll, J. H. Müller, J. Nunn, E. S. Polzik, J. G. Rarity, H. De Riedmatten, W. Rosenfeld, A. J. Shields, N. Sköld, R. M. Stevenson, R. Thew, I. A. Walmsley, M. C. Weber, H. Weinfurter, J. Wrachtrup, and R. J. Young Eur. Phys. J. D 58, 1 (2010).
- Müstecaplioğlu (2011a) Ö. E. Müstecaplioğlu, Phys. Rev. A 83, 023805 (2011a).
- Benjamin et al. (2009) S. Benjamin, B. Lovett, and J. Smith, Laser & Photonics Rev. 3, 556 (2009).
- Toyli et al. (2010) D. M. Toyli, C. D. Weis, G. D. Fuchs, T. Schenkel, and D. D. Awschalom, Nano Lett. 10, 3168 (2010).
- Bernien et al. (2012) H. Bernien, L. Childress, L. Robledo, M. Markham, D. Twitchen, and R. Hanson, Phys. Rev. Lett. 108, 043604 (2012).
- Gruber et al. (1997) A. Gruber, A. Drabenstedt, C. Tietz, L. Fleury, J. Wrachtrup, and C. vonBorczyskowski, Science 276, 2012 (1997).
- Kagami et al. (2011) S. Kagami, Y. Shikano, and K. Asahi, Physica E 43, 761 (2011).
- Fuchs et al. (2011) G. D. Fuchs, G. Burkard, P. V. Klimov, and D. D. Awschalom, Nature Physics 7, 789 (2011).
- Kempe (2003) J. Kempe, Contemp. Phys. 44, 307 (2003).
- Venegas-Andraca (2012) S. E. Venegas-Andraca, Quant. Inf. Proc. 11, 1015 (2012).
- Ambainis (2003) A. Ambainis, Int. J. Quant. Inf. 01, 507 (2003).
- Shikano (2012) Y. Shikano, Quant. Inf. Proc. 11, 1013 (2012).
- Kitagawa (2012) T. Kitagawa, Quant. Inf. Proc. 11, 1107 (2012).
- Shikano (2013) Y. Shikano, J. Comp. Theor. Nanosci.10,1558 (2013).
- Lovett et al. (2010) N. B. Lovett, S. Cooper, M. Everitt, M. Trevers, and V. Kendon, Phys. Rev. A 81, 042330 (2010).
- Kurzynski and Wojcik (2011) P. Kurzynski and A. Wojcik, Phys. Rev. A 83, 062315 (2011).
- Müstecaplioğlu (2011b) Ö. E. Müstecaplioğlu, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (2011b), vol. 8069, p. 9.
- Sanders et al. (2003) B. C. Sanders, S. D. Bartlett, B. Tregenna, and P. L. Knight, Phys. Rev. A 67, 042305 (2003).
- Agarwal and Pathak (2005) G. S. Agarwal and P. K. Pathak, Phys. Rev. A 72, 033815 (2005).
- Xue and Sanders (2008) P. Xue and B. C. Sanders, New J. Phys. 10, 053025 (2008).
- Xue et al. (2008) P. Xue, B. C. Sanders, A. Blais, and K. Lalumiere, Phys. Rev. A 78, 042334 (2008).
- Lehman et al. (2011) L. Lehman, V. Zatloukal, G. K. Brennen, J. K. Pachos, and Z. Wang, Phys. Rev. Lett. 106, 230404 (2011).
- Chandrashekar (2011) C. M. Chandrashekar, Phys. Rev. A 83, 022320 (2011).
- Goyal et al. (2012) S. K. Goyal, F. S. Roux, A. Forbes, and T. Konrad, arXiv:1211.1705 (2012).
- Sansoni et al. (2012) L. Sansoni, F. Sciarrino, G. Vallone, P. Mataloni, A. Crespi, R. Ramponi, and R. Osellame, Phys. Rev. Lett. 108, 010502 (2012).
- Schmitz et al. (2009) H. Schmitz, R. Matjeschk, C. Schneider, J. Glueckert, M. Enderlein, T. Huber, and T. Schaetz, Phys. Rev. Lett. 103, 090504 (2009).
- Du et al. (2003) J. F. Du, H. Li, X. D. Xu, M. J. Shi, J. H. Wu, X. Y. Zhou, and R. D. Han, Phys. Rev. A 67, 042316 (2003).
- Perets et al. (2008) H. B. Perets, Y. Lahini, F. Pozzi, M. Sorel, R. Morandotti, and Y. Silberberg, Phys. Rev. Lett. 100, 170506 (2008).
- Broome et al. (2010) M. A. Broome, A. Fedrizzi, B. P. Lanyon, I. Kassal, A. Aspuru-Guzik, and A. G. White, Phys. Rev. Lett. 104, 153602 (2010).
- Kitagawa et al. (2012) T. Kitagawa, M. A. Broome, A. Fedrizzi, M. S. Rudner, E. Berg, I. Kassal, A. Aspuru-Guzik, E. Demler, and A. G. White, Nature Commun. 3 (2012).
- Marcos et al. (2010) D. Marcos, M. Wubs, J. M. Taylor, R. Aguado, M. D. Lukin, and A. S. Sorensen, Phys. Rev. Lett. 105, 210501 (2010).
- Vanoort et al. (1988) E. Vanoort, N. Manson, and M. Glasbeek, J. Phys. C 21, 4385 (1988).
- Redman et al. (1991) D. A. Redman, S. Brown, R. H. Sands, and S. C. Rand, Phys. Rev. Lett. 67, 3420 (1991).
- Jelezko and Wrachtrup (2006) F. Jelezko and J. Wrachtrup, Phys. Stat. Sol. (a) 203, 3207 (2006).
- Smeltzer et al. (2011) B. Smeltzer, L. Childress, and A. Gali, New J. Phys. 13, 025021 (2011).
- Loubser and Wyk (1978) J. H. N. Loubser and J. A. van Wyk, Rep. Prog. Phys. 41, 1201 (1978).
- Doherty et al. (2012) M. W. Doherty, F. Dolde, H. Fedder, F. Jelezko, J. Wrachtrup, N. B. Manson, and L. C. L. Hollenberg, Phys. Rev. B 85, 205203 (2012).
- Manson et al. (2006) N. B. Manson, J. P. Harrison, and M. J. Sellars, Phys. Rev. B 74, 104303 (2006).
- Ovchinnikov and Val’kov (2004) S. S. G. Ovchinnikov and V. V. Val’kov, Hubbard Operators in the Theory of Strongly Correlated Electrons (Imperial College Press, 2004).
- Orlando et al. (1999) T. P. Orlando, J. E. Mooij, L. Tian, C. H. van der Wal, L. S. Levitov, S. Lloyd, and J. J. Mazo, Phys. Rev. B 60, 15398 (1999).
- Mooij et al. (1999) J. E. Mooij, T. P. Orlando, L. Levitov, L. Tian, C. H. van der Wal, and S. Lloyd, Science 285, 1036 (1999).
- Forn-Díaz (2010) P. Forn-Díaz, Ph.D. thesis, Delft University of Technology (2010).
- Paauw et al. (2009) F. G. Paauw, A. Fedorov, C. J. P. M. Harmans, and J. E. Mooij, Phys. Rev. Lett. 102, 090501 (2009).
- Fedorov et al. (2010) A. Fedorov, A. K. Feofanov, P. Macha, P. Forn-Diaz, C. J. P. M. Harmans, and J. E. Mooij, Phys. Rev. Lett. 105, 060503 (2010).
- Zhu et al. (2010) X. Zhu, A. Kemp, S. Saito, and K. Semba, Appl. Phys. Lett. 97, 102503 (2010).
- Zhu et al. (2011) X. Zhu, S. Saito, A. Kemp, K. Kakuyanagi, S.-i. Karimoto, H. Nakano, W. J. Munro, Y. Tokura, M. S. Everitt, K. Nemoto, M. Kasu, N. Mizuochi, and K. Semba, Nature 478, 221 (2011).
- Li and Shen (2004) P. Li and S. Q. Shen, New J. Phys. 6, 160 (2004).
- Shen and Zhang (2002) S. Q. Shen and G. M. Zhang, Europhys. Lett. 57, 274 (2002).
- Bauer et al. (2012) B. Bauer, P. Corboz, A. M. Läuchli, L. Messio, K. Penc, M. Troyer, and F. Mila, Phys. Rev. B 85, 125116 (2012).
- Mollow and Glauber (1967) B. R. Mollow and R. J. Glauber, Phys. Rev. 160, 1076 (1967).
- D’Ariano et al. (1999) G. M. D’Ariano, M. G. A. Paris, and M. F. Sacchi, Nuovo Cimento 114, 339 (1999).
- Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
- Morris and Shore (1983) J. R. Morris and B. W. Shore, Phys. Rev. A 27, 906 (1983).
- Fröhlich (1950) H. Fröhlich, Phys. Rev. 79, 845 (1950).
- Schrieffer and Wolff (1966) J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491 (1966).
- Nakajima (1955) S. Nakajima, Advances in Physics 4, 363 (1955).
- Blais et al. (2007) A. Blais, J. Gambetta, A. Wallraff, D. I. Schuster, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Phys. Rev. A 75, 032329 (2007).
- Gardiner and Zoller (2004) C. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004).
- Johansson et al. (2012) J. R. Johansson, P. D. Nation, and F. Nori, Comp. Phys. Comm. 183, 1760 (2012).
- Shumovsky and Müstecaplioğlu (1995) A. Shumovsky and Ö. E. Müstecaplioğlu, Phys. Lett. A 209, 88 (1995).
- Mandel and Wolf (1995) L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995).
- Steffen et al. (2010) M. Steffen, S. Kumar, D. P. DiVincenzo, J. R. Rozen, G. A. Keefe, M. B. Rothwell, and M. B. Ketchen, Phys. Rev. Lett. 105, 100502 (2010).
- Balasubramanian et al. (2009) G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, et al., Nature Mat. 8, 383 (2009).
- Yoshihara et al. (2006) F. Yoshihara, K. Harrabi, A. O. Niskanen, Y. Nakamura, and J. S. Tsai, Phys. Rev. Lett. 97, 167001 (2006).
- Kakuyanagi et al. (2007) K. Kakuyanagi, T. Meno, S. Saito, H. Nakano, K. Semba, H. Takayanagi, F. Deppe, and A. Shnirman, Phys. Rev. Lett. 98, 047004 (2007).
- Holevo (1984) A. S. Holevo, Springer Lecture Notes Math. 1055, 153 (1984).
- Berry and Wiseman (2000) D. W. Berry and H. M. Wiseman, Phys. Rev. Lett. 85, 5098 (2000).
- Bagan et al. (2008) E. Bagan, A. Monras, and R. Munoz-Tapia, Phys. Rev. A 78, 043829 (2008).
- Kurucz et al. (2011) Z. Kurucz, J. H. Wesenberg, and K. Mølmer, Phys. Rev. A 83, 053852 (2011).
- Gell-Mann (1953) M. Gell-Mann, Phys. Rev. 92, 833 (1953).
- Sun et al. (2003) C. P. Sun, Y. Li, and X. F. Liu, Phys. Rev. Lett. 91, 147903 (2003).