Co-Design quantum simulation of nanoscale NMR

PHYSICAL REVIEW RESEARCH 4, 043089 (2022) Co-Design quantum simulation of nanoscale NMR Manuel G. Algaba ,1,*,† Mario Ponce-Martinez ,1,2,*,‡ Carlos Munuera-Javaloy ,3Vicente Pina-Canelles ,1 Manish J. Thapa ,1Bruno G. Taketani ,1Martin Leib ,1Inés de Vega ,1,2Jorge Casanova ,3,4and Hermanni Heimonen5 1IQM Quantum Computers, Nymphenburgerstr. 86, 80636 Munich, Germany 2Department of Physics and Arnold Sommerfeld Center for Theoretical Physics, Ludwig-Maximilians-Universität München, Theresienstrasse 37, 80333 Munich, Germany 3Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain 4IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009 Bilbao, Spain 5IQM Quantum Computers, Keilaranta 19, FI-02150 Espoo, Finland (Received 14 February 2022; revised 5 August 2022; accepted 28 September 2022; published 8 November 2022) Quantum computers have the potential to efficiently simulate the dynamics of nanoscale NMR systems. In this work, we demonstrate that a noisy intermediate-scale quantum computer can be used to simulate and predict nanoscale NMR resonances. In order to minimize the required gate fidelities, we propose a superconducting application-specific Co-Design quantum processor that reduces the number of SWAP gates by over 90% for chips with more than 20 qubits. The processor consists of transmon qubits capacitively coupled via tunable couplers to a central co-planar waveguide resonator with a quantum circuit refrigerator (QCR) for fast resonator reset. The QCR implements the nonunitary quantum operations required to simulate nuclear hyperpolarization scenarios. DOI: 10.1103/PhysRevResearch.4.043089 I. INTRODUCTION Computer simulations are the backbone of scientific research and technological development. Quantum computers promise in the long term to enable simulations of systems that are intractable to even the largest supercomputers [1,2]. Currently, scientists have access to so-called noisy intermediate-scale quantum (NISQ) computers [3], that present limited qubit counts without error correction. While applications of error-corrected quantum computers are well established, use cases where NISQ devices might achieve quantum advantage are still elusive [4]. In the search for these early applications, the problem must fit the hardware, and the hardware must enable implementation with minimal overheads. Application-specific integrated chips (ASICs) are highly specialized processors optimized for specific problems when execution speed, power efficiency, or miniaturization is of utmost importance [5]. A prominent example where computational speed and energy efficiency are optimised through the use of ASICs is training of artificial neural networks using tensor processing units [6,7]. Building a general-purpose *Both authors contributed equally to this work. †Corresponding author: [email protected] ‡Corresponding author: [email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. quantum computer capable of rivaling the most powerful classical computers has proven to be a difficult task, so it is likely that the first devices reaching useful quantum advantage will use quantum ASICs, also called Co-Design quantum computers. A good example of a problem with suitable structure for simulation by quantum computers is nanoscale nuclear magnetic resonance (NMR) [8]. The problem can be described by a number of mutually interacting spins, which natively map to the qubits of a quantum computer, thereby circumventing the overheads in mapping the problem to qubits, such as in the case of fermions [9]. In general, fast and reliable quantum simulations of interacting spin systems would improve the interpretability of solid-state NMR and electron spin resonance (ESR) spectra, where advanced numerical techniques present very limited performance [10]. This shows the potential of quantum computers with a moderate number of qubits to shed light on the dynamics of these important systems. A Co-Design quantum computer that minimizes algorithm implementation overheads could be the first method to access these simulations. Note that, other NMR problems, such as zero-field NMR [11] and Hamiltonian learning [12], have already attracted research on how quantum computers can be used to tackle them and methods based on Bayesian computation [13] and generative models [14] have been developed for computing NMR spectra as well. NMR techniques have a profound impact in research areas such as material science, chemistry, biology, and medicine [15]. Recently they have approached the nanoscale through solid-state quantum sensors such as the nitrogen vacancy (NV) center in diamond [16]. This is a particularly powerful 2643-1564/2022/4(4)/043089(20) 043089-1 Published by the American Physical Society MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) FIG. 1. NV center with a microwave drive interacting with two mutually interacting 13Cnuclei in a magnetic field  BZ, corresponding to the Hamiltonian in Eq. (1)forM=1andN=2. quantum device, as it enables detection and control of nearby nuclear spins with nanoscale resolution [17]. Applications of the device are, e.g., the precise determination of the structure and dynamics of nuclear ensembles such as proteins [18], finding inter-label distances (via, e.g., Bayesian analysis of the NV response) in electronically labeled biomolecules [19], and the exploration of bespoke microwave (MW) sequences that efficiently transfer NV center polarization to the nuclear environment. Hyperpolarization (i.e., polarization beyond that of a thermal state in a magnetic field) of nuclear spins in diamond presents the potential to develop new and safer contrast agents for magnetic resonance imaging. This problem, which we aim to address through simulation by a quantum computer, could lead to improved detection of different malformations in tissues—such as heart or brain—without the need to deliver ionizing radiation, in contrast to other techniques [20]. This manuscript describes a Co-Design process for a quantum chip able to efficiently simulate nanoscale NMR scenarios. It is structured in three main parts, each of which is a crucial step in the Co-Design process: (1) identifying the problem (Sec. II), which here is simulating a nanoscale NMR system for hyperpolarizing nuclear spins; (2) choosing an algorithm for the nanoscale NMR problem and showing that a star-topology chip implements it with minimal overhead (Sec. III); and (3) Co-Designing the corresponding quantum chip using a central resonator bus (Sec. IV). The sections are followed by results and discussions (Sec. V) and an outlook (Sec. VI). II. NANOSCALE NMR: HYPERPOLARIZATION Let us consider a system consisting of Mnitrogen-vacancy (NV) centers and Ncarbon-13 isotopes in the presence of a driving field and an external magnetic field  BZ. NV centers and nuclei are all effectively described as spin-1/2systems. The representation of such a system for M=1, N=2is shown in Fig. 1. For simplicity, we consider the NV centers aligned with the external magnetic field, leading to the following Hamiltonian: H= M  j=1 δjσz j− N  k=1 ωc k· Ik+ M  j=1 N  k=1 σz j 2 Ajk · Ik + N  k>k gkkIz kIz k−1 4(I+ kI− k+I− kI+ k) + M  j>j hjjσz jσz j−2(σ+ jσ− j+σ− jσ+ j)+Hdr.(1) In Eq. (1), we find the spin operators in the joint Hilbert space C2(M+N)of NV centers and nuclei: σμ j=1⊗···⊗1⊗ jthpos.  σμ⊗1⊗···⊗1   Mfactors ⊗1⊗...⊗1   Nfactors , Iμ k=1⊗...⊗1   Mfactors ⊗1⊗···⊗1⊗ (M+k)thpos.  1 2σμ⊗1⊗···⊗1   Nfactors , where (σμ)2×2,μ∈{x,y,z}is the corresponding 2 ×2 Pauli matrix on the jth NV center and the kth nucleus respectively, and 1is the 2 ×2 identity matrix. Accordingly, σ± j= σx j±iσy j 2(I± k=Ix k±iIy k)arethe jth NV center (kth nucleus) ladder operators. The term δjis the detuning of the jth NV center with respect to the microwave drive Hdr. The hyperfine coupling vector  Ajk represents the coupling between the jth NV center and the kth nucleus, while  ωc k=γc BZ−1 2M j=1 Ajk is the modified Larmor frequency of the kth nucleus with the 13Cgyromagnetic ratio γc≈(2π)×10.7MHz/T, gkkis the coupling between the kth and kth nuclei, and hjjis the coupling between the jth and jth NV centers. Note that, Eq. (1) is expressed in a rotating frame with respect to the free NV Hamiltonian, while Hdr represents an external driving tuned near resonance with a certain NV energy transition. The derivation of Eq. (1) can be found in Appendix A. In order to hyperpolarize a diamond sample at room temperature, the NV centers are first optically polarized employing laser light, and then their state is transferred to the surrounding nuclei with the aid of a tailored microwave radiation scheme. The initial state of the nuclei in a roomtemperature sample is well described by a fully mixed state due to the small energy splitting of the nuclear spins. By reinitializing the NV centers and repeating this procedure, the polarization transferred into the sample can be amplified. In this paper, we will consider the quantum simulation of the polarization transfer mechanism and study two different driving schemes acting on the NV centers in a room-temperature diamond. The first driving scheme is a continuous driving whose Hamiltonian in the rotating frame mentioned earlier is Hdr =  2σφ, where σφ=e−iφ|10|+eiφ|01|=e−iφσ−+eiφσ+, φa phase, the Rabi frequency and the kets |1and |0 are the eigenvectors of the operator σzwith eigenvalues ±1 respectively. The set {|0,|1} is called the computational basis of the state space of a two level system, and will be 043089-2 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) our standard choice for a basis, |0≡(1,0)tand |1≡(0,1)t. NV-nucleus polarization transfer is achieved when the Rabi frequency matches the modified nuclear Larmor frequency (i.e., when =| ωc|), leading to the Hartmann-Hahn doubleresonance condition [21]. For a single NV center and nucleus, the Hamiltonian in Eq. (1) reduces, in an interaction picture, to HI=A⊥ 4(|+−|I++|−+|I−), where ±| = 0|±1|, which shows a polarization transfer mechanism with the effective transfer rate A⊥ 4(a detailed derivation can be found in Appendix B). The second type of driving we consider is a pulsed-driving scheme, Hdr =(t) 2σφ, where (t)isatrainofπ-pulses, such as the Carr-Purcell-Meiboom-Gill sequence [22,23]orthe XY8 sequence [24,25]. We consider pulses with a negligible width compared to the time spacing τbetween the π-pulses. If τis selected such that τ=nπ | ωc|(nbeing an arbitrary integer number) and the pulses are evenly spaced one finds that, in an interaction picture, for a single nucleus and NV center, the Hamiltonian reduces to HI=αA⊥σzIx, where αis a factor that depends on the integer n(see Appendix B). A phase imprinted on the pulse sequence through a time delay turns the interaction into HI=αA⊥σzIy. By combining both sequences with the appropriate rotations over the NV center, the polarization transfer interaction HI=−αA⊥ 4(σ+I−+σ−I+) is achieved (see Appendix Band Ref. [26] for more details). Regarding common error sources, NV centers located at different positions in the diamond lattice experience stress conditions that lead to local energy deviations from the zerofield splitting. The corresponding term in Eq. (1)isthe detuning δj. Another common type of imperfection appears due to unavoidable fluctuations of the Rabi frequency of the driving. This fluctuation can be modelled as an OrnsteinUhlenbeck (OU) process [27], which has been shown to be an accurate description for NV centers [28]. It is a Gaussian process of the following form [29]: X(t+t)=X(t)e−t/τ +cτ 2(1 −e−2t/τ )1/2 N(t),(2) where tis the time step, τthe correlation time, cthe diffusion constant of the process, and N(t) a temporally uncorrelated normally distributed random variable. It is a dimensionless term, which yields an effective Rabi frequency of (1 +X). Neither of the system error types lead to considerable overheads in a simulation on a quantum computer. Finally, 13Cnuclear spin decay is not a relevant error source on the time scale of the protocol, since it is of the order of seconds [30], while the hyperpolarization process operates in the order of microseconds. III. CO-DESIGN ALGORITHM In this section, we provide an in-depth description of our Co-Design algorithm, starting with the choice of a simulation technique, followed by a short listing of hardware assumptions related to the allowed qubit operations (gates and resets), as well as the noise and errors present in the physical NMR system and in the quantum computer. Subsequently, the algorithm components are introduced. We end the section with a discussion on layout and gate-level optimization. The highlevel structure of the simulation protocol is shown in Fig. 2(a). A. Simulation technique The best established digital quantum simulation technique is based on decomposing the time-evolution operator into single-qubit and two-qubit gates through the Lie-TrotterSuzuki formula [31], known as Trotterization. To simulate our problem on a quantum computer, we base our strategy on regular Trotterization [2] but we also explore the randomized Trotterization method qDRIFT [32] in Appendix C. Other, more NISQ-specific, simulation techniques such as the variational quantum simulator [33], the quantum assisted simulator [34], numerical quantum circuit synthesis [35], and a plethora of other quantum algorithms [4] can also be used as simulation methods. One advantage of Trotterization over some of these NISQ methods is that it closely follows the real time evolution for each time step. This is particularly important for pulseddriving schemes, where the free evolution in between different pulses always starts with a different initial state. Variational and quantum assisted methods would then require that each interpulse evolution is solved independently, making them impractical for the problem. A second advantage of Trotterization is that its complexity and precision are straightforward to analyze. The Trotterization procedure can also be expanded to higher orders, and symmetrized expansions converge more rapidly and reduce the error with respect to the continuum time limit [36]. B. Hardware assumptions 1. Native gates The hardware for the quantum simulation plays a major role in choosing the optimal quantum algorithm and its specific implementation. In our case, we consider a quantum computer based on superconducting qubits with the following native single-qubit gate set: Rxy(φ,θ)=e−i(cos φX+sin φY)θ 2and (3) Rz(θ)=e−iZ θ 2,(4) where X,Y, and Zare Pauli operators on the superconducting transmon qubits. The Rxy(φ,θ) can physically be implemented through a microwave drive [37]. The gate Rz(θ)onthe other hand does not need to be implemented directly, but can be performed virtually by tuning the phase of the subsequent gates applied on the qubit [38]. This reduces the number of single-qubit gates (SQGs) that need to be implemented. The native two-qubit gate (TQG) that arises from the superconducting system Hamiltonian shown in Sec. IV and Appendix G, is a continuously parameterized controlled-Z (CZ) interaction [39], which can be transformed through local virtual Rzrotations into the form of a ZZ interaction: UZZ(φ)=⎛ ⎜ ⎜ ⎝ e−iφ00 0 0eiφ00 00eiφ0 000e−iφ ⎞ ⎟ ⎟ ⎠.(5) 043089-3 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) FIG. 2. (a) Sketch of the overall operation of the simulation algorithm for one NV center and two nuclei, with continuous driving; (b) corresponding gate sequence of one Trotter step on a star-topology chip for noninteracting nuclei. HSQG refers to the single-qubit-gate component of the Hamiltonian, Ax,y,z j,1parameters are the various coupling strengths of the simulated system, and the Xrnd gates refer to X gates applied with a 50% probability to prepare an effective fully-mixed state. The initial state preparation can also be performed using the alternative random-phase approximation-inspired method. NV init is an initial-state preparation using single-qubit gates to the state required by the driving scheme. Details of the circuit components can be found in Appendix Dalong with a figure representing the pulsed driving case. Even though the ZZ-interaction and the controlled-Zinteractions appear different, their physical implementation is identical since they are related through virtual Rzrotations which come at no additional cost. Section IV goes into more depth on the two-qubit-gate implementation on our Co-Design quantum chip. 2. Qubit reset In the hyperpolarization process, the state of the NV needs to be re-initialized after each cycle. It is therefore necessary to be able to reset the state of the qubit representing the NV center in the quantum computer. A qubit reset operation can be defined by two Kraus operators: Kreset 1=10 00 ,Kreset 2=01 00 .(6) On superconducting hardware this can be realized through connecting a quantum circuit refrigerator (QCR) to each circuit element that needs to be reset [40–43]. Different reset schemes are discussed in Sec. IV B. 3. Noise and errors In this paper, we show that the simulation can tolerate the noise of the quantum processing unit (QPU), and that the simulation does not require large overheads to implement imperfections present in the nanoscale-NMR system, as discussed in Sec. II. We will refer by systemimper fections to effects in the nanoscale NMR system only, while the QPU is affected by noise, referring to the effect of the environment on the qubits, and errors, referring to inaccuracies of gates. In our simulation of the algorithm, we use the most common noise models for superconducting transmon qubits [37], namely, an amplitude damping channel modelled by the Kraus operators: Kamp 1(t)=|00|+1−p(t)|11|=10 0√1−p(t), Kamp 2(t)=p(t)|01|=0√p(t) 00 ,(7) with p(t)=1−exp(−t/T1) and T1=60 μs, and a pure dephasing channel represented by the Kraus operators: Kdeph 1(t)=10 0√1−p(t),Kdeph 2(t)=10 0√p(t), (8) with p(t)=1−exp(−(t)) and (t) given by the expression (t)=t2 2∞ 0dωI(ω) cotanh( βω 2)sinc 2(ωt 2) where βis the inverse temperature of the environment. We chose the spectral function I(ω) to be of the type 1/f[37], and T2=60 μs.Additionally, each gate operation is assumed to be calibrated up to a two-qubit-gate (TQG) error εTQG ∈[10−4,10−2], with the induced effective noise modelled by a depolarizing channel 043089-4 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) defined for single-qubit gates by the Kraus operators: Kdepol 1=1−pI, Kdepol 2=p/3X, Kdepol 3=p/3Y, Kdepol 4=p/3Z, (9) and for two-qubit gates by an analogous expression with the tensor products of two Pauli matrices and the coefficients √1−pfor the identity and √p/15 for the other operators. Single-qubit-gate (SQG) errors εSQG areassumedtobeone order of magnitude lower than TQG errors. C. Algorithm components Our simulation of the nanoscale NMR problem follows the general structure shown in Fig. 2(a). It starts by initializing the states of all qubits, according to whether they represent a nucleus or a NV center, then evolving them using Trotter steps, followed by reset and re-initialization of the qubits representing NV centers. The cycle of time evolution and re-initialization is then repeated as many times as the protocol calls for. Finally the qubits are measured, and the polarization of the NV centers and nuclei are extracted as the expectation values of the qubit representing each element. Figure 2(a) shows the circuits for the case of continuous driving, while the details of pulsed driving schemes are shown in Fig. 10(a) in Appendix B. In the following, we go through these steps in more detail for the case of a single NV center. 1. Initial state preparation To enable the polarization transfer, it is necessary to prepare the NV center in a specific initial state that depends on the driving scheme. For the continuous-driving scheme, it is the |+ or |− state, and for the pulsed-driving scheme, it is one of the two computational basis states, |0or |1. For a diamond at room temperature, the initial state of the nuclear spins is well described by a fully mixed state ρmixed =1⊗N 2N, where 1⊗Nis the 2N×2Nidentity matrix. The state can be approximated by running the algorithm several times, each time with a different initial state obtained by applying Xgates randomly on the qubits representing nuclei. A faster alternative to this sampling is the randomphase-approximation-inspired method, described in Ref. [44], and introduced into quantum computing in Ref. [45]. In this method, the qubits are all prepared in an equal superposition by applying Hadamard gates, and then the phases are randomized through the application of random phase gates. The method effectively reduces the prefactor in the scaling of the sampling error [45]. 2. Time evolution We choose to implement the time evolution generated by the Hamiltonian in Eq. (1) through Trotterization. For that, the Hamiltonian is rewritten in terms of qubit Pauli operators and arranged into noncommuting terms for an optimal Trotter splitting. The resulting circuit, which performs one Trotter step of the evolution in the continuous driving case, is depicted in Fig. 2(b). It consists of a set of initial single-qubit gates, including the ones corresponding to the driving and the detuning of the NV center, followed by three two-qubit gates per nucleus. There are three types of interaction terms, of the form XZ,YZ and ZZ, when no internuclear interactions are considered. With interactions there are a total of five interaction terms. Our native gate set only includes one type of two-qubit interaction as explained in Sec. III B 1. Therefore some SQGs need to be applied in order to convert the interaction terms into the right form, as discussed in Appendix D. Under specific circumstances, some TQGs can be removed by rotating the Hamiltonian into a more suitable basis as shown in Appendix E. 3. Cycles and reset The dynamics of the system is known to produce an exchange of polarization between the NV center and the nuclei. This exchange is oscillatory, and therefore choosing a proper stopping time is important in order to achieve an effective polarization transfer from the NV center to the nuclei. In practice, a suboptimal transfer time can suffice, and the protocol is then repeated several times by resetting the NV center to its initial state and letting the system evolve under the drive again. Due to the re-initializations the full evolution of the system is nonunitary and a net gain of polarization of the system is enabled. This structure is represented in the quantum circuit in Fig. 2(a) by the repeated Trotter evolution, followed by reset operations on the qubit representing the NV center, and a single-qubit gate to prepare the initial state of the driving protocol. D. Layout optimization When implementing a quantum algorithm on a superconducting QPU, the planar qubit connectivity forces us to solve the qubit-routing problem by introducing additional SWAP gates to connect distant qubits. In this section, we study the advantages of an optimized chip topology, a star topology, over a square-grid array of qubits in terms of reducing the number of SWAP gates that must be inserted to run the algorithm in Fig. 2on the device. Different topologies will imply different counts of SWAPs added on top of the gates arising from the algorithm itself, as shown in Fig. 3.OnaNISQ device, this implies different computational precision for the same gate error magnitudes. We choose the SWAP count as our metric to compare different topologies, as commonly gates have fidelities limited by calibration. The errors could be due to crosstalk, leakage, or filtering causing disturbances to the control signals. Under this scenario we want to minimize the gate count. On the other hand, for a highly tuned up device whose gates are limited by qubit coherence times, it would be optimal to minimize the circuit depth instead of the TQG count. Assuming the gate errors are independent, the total error will be bounded by: εgates =1−(1 −εTQG )NTQG (1 −εSQG)NSQG ,(10) where NTQG is the number of two-qubit gates, NSQG the number of single-qubit gates, and εSQG is the SQG error. 043089-5 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) FIG. 3. (a) Three steps of the SWAP patterns in a five-qubit linear chain displayed from top to bottom. Green (blue) arrows represent the SWAP pattern for the case with (without) internuclear interactions. The green pattern is known as the “odd-even” SWAP pattern. The numbers are expressed according to the blue pattern, where label 0 represents the position of the NV center. (b) Star chip topology with the SWAP pattern for the interaction with internuclear interactions. Consequently, reducing the gate count, especially NTQG, has an exponential effect on the precision of the computations, underlining the effect of minimizing the SWAP gate overhead. As SWAP gates are not native to the hardware, but must be compiled out of three CZ gates, their effective error rate is also much higher than those of native gates. 1. Square grid A common choice in superconducting quantum chips is the square grid of qubits. It has high connectivity and is suitable for performing the surface code error correction when scaled to large enough qubit counts with fast measurement and feedback [46]. The qubit routing problem on a square grid can be tackled using various numerical approaches [47–50]. However, these methods are inefficient. In our case, a tailored SWAP routing method, shown in Fig. 3(a), has been chosen and developed in Appendix Fthat can be shown to be well suited from two perspectives. First, a comparison against the cited numerical approaches (shown in Appendix F) reveals FIG. 4. The (top) panel shows the percentage of SWAP gates saved by using a star topology instead of a square grid for nqubits for the cases with and without internuclear interactions. The (bottom) panel shows the total TQG count against the qubit count in the interacting case for the square grid and the star architecture. that our routing method is better in terms of number of gates. Second, it is completely deterministic and does not rely on expensive numerical optimization methods. It can also be shown not to be far from optimal: on a square grid each qubit has at most four nearest neighbors, implying that any SWAP operation provides at most three new neighbors. For an all-to-all (ATA) interacting Hamiltonian there are n2/2 interactions, to leading order, for a simulation performed on nqubits (corresponding to Nnuclei and one NV center). This implies a lower bound of at least n2/6 SWAPs for any SWAP pattern on the square grid topology. Our SWAP pattern with n2/2 SWAPs, discussed in Appendix F, is thus not far from optimal. 2. Star architecture A star topology allows to implement the simulation of the simplified case without internuclear interactions directly, without any SWAP gates. With internuclear interactions considered, we still find a reduction in SWAP gates as compared to the square grid topology, as shown in Fig. 3(b). This reduction comes from the SWAP routing we implement, that consists of making the qubit 0 in Fig. 3(b) interact with all the external qubits and then swap its state with that of qubit 1 and repeat this process until all interactions have been performed. This allows us to use only n−1 SWAP gates. The percentage of SWAP gates that can be saved can be observed in Fig. 4. However, this improvement in the number of gates comes with a price to pay in the depth of the algorithm. We can only 043089-6 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) TABLE I. Overheads introduced by the decomposition of UZZ(φ) gates into different examples of native TQGs in superconducting devices. The single-qubit-gate (SQG) count includes only Rxy rotations, as the Rzrotations can be implemented virtually. UZZ (φ)UZZ (−π/4) CZ(π)CNOT TQGs 1 2 2 2 SQGs 0 5 3 1 do one TQG at a time in the star chip and we have 3 2n(n−1) TQGs from simulating the physical interactions and 3(n−2) TQGs from the SWAPs. This yields a depth for the TQGs of 3 2n2+3 2n−6 in a star chip, while for a square grid it is 6n. Such depth increase comes from the reduction in parallelization, since all gates now act via the central qubit. On the other hand, less parallelization reduces the types of possible crosstalk errors. Adding connections between external qubits reduces the depth of the circuit, since the main cause of circuit depth is the fact that the interaction of two external qubits needs to be done exclusively by the central qubit. Further studies are required to see if the addition of more external layers to this topology (such as in a spider web) can lead to better compromises between depth and gate count, especially for simulating systems with clusters of strongly interacting nuclei. E. Gate-level optimization The two-qubit interactions that appear in the algorithm are the XZ,YZ and ZZ interactions, as shown in Sec. III C 2 and Fig. 2(b). When compiling the algorithm into the native gates of the device, all these interactions must be implemented in terms of some available gate set. We study in Table I the overhead introduced by decomposing these interactions into different examples of native TQGs of superconducting devices; namely, the parametrizable and fixed-phase UZZ gate, the fixed-phase controlled-Zgate CZ, and the CNOT gate. The CNOT gate is usually performed by making use of the crossresonance gate [37,51], which introduces an UXZ interaction, making it equivalent to the UZZ for the purpose of this algorithm. We assume that the SQGs that can be implemented are the Rxy and the Rzgates. These numbers can be further reduced if the first and last SQGs introduced by this compilation are combined with the adjacent SQGs in the algorithm. The conclusion is that fixed-angle gates will double the number of TQGs that need to be physically performed. In Ref. [52], the improvements coming from the reduction of the gate count are compared to the new errors introduced by the interpolation of the calibrated phases. For two instances of a Quantum Approximate Optimization Algorithm (QAOA) [53], it is shown that the performance is better when using parametrized TQGs. The gate sequences for some of the gate decompositions are shown in Fig. 5. IV. CO-DESIGN HARDWARE A star-architecture chip has fundamental scaling issues using a transmon as the central qubit as the number of neighFIG. 5. (a) Gate decomposition of e−iφZZ in terms of the fixedphase UZZ (π 4) gate, (b) the CNOT. bors grows. Every neighbor added to the center qubit would decrease its charging energy Ec. To keep the qubit frequency constant and anharmonicity in the transmon regime, the ratio of the qubit’s Josephson energy to its charging energy, Ej/Ec, must remain unaffected. Therefore we cannot afford to change its charging energy. This leads to a trade-off between the number of coupled qubits and their coupling strength to the central element. The spirit of Co-Design calls for replacing the central transmon with another object that enables this scaling in size. A resonator has no Josephson energy Ej,sotheEj/Ecratio is not altered by adding more capacitive couplings to the resonator. Only small corrections to its frequency are introduced by adding coupled qubits. As a distributed element, a co-planar waveguide resonator also has physically more space for couplings than a central transmon qubit. By elongating the resonator and choosing the mode with the target frequency, the number of qubits coupled to it can further be increased. These properties make a resonator a favourable component in the center of the chip. In the device in Fig. 6(a), the qubits are capacitively coupled to the resonator via tunable couplers [39,54,55]inthe proximity of a voltage maximum of a standing wave in the resonator. As the resonator is elongated, we must use higher harmonic excitations of the resonator to keep the frequency around the operational frequency of the qubits. Tunable couplers avoid the frequency crowding issues related to direct coupling [56,57], and the linear resonator has higher connectivity in the center than ring resonator structures with quasi-all-to-all connectivities [58]. A linear resonator cannot in general be used as a qubit, since a microwave drive on it will not only populate the {0|,1|} subspace, but also higher excited states. However, the effective interactions mediated via the tuneable coupler in Fig. 6(a) are of the type a†aZ and (a+a†)X+(a−a†)Y where aand a†[37] are the resonator creation and annihilation operators. These types of interactions conserve excitation number, so when at most one excitation is in the qubitresonator system, the resonator cannot be populated beyond its first excited state through interaction with a qubit mediated a tuneable coupler. CZ and iSWAP gates between the resonator and a qubit can be performed using the two interactions, and the theory is developed more fully in Sec. IV A. Then, a resonator together with an external qubit can be used as an effective central qubit in the following way. 043089-7 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) FIG. 6. (a) Central λ/4 resonator with 6 qubits coupled via tunable couplers. The resonator is also coupled to a quantum circuit refrigerator enabling fast reset. The device acts effectively as a six qubit star-architecture chip. (b) Electrical diagram of transmon qubit (left) coupled to a resonator mode (right) via a tunable coupler (center). The qubit has frequency ωq, coupler ωc, and resonator ωr.The qubit and resonator have a direct capacitance Cqr and capacitances Cqc and Crc respectively to the coupler. 1. Prepare all qubits and the resonator in their ground states 2. Select one qubit to form the effective central qubit together with the resonator 3. Prepare an arbitrary state in the selected qubit 4. Perform an iSWAP operation from the selected qubit to the resonator initially in the ground state 5. Perform CZ gates between the resonator and any other qubits 6. Perform an iSWAP operation back from the resonator to the selected qubit for measurement The theoretically most straightforward protocol would be to perform a SWAP gate from the qubit to the resonator. The iSWAP, on the other hand, is a native gate that can directly be implemented on the hardware in Fig. 6(b).TheiSWAP gate between the resonator and the qubit is represented by the unitary operator: UiSWAP =⎛ ⎜ ⎝ 10 00 00−i0 0−i00 00 01 ⎞ ⎟ ⎠.(11) TABLE II. Parameters of star-architecture chip. Parameter Symbol Value Resonator frequency ωr2π×4.3 GHz Qubit anharmonicity αq−2π×0.187 GHz Coupler anharmonicity αc−2π×0.110 GHz Resonator-coupler coupling grc 2π×98.5MHz Qubit-coupler coupling gqc 2π×101.8MHz Resonator-qubit coupling grq 2π×8.9MHz Resonator relaxation Tr 160 μs Qubit relaxation Tq 160 μs Coupler relaxation Tc 130 μs Resonator dephasing Tr 260 μs Qubit dephasing Tq 260 μs Coupler dephasing Tc 230 μs Since the CZ gates performing the computation following the iSWAP are diagonal in the computational basis, the phase introduced by the iSWAP is uninvolved in the gate. This enables substituting the SWAP gate by an iSWAP gate in the protocol to further minimize the gate count. A. Gate theory and simulations Here we demonstrate that in our star architecture CZ and iSWAP-type gates between any of the qubits and the {0|,1|} subspace of a chosen resonator mode can be implemented. The operational principles of these gates are very similar to those between two qubits coupled with a tunable coupler [39,54,55,59]. The main limitation of our architecture (where one transmon is replaced by a resonator) is that iSWAP operations can only be performed in the zeroand single-excitation subspace of the two-qubit computational basis. 1. Conditional-Z gate The CZ operation between the resonator and the qubit is described by the unitary operator: CZ(φ)=⎛ ⎜ ⎜ ⎝ 100 0 010 0 001 0 000e−iφ ⎞ ⎟ ⎟ ⎠.(12) This gate is equivalent to the UZZ(φ) gate in Eq. (5)up to two Rzrotations. To operate a CZ gate, we initialize the resonator-coupler-qubit set up shown in Fig. 6(b) at the idling configuration with zero effecting coupling between the qubit and resonator. Note that the coupler is also a transmon that shows a higher sensitivity to the magnetic flux than regular qubits. We next apply a flux pulse that lowers the coupler frequency, turning on the effective coupling between the resonator and the qubit. Depending on the flux pulse shape, the state collects conditional phase φand possibly experiences population oscillations between computational and noncomputational states, as a function of the time spent at the gate-operation frequency. We optimize the pulse amplitude and duration such that after the flux pulse the CZ gate fidelity is maximized. Details of the gate theory can be found in Appendix Gand the considered device parameters in Table II. 043089-8 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) FIG. 7. (a) CZ gate error landscape averaged over random initial states. Contours with a low error are highlighted with a dashed line. (b) iSWAP gate error landscape obtained by averaging over a number of random initial states in the zeroand one-excitation manifolds. Both plots are produced using system parameters shown in Table II. In Fig. 7(a), we operate our CZ gate by tuning the coupler frequency using a flattop Gaussian shaped flux pulse. The width of our Gaussian filter was fixed at 3 ns. Applying such a flux pulse to coupler results in a coupler frequency shift by ωshift cfrom the idling configuration. Then by appropriately tuning ωshift cand the gate time τ, one locates the optimal pulse configuration that minimizes the CZ(π) gate error εCZ =1− (tr√ρσ√ρ)2, where σis the target density matrix obtained after propagating some initial state ||with the ideal unitary of Eq. (12) and ρthe final density matrix obtained after propagating ||with the Lindbladian corresponding to our system defined in Eq. (G1). For our device parameters, the maximal decoherence limited CZ gate error averaged over a number of random initial states is 1.6×10−3. Note that the system parameters in Table II were chosen such that they allow for the possibility to find a good idling configuration, where the residual CZ interaction vanishes before the gate operation. In our simulations, we have included environmental noise, such as amplitude damping and pure dephasing and treated them using a Lindblad master equation solver in QUTIP [60,61]. 2. iSWAP gate Just as the CZ gate, the iSWAP gate can be natively realized in superconducting quantum computing architecture [37]. With our device, we can perform high-fidelity iSWAP gates between zeroand single-excitation computational states. The two-photon state |1r⊗|1, where |1rdenotes the first excited state of the resonator, must be excluded because it resonantly interacts with the state |2r⊗|0inducing a population exchange between the states. Hence the resulting operation in this subspace does not match the action of the targeted iSWAP operation. The capacitive coupling between the elements of the electrical circuit shown in Fig. 6(b) gives rise to an effective XY interaction between the qubit and resonator under the rotating wave approximation. Such an interaction conserves excitation number. With only the qubit or resonator (or neither) initially populated, we stay within the single excitation subspace of the joint system, thereby minimizing leakage of quantum population into the higher excited states of the resonator. The XY interaction can be turned on by first tuning the qubit in resonance with the resonator, and then applying a flux-pulse to the coupler to turn on the coupling, similar to the CZ gate operation. Figure 7(b) shows iSWAP gate error landscape for the same device parameters (given in Table II). The optimal average iSWAP gate error εiSWAP obtained for our device is 1.7×10−3. This result is obtained by averaging over a number of random initial states within the zeroand one-excitation manifolds. The results of our two-qubit-gate simulations demonstrate that our star architecture supports operating gates with similar fidelities as regular transmon qubits coupled together. The increased local connectivity of the device reduces the need for SWAP gates to simulate the nanoscale NMR problem (and others with a similar structure) and consequently in the end improves simulation fidelities. B. Reset The hyperpolarization protocol described in Sec. II needs regular re-initializations of the state of the NV center. The Co-Design hardware for simulating the protocol must therefore support this operation within qubit lifetimes. This is a hardware challenge, but one with solutions in sight. In particular, the quantum circuit refrigerator (QCR) has been used to perform the reset in tens of nanoseconds [40–43], which is a similar timescale to gate operations. The advantage of using a QCR for the reset is the possibility to reset the central resonator directly, without the need transfer the resonator population back to the central qubit using an iSWAP gate. Alternatively, a fast reset is possible through applying a flux drive to a qubit to SWAP its state with its measurement line [62]. This scheme has the advantage of not requiring any additional hardware not already present on the chip, but comes with a small cost in the circuit depth, as the state of the resonator must be transported using an iSWAP gate into the 043089-9 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) ( a )( b ) FIG. 11. (a) Coefficients vectors of the first qubit  A1, ωc 1before the rotation, with projection over the three axis and (b) coefficient vectors of the first qubit  Arot 1, ωc,rot 1after the rotation, being  Arot 1in the Zaxis. the continuous-driving case, due to the two different timedependent processes involved in the Trotter decomposition: the free dynamics of the spins and the sequence of pulses. The most crucial point to be aware of is the interplay between Trotter steps and interpulse spacing. The number of interpulse evolutions, i.e. number of pulses minus one, bounds from below the minimum number of Trotter steps for the simulation. Clearly, at least one Trotter step is needed for each interpulse evolution. Taking this interplay into account, the most straightforward setup is to choose a frequency which will determine the spacing of the pulse sequence, and to identify each interpulse evolution with a single Trotter step. If the achieved precision is not high enough, more Trotter steps can be added for each interpulse evolution. Each π-pulse itself is simply implemented as an Xor Ygate on the qubit representing the NV center. The OU-distributed Rabi frequency fluctuations present in nanoscale NMR systems are then simulated by overand under-rotations of the Xand Ygates. APPENDIX E: ROTATIONAL OPTIMIZATION In principle, we had a Hamiltonian with terms of the type ZX,ZY,and ZZ for the case of no internuclear interactions. However, we can rotate the basis so the Hamiltonian loses the ZX and ZY terms, allowing to reduce the number of TQGs. To make up for this rotation, we need to introduce different constants  Arot ifor the problem and rotate the vector state we obtain at the end before measuring it. The rotations that we will consider are only one-qubit rotations on nuclei qubits and we are applying this just to the case with no internuclear interactions. Therefore we can consider the effect of this rotation on only one qubit representing an arbitrary nucleus. We will exemplify this procedure using nucleus 1. If we want to obtain the mean value of σzacting on the nucleus: σz=Tr(ρ(tf)σz)=Tr(U(0,tf)ρ(0)U†(0,tf)σz),(E1) where U(0,tf) represents the evolution operator from t=0 to t=tf. The density matrix ρ(0) contains the state of the NV center (which is in the |+ state at t=0) and nucleus 1, i.e., ρ(0) = |++|NV ⊗11 2. Our intention is to obtain an expression of this mean value in terms of the rotated evolution operators and later, we will find the appropriate rotation to be performed. Then, taking into account that the trace is invariant under a rotation R=1NV ⊗R1, we get: σz=Tr(RU (0,tf)ρ(0)U†(0,tf)σzR†) =Tr(RU (0,tf)R†Rρ(0)R†RU †(0,tf)R†RσzR†).(E2) This can be expressed as: σz=Tr(Urot(0,tf)ρrot(0)U† rot(0,tf)RσzR†).(E3) The density matrix of the nucleus is the identity. Thus any rotation on nuclei qubits leaves the density matrix unaffected, leading to: σz=Tr(Urot(0,tf)ρ(0)U† rot(0,tf)RσzR†).(E4) Then we need to rotate the system previous to the measurement. By using the invariance of the trace under cyclic permutations, we get: σz=Tr(R†Urot(0,tf)ρ(0)U† rot(0,tf)Rσz),(E5) which is equivalent to introducing a counter-rotation in the circuit before measurement. Now let us focus on the specific rotation we have to implement. Since the constants multiplying the Pauli matrices in the Hamiltonian are  A1 2and  ωc 1= A1 2−γcBz ez(for nucleus 1), we can rotate the basis to obtain a representation in which the vectors have only zcomponent for  A1and thus, XZ and YZ terms are removed. The vectors before and after the needed rotation can be seen in Fig. 11. To compute the new vectors (and thus the new coefficients for the gates of our algorithm), we can use Rodrigues’ rotation formula to rotate a vector  van angle θaround a unitary 043089-16 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) TABLE III. Gate count for one Trotter step and for one cycle for different topologies with and without internuclear interactions. All-To-All Star topology Square grid Nnonint TQG n−1n−14n−4 Nnonint SQG 5 2n+25 2n+221 2n−47 2 Nint TQG 3 2n2−3 2n3 2n2+3 2n−63n2−6n+3 Nint SQG 4n2−9 2n+7 24n2+7 2n−25 28n2−33 2n+11 2 axis ˆ k:  vrot = vcos θ+(ˆ k× v)sinθ+ˆ k(ˆ k· v)(1 −cos θ),(E6) being in our case, θ=arccos (Az 1/| A1|) and ˆ k= (cos(φ),sin(φ),0), with φ=−π 2+φxy =−π 2+arctan (Ay 1/Ax 1). For implementing the counter-rotation of this in the quantum circuit, we use: R† 1=eiθ 2(cos(φ)X−sin(φ)Y).(E7) APPENDIX F: SWAP ROUTING Our qubit routing method consists of mapping the square grid to a linear chain with qubits labeled from 0 to n. Then, in the simplified case of no internuclear interactions, the optimal SWAP method for the one-to-all interaction case on a linear chain can be used. For a single NV center the protocol goes as follows. (1) Initialize the state of the NV center in the second qubit. (2) Perform interactions with the first and third qubits. (3) SWAP the NV center qubit to the right. (4) Perform interaction with right qubit. (5) Repeat steps 3-4 until all interactions have been achieved. The pattern is seen in Fig. 3(a) denoted by the intense blue arrows. With internuclear interactions we need to perform a swap pattern that enables all-to-all interactions. The so-called odd-even mapping in Fig. 3(a) is an efficient one [69] represented by green arrows in Fig. 3(a). This consists of swapping first all the even qubits with their right neighbors and then swapping all the odd qubits with their right neighbors. This way, we will obtain all-to-all interactions with 1 2(n−1)(n−2) SWAP gates and a total TQG depth of 6n.A summary of the TQG counts is shown in Table III. To motivate the creation of a chip with a star topology and the use of an alternative linearized SWAP routing for a square grid instead of standard numerical approaches, a comparison between all the cases is provided in Fig. 12. A reduction in the number of SWAPs can be noticed for both the linear chain approach and the star-topology chip against standard numerical approaches for a square grid. APPENDIX G: QUBIT-RESONATOR GATE THEORY In the following discussion, we consider gate operation between the resonator and one of the qubits, and neglect any effects that arise from the interactions with spectator qubits and other resonator modes. The time dynamics in such a system are determined by the Hamiltonian: H=H0+Hrc +Hqc +Hrq,(G1) where the uncoupled part of the total Hamiltonian H0=Hr+ Hc+Hqis: Hr=¯hωrb† rbr, Hc=¯hωcb† cbc+¯h 2αcb† cb† cbcbc, Hq=¯hωqb† qbq+¯h 2αqb† qb† qbqbq, (G2) FIG. 12. (a) Comparison of the required number of SWAPs for simulating the proposed system with no internuclear interactions for each Trotter step. Numerical approaches from references are applied to a square grid. (b) Equivalent comparison with internuclear interactions. Zulehner et al. and Saeedi et al. do not improve the linear chain approach for few qubits and are intractable for larger numbers of qubits and thus are not displayed. 043089-17 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) where bλand ωλare the annihilation operator and fundamental frequency for the mode λ={r,c,q}, respectively, and αγ is the anharmonicity of the mode γ={q,c}. The interaction component of the Hamiltonian is: Hλμ =−¯hgλμ(b† λ−bλ)(b† μ−bμ),(G3) where λμ ={rc,qc,rq}, and gλμ denote resonator-coupler, qubit-coupler and resonator-qubit coupling frequencies. With the Hamiltonian of Eq. (G1), we are now in a position to perform simulations of two-qubit gates by propagating a suitably chosen initial state. Before the gate operation, we choose the idling frequencies for the qubit, resonator, and the coupler such that the CZ coupling rate ζis minimized. This CZ coupling rate is defined as: ζ=ω101 −ω100 −ω001 +ω000,(G4) where ωnr0nqcorresponds to the eigenenergy of Hamiltonian in Eq. (G1) with nrexcitations in resonator and nqexcitations in qubit with coupler being in the ground state. The point of minimal |ζ|is also known as the idling configuration, which we found to be at [ωr,ω c,ω q]/(2π)=[4.30,6.14,4.47] GHz for the parameters given in Table II. The CZ gate is operated by sending a flux pulse that modifies the coupler frequency ωc, which then in the coupled basis modifies the frequencies ω101,ω 100,ω 001 and ω000. This makes ζnonzero, so the system collects a CZ phase. [1] R. Feynman, Simulating physics with computers, in Feynman and Computation (CRC Press, 2018), pp. 133–153. [2] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996). [3] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018). [4] K. Bharti, A. Cervera-Lierta, T. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. Kottmann, T. Menke et al., Noisy intermediate-scale quantum (NISQ) algorithms, Rev. Mod. Phys. 94, 015004 (2022). [5]M.J.S.Smith,Application-Specific Integrated Circuits,Vol.7 (Addison-Wesley Reading, MA, 1997). [6] K.-C. Hsu and H.-W. Tseng, GPTPU: Accelerating applications using edge tensor processing units, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Association for Computing Machinery, New York, NY, 2021), pp. 1–4, doi: 10.1145/3458817.3476177. [7] T. Lu, T. Marin, Y. Zhuo, Y.-F. Chen, and C. Ma, Accelerating MRI reconstruction on TPUs, in 2020 IEEE High Performance Extreme Computing Conference (HPEC) (IEEE, 2020), pp. 1–9, doi: 10.1109/HPEC43674.2020.9286192. [8] T. Staudacher, F. Shi, S. Pezzagna, J. Meijer, J. Du, C. A. Meriles, F. Reinhard, and J. Wrachtrup, Nuclear magnetic resonance spectroscopy on a (5-Nanometer)3 sample volume, Science 339, 561 (2013). [9] M. Nielsen and I. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, 2000). [10] H. Hogben, M. Krzystyniak, G. Charnock, P. Hore, and I. Kuprov, Spinach a software library for simulation of spin dynamics in large spin systems, J. Magn. Reson. 208, 179 (2011). [11] K. Seetharam, D. Biswas, C. Noel, A. Risinger, D. Zhu, O. Katz, S. Chattopadhyay, M. Cetina, C. Monroe, E. Demler et al., Digital quantum simulation of NMR experiments, arXiv:2109.13298. [12] T. E. O’Brien, L. B. Ioffe, Y. Su, D. Fushman, H. Neven, R. Babbush, and V. Smelyanskiy, Quantum computation of molecular structure using data from challenging-toclassically-simulate nuclear magnetic resonance experiments, arXiv:2109.02163. [13] D. Sels, H. Dashti, S. Mora, O. Demler, and E. Demler, Quantum approximate Bayesian computation for NMR model inference, Nat. Mach. Intell 2, 396 (2020). [14] D. Sels and E. Demler, Quantum generative model for sampling many-body spectral functions, Phys.Rev.B103, 014301 (2021). [15] M. H. Levitt, Spin Dynamics: Basics of Nuclear Magnetic Resonance (Wiley, 2013). [16] M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. L. Hollenberg, The nitrogen-vacancy colour centre in diamond, Phys. Rep. 528, 1 (2013). [17] M. H. Abobeih, J. Randall, C. E. Bradley, H. P. Bartling, M. A. Bakker, M. J. Degen, M. Markham, D. J. Twitchen, and T. H. Taminiau, Atomic-scale imaging of a 27-nuclearspin cluster using a quantum sensor, Nature (London) 576, 411 (2019). [18] C. Munuera-Javaloy, R. Puebla, and J. Casanova, Dynamical decoupling methods in nanoscale NMR, Europhys. Lett. 134, 30001 (2021). [19] C. Munuera-Javaloy, R. Puebla, B. D’Anjou, M. B. Plenio, and J. Casanova, Detection of molecular transitions with nitrogenvacancy centers and electron-spin labels, arXiv:2110.14255. [20] A. Ajoy, K. Liu, R. Nazaryan, X. Lv, P. Zangara, B. Safvati, G. Wang, D. Arnold, G. Li, A. Lin et al., Orientation-independent room temperature optical 13c hyperpolarization in powdered diamond, Sci. Adv. 4, eaar5492 (2018). [21] S. Hartmann and E. Hahn, Nuclear double resonance in the rotating frame, Phys. Rev. 128, 2042 (1962). [22] H. Carr and E. Purcell, Effects of diffusion on free precession in nuclear magnetic resonance experiments, Phys. Rev. 94, 630 (1954). [23] S. Meiboom and D. Gill, Modified spin-echo method for measuring nuclear relaxation times, Rev. Sci. Instrum. 29, 688 (1958). [24] A. Maudsley, Modified Carr-Purcell-Meiboom-Gill sequence for NMR Fourier imaging applications, J. Magn. Reson. (1969) 69, 488 (1986). [25] T. Gullion, D. Baker, and M. Conradi, New, compensated CarrPurcell sequences, J. Magn. Reson. (1969) 89, 479 (1990). [26] J. Casanova, Z. Y. Wang, and M. B. Plenio, Noise-Resilient Quantum Computing with a Nitrogen-Vacancy Center and Nuclear Spins, Phys.Rev.Lett.117, 130502 (2016). [27] G. Uhlenbeck and L. Ornstein, On the theory of the Brownian motion, Phys. Rev. 36, 823 (1930). [28] J. Cai, B. Naydenov, R. Pfeiffer, L. McGuinness, K. Jahnke, F. Jelezko, M. Plenio, and A. Retzker, Robust dynamical 043089-18 CO-DESIGN QUANTUM SIMULATION OF NANOSCALE NMR PHYSICAL REVIEW RESEARCH 4, 043089 (2022) decoupling with concatenated continuous driving, New J. Phys. 14, 113023 (2012). [29] D. T. Gillespie, Exact numerical simulation of the OrnsteinUhlenbeck process and its integral, Phys. Rev. E 54, 2084 (1996). [30] A. Ajoy, B. Safvati, R. Nazaryan et al., Hyperpolarized relaxometry based nuclear T1 noise spectroscopy in diamond, Nat. Commun. 10, 5160 (2019). [31] M. Suzuki, Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems, Commun. Math. Phys. 51, 183 (1976). [32] E. Campbell, Random Compiler for Fast Hamiltonian Simulation, Phys.Rev.Lett.123, 070503 (2019). [33] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019). [34] K. Bharti and T. Haug, Quantum-assisted simulator, Phys. Rev. A104, 042418 (2021). [35] E. Younis, K. Sen, K. Yelick, and C. Iancu, Qfast: Conflating search and numerical optimization for scalable quantum circuit synthesis, arXiv:2103.07093. [36] N. Hatano and M. Suzuki, Finding Exponential Product FormulasofHigherOrders(Springer, Berlin, Heidelberg, 2005), pp. 37–68. [37] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Appl. Phys. Rev. 6, 021318 (2019). [38] D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Efficient Z gates for quantum computing, Phys. Rev. A96, 022330 (2017). [39] F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Tunable Coupling Scheme for Implementing High-Fidelity Two-Qubit Gates, Phys. Rev. Appl. 10, 054062 (2018). [40] K. Y. Tan, M. Partanen, R. E. Lake, J. Govenius, S. Masuda, and M. Möttönen, Quantum-circuit refrigerator, Nat. Commun. 8, 15189 (2017). [41] M. Silveri, H. Grabert, S. Masuda, K. Y. Tan, and M. Möttönen, Theory of quantum-circuit refrigeration by photon-assisted electron tunneling, Phys. Rev. B 96, 094524 (2017). [42] H. Hsu, M. Silveri, A. Gunyhó, J. Goetz, G. Catelani, and M. Möttönen, Tunable refrigerator for nonlinear quantum electric circuits, Phys.Rev.B101, 235422 (2020). [43] V. Sevriuk, K. Y. Tan, E. Hyyppä, M. Silveri, M. Partanen, M. Jenei, S. Masuda, J. Goetz, V. Vesterinen, L. Grönberg et al., Fast control of dissipation in a superconducting resonator, Appl. Phys. Lett. 115, 082601 (2019). [44] M. Celio, New Method to Calculate the Muon Polarization Function, Phys. Rev. Lett. 56, 2720 (1986). [45] S. McArdle, Learning from physics experiments with quantum computers: Applications in muon spectroscopy, PRX Quantum 2, 020349 (2021). [46] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Surface codes: Towards practical large-scale quantum computation, Phys.Rev.A86, 032324 (2012). [47] Y. Hirata, M. Nakanishi, S. Yamashita, and Y. Nakashima, An efficient method to convert arbitrary quantum circuits to ones on a linear nearest neighbor architecture, in 2009 Third International Conference on Quantum, Nano and Micro Technologies, Cancun, Mexico (IEEE, 2009), pp. 26–33, doi: 10.1109/ICQNM.2009.25. [48] G. Li, Y. Ding, and Y. Xie, Tackling the qubit mapping problem for NISQ-era quantum devices, in Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems (ACM, New York, NY, USA, 2019), pp. 1001–1014, doi: 10.1145/3297858.3304023. [49] M. Saeedi, R. Wille, and R. Drechsler, Synthesis of quantum circuits for linear nearest neighbor architectures, Quant. Info. Proc. 10, 355 (2011). [50] A. Zulehner, A. Paler, and R. Wille, An efficient methodology for mapping quantum circuits to the IBM QX Architectures, in IEEE transactions on computer-aided design of integrated circuits and systems (IEEE, 2019), Vol. 38, pp. 1226–1236, doi: 10.1109/TCAD.2018.2846658. [51] E. Magesan and J. M. Gambetta, Effective Hamiltonian models of the cross-resonance gate, Phys.Rev.A101, 052308 (2020). [52] N. Lacroix, C. Hellings, C. K. Andersen, A. D. Paolo, A. Remm, S. Lazar, S. Krinner, G. J. Norris, M. Gabureac, J. Heinsoo, A. Blais, C. Eichler, and A. Wallraff, Improving the performance of deep quantum optimization algorithms with continuous gate sets, PRX Quantum 1, 020304 (2020). [53] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv:1411.4028. [54] M. Mariantoni, F. Deppe, A. Marx, R. Gross, F. K. Wilhelm, and E. Solano, Two-resonator circuit quantum electrodynamics: A superconducting quantum switch, Phys.Rev.B78, 104508 (2008). [55] B. Foxen, C. Neill, A. Dunsworth, P. Roushan, B. Chiaro, A. Megrant, J. Kelly, Z. Chen, K. Satzinger, R. Barends et al., Demonstrating a Continuous Set of Two-qubit Gates for Near-term Quantum Algorithms, Phys. Rev. Lett. 125, 120504 (2020). [56] C. Song, K. Xu, W. Liu, C. P. Yang, S. B. Zheng, H. Deng, Q. Xie, K. Huang, Q. Guo, L. Zhang, P. Zhang, D. Xu, D. Zheng, X. Zhu, H. Wang, Y. A. Chen, C. Y. Lu, S. Han, and J. W. Pan, 10-Qubit Entanglement and Parallel Logic Operations with a Superconducting Circuit, Phys.Rev.Lett.119, 180511 (2017). [57] C. Song, K. Xu, H. Li, Y. Zhang, X. Zhang, W. Liu, Q. Guo, Z. Wang, W. Ren, J. Hao et al., Generation of multicomponent atomic Schrödinger cat states of up to 20 qubits, Science 365, 574 (2019). [58] S. Hazra, A. Bhattacharjee, M. Chand, K. V. Salunkhe, S. Gopalakrishnan, M. P. Patankar, and R. Vijay, Ring-ResonatorBased Coupling Architecture for Enhanced Connectivity in a Superconducting Multiqubit Network, Phys. Rev. Appl. 16, 024018 (2021). [59] J. Chu and F. Yan, Coupler-Assisted Controlled-Phase Gate with Enhanced Adiabaticity, Phys. Rev. Appl. 16, 054020 (2021). [60] J. Johansson, P. Nation, and F. Nori, Qutip: An open-source Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760 (2012). [61] J. Johansson, P. Nation, and F. Nori, Qutip 2: A Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013). [62] Y. Zhou, Z. Zhang, Z. Yin, S. Huai, X. Gu, X. Xu, J. Allcock, F. Liu, G. Xi, Q. Yu et al., Rapid and unconditional 043089-19 MANUEL G. ALGABA et al. PHYSICAL REVIEW RESEARCH 4, 043089 (2022) parametric reset protocol for tunable superconducting qubits, arXiv:2103.11315. [63] S. Endo, S. C. Benjamin, and Y. Li, Practical Quantum Error Mitigation for Near-Future Applications, Phys.Rev.X8, 031027 (2018). [64] Z. Cai, Multi-exponential error extrapolation and combining error mitigation techniques for NISQ applications, npj Quantum Inf. 7, 80 (2021). [65] M. Krebsbach, B. Trauzettel, and A. Calzona, Optimization of Richardson extrapolation for quantum error mitigation, arXiv:2201.08080. [66] A. M. Childs, A. Ostrander, and Y. Su, Faster quantum simulation by randomization, Quantum 3, 182 (2019). [67] C. F. Chen, H. Y. Huang, R. Kueng, and J. A. Tropp, Concentration for random product formulas, PRX Quantum 2, 040305 (2021). [68] Y. Ouyang, D. R. White, and E. T. Campbell, Compilation by stochastic Hamiltonian sparsification, Quantum 4, 235 (2020). [69] A. Cowtan, S. Dilkes, R. Duncan, A. Krajenbrink, W. Simmons, and S. Sivarajah, On the Qubit Routing Problem, in 14th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 135 (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 2019), pp. 51–532, doi: 10.4230/LIPIcs.TQC.2019.5. 043089-20