{'fraction_non_alphanumeric': 0.04160570160947751, 'fraction_numerical': 0.015245197526785293, 'mean_word_length': 4.288281542493448, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Second-order shaped pulses for solid-state quantum computation

We present the constructon and detailed analysis of highly-optimized self-refocusing pulse shapes for several rotation angles. We characterize the constructed pulses by the coefficients appearing in the Magnus expansion up to second order. This allows a semi-analytical analysis of the performance of the constructed shapes in sequences and composite pulses by computing the corresponding leading-order error operators. Higher orders can be analyzed with the numerical technique suggested by us previously. We illustrate the technique by analysing several composite pulses designed to protect against pulse amplitude errors, and on decoupling sequences for potentially long chains of qubits with on-site and nearest-neighbor couplings.


I. INTRODUCTION
The implementation of quantum algorithms using NMR on molecules in liquid 1 , solid 2 , and liquid crystals 3 has demonstrated in principle that pulse-based control methods can be useful for quantum information processing 4 (QIP). The technique has long been a staple in NMR spectroscopy, where complex molecules like proteins are analyzed with the help of long sequences of precisely designed radiofrequency pulses 5 . Related techniques for coherent manipulation of many-body quantum systems have emerged as an important tool in many areas of science and technology.
A useful quantum computer should contain hundreds, if not thousands of qubits. The only hope of scaling to such system sizes is with the help of multiple levels of quantum error correction (QEC). For this to work, the benefits due to each additional level of encoding should outweigh the corresponding overhead of additional errors. This leads to various threshold theorems 6,7 , estimating the maximum error rate for which such concatenated encoding can be beneficial. The corresponding thresholds are rather stringent, meaning that for scalability one needs very accurate elementary gates.
Even for relatively small n-body systems, the number of states scales exponentially with n, and the accuracy required for QIP is high. As demonstrated in several recent experiments in NMR QIP, required accuracy can be reached with the help of strongly-modulating pulses, where entire single-or few-qubit gates are designed numerically for a given molecule 8,9 [also see 10,11,12,13]. While the technique indeed offers unprecedentedly accurate, fast gates (which also helps to avoid relaxation), it obviously cannot be generalized to larger systems.
In contrast, the traditional pulse and sequence design rely on the Magnus (cumulant) expansion 5 . The expansion is done around the evolution in the applied controlling fields, while the chemical shifts 14,15,16 (resonantfrequency offsets) and inter-qubit couplings are treated perturbatively. The main advantage of the Magnus expansion is its locality. Namely, when local qubit cou-plings are dominant, the control fields accurate to a given order can be designed by analyzing relatively small clusters. The results remain exact independent of the system size, or even in the limit of infinite system. One can thus characterize pulse-based method for designing control fields as scalable to large system sizes. 17 .
A scheme to systematically construct high-order selfrefocusing pulses and pulse sequences was developed by the authors in Ref. 18. Specifically, we constructed "soft" NMR-style 14,19 second-order self-refocusing inversion (π) pulses and several high-order sequences based on such pulses for refocusing qubits arranged in spin chains with on-site chemical shifts and XXZ nearest-neighbor coupling. The main technical advance which enabled the calculations 18 was the efficient numerical algorithm for computing high-order terms of the Magnus expansion. The algorithm is based on the usual time-dependent perturbation theory; the direct computation of multiple integrals entering higher-order cumulants would be totally impossible.
In this paper we present highly-optimized selfrefocusing pulse shapes for rotation angles φ 0 other than 180 • . For such pulses, we extend the results of Ref. 20 and construct the analytical expansion of the evolution operator for an arbitrary coupled qubit. While the expansion is more complicated than that for the inversion pulses with φ 0 = 180 • , to second order, it is still characterized by only three coefficients, two of which we suppress by pulse shaping. This allows us to compute the error operators associated with a given control sequence semi-analytically, by evaluating the leading order terms in the corresponding products of the evolution operators. We illustrate the technique on several newly-constructed decoupling sequences for a chain of qubits with on-site and nearest-neighbor couplings, as well as with the composite pulses protecting against amplitude errors.

A. Dynamical Decoupling
In principle, the simplest type of control pulses consist of short, intense bursts of coherent, resonant electromagnetic radiation, popularly known as "hard" or "bangbang" pulses. In this limit, for the duration τ of the pulse one can totally ignore all other couplings of the qubit(s). Then, a pulse sequence can be viewed as a series of free evolution intervals [unitary evolution operators U f (t) = exp(−itH S ), where H S is the system Hamiltonian] intercalated with pulse operators. For example, with single-qubit control Hamiltonian, where σ µ a , µ = x, y, z are the Pauli matrices for a-th qubit and V µ a (t) are the corresponding control fields (or, more precisely, the envelopes of the control fields applied at the resonant frequency of the corresponding qubit), the pulse operator is the product of those for individual qubits, P = a P a (φ (a) 0 ,n a ), Heren is the unit vector that determines the spin rotation axis and φ 0 is the corresponding angle, The corresponding pulse algebra is straightforward. For example, for a single spin with the chemical shift Hamiltonian, the sequence of two equally-spaced inversion pulses (φ 0 = ±π) in x-direction is equivalent to a unitary, which is, of course, the formal identity behind the wellknown spin-echo sequence 21 .
In reality, the pulse duration τ is always finite. Thus, during the pulse application, the rotation actually happens around the axis determined by the sum of the system and control Hamiltonians; one gets finite-pulseduration errors. Generically, such errors scale linearly with pulse duration. These errors are especially dangerous in systems where one cannot reduce the pulse duration indefinitely, e.g., because of the need for spectral addressing in homonuclear NMR, or in order to avoid exciting levels outside the qubit state in Josephson phase qubits 22 . Yet, even in systems with optical qubit coupling where τ could be in the sub-picosecond range, the finite-pulse-duration errors can be significant if one targets the accuracy necessary for achieving the scalability thresholds 23 .
The finite-pulse-duration errors can be significantly reduced by suppressing the leading-order error operators. The errors would scale with a higher power of the pulse duration, which makes them much more manageable. This can be achieved either by designing special sequences (which typically doubles or quadruples the number of pulses), or by using specially designed selfrefocusing pulse shapes 14 . Typically, the latter strategy is more efficient 18 ; besides, the self-refocusing pulses can often be used as drop-in replacements for corresponding δ-pulses 18,20 .

B. Model
We consider the following simplified Hamiltonian with the first (main) term due to individual control fields, where σ µ n , µ = x, y, z, are the usual Pauli matrices for the n-th qubit (spin) of the 1D chain. The other terms include the native Hamiltonian of the system describing the constant qubit couplings and the interactions between the qubits, and the coupling with the oscillator thermal bath, In Eq. (8), A µ n ≡ A µ n (p i , q i ) account for the possibility of a direct coupling of the controlling fields V µ n with the bath variables q i , p i , while B µ n ≡ B µ n (p i , q i ) describe the usual coupling of the spins with the oscillator bath. Already in the linear response approximation, the bath couplings (8) produce a frequency-dependent renormalization of the control Hamiltonian H C (t) [Eq. (6)], as well as the thermal bath heating via the dissipative part of the corresponding response function. Both effects become more of a problem with increased spectral width of the controlling signals V µ n . In this work we do not specify the explicit form of the coupling H V (t). Instead, we minimize the spectral width of the constructed pulses.
While the Hamiltonian (7) is a generic spin Hamiltonian, we will also consider specifically the Hamiltonian of XXZ model with additional on-site fields, n,n ′ J z n,n ′ σ z n σ z n ′ + J ⊥ n,n ′ (σ x n σ x n ′ + σ y n σ y n ′ )

C. Magnus expansion
In a qubit-only system with the Hamiltonian the effect of the applied fields is fully described by the evolution operator U (t), where T is the Dyson time ordering operator. For pulses with finite width, any desired unitary transformation can only be implemented approximately. A widely used framework to design pulses to effect a desired unitary transformation (or, equivalently, remove the effect of undesired terms in the Hamiltonian) is the Magnus 24 expansion and the average Hamiltonian theory 25,26 . The expansion is done with respect to the evolution due to the control fields alone, by defining the system Hamiltonian in the interaction representation (the "rotating-frame Hamiltonian"), For a periodic control field, H 1 (t + τ c ) = H 1 (t), such that the zeroth-order driven evolution is also periodic, U 0 (t + τ c ) = U 0 (t), one has the following expansion in powers of τ c : Generally, the termH (k−1) contains a k-fold integration of the commutators of the rotating-frame Hamiltoniañ H i ≡H S (t i ) at different time moments t i and scales as H (k) τ c ∝ τ c H S k . Note that for small enough τ c , the expansion parameter remains small even for long evolution time.
The Magnus expansion thus offers a basis for constructing successful approximations towards the desired unitary evolution. With the simpler problem of decoupling, the goal is to have no evolution. A K-th order refocusing sequence can be defined as that where there is no evolution to K-th order, that is, Respectively, at time nτ c , the error in the unitary evolution operator would scale as U (nτ c ) − 1 1 ∝ n τ c H S K+1 , and the corresponding fidelity F (t) differs from unity by A crucial advantage of the cumulant expansion is that the cumulants do not contain the disconnected terms arising from different parts of the system (cluster theorem 27,28 ). For an arbitrary lattice model of the form (9), with bonds representing the qubit interactions, the terms contributing to k-th order can be represented graphically as connected clusters involving up to k lattice bonds; for a chain of qubit such clusters cannot have more than n = k + 1 vertices. Thus, to obtain the exact form of the expansion up to and including K-th order, one needs to analyze all distinct chain clusters with up to K + 1 vertices.

D. Time dependent perturbation theory
While the Magnus expansion is conceptually straightforward, it is cumbersome to implement and, most importantly, the repeated integrations are very expensive computationally already at the second order, see Eq. (17). An alternative strategy for evaluating high-order terms of the Magnus expansion was suggested by the present authors in Ref. 18. Instead of working with the cumulants, the technique is based on the time-dependent perturbation theory (TDPT). The expansion is done around the non-perturbed evolution due to control fields alone, see Eq. (12). However, for actual computations, it is more convenient to use the differential equatioṅ The slow evolution operator obeys the equatioṅ which can be iterated to construct the standard expansion R(t) = I + R 1 (t) + R 2 (t) + . . . in powers of (t H S ), The successive terms can be evaluated by solving, at each step, a set of coupled first order ODE's simultaneously. For a finite system of n qubits and a given maximum order K of the expansion, one needs to solve Eqs. (22) and (25) with 1 ≤ k ≤ K. These are total of (K + 1) coupled systems of first order ordinary differential equations for the 2 n × 2 n matrices U 0 , R 1 , R 2 , . . . , R K , and can be integrated efficiently. Computationally, this is a much less challenging job than that of evaluating repeated integrals (17), (18), and higher order terms. For a given system, solving the full set of equations (24) is simpler by a factor of at least (K + 1). However, it is the analysis of the perturbative expansion that is the key for achieving the scalability of the results. Given the set of computed R k (τ c ) ≡ R k , the standard Magnus expansion can be readily obtained by taking the logarithm of the series, −iH (2) Obviously, the order-K universal self-refocusing condition (20) is formally equivalent to The matrices R k in the latter condition are much easier to evaluate numerically using Eqs. (22), (25). Importantly, the benefits of the cluster theorem are retained: to K-th order only clusters with up to K + 1 vertices need to be analyzed.

III. PULSE DESIGN AND ANALYSIS
A. Pulse design using TDPT The shapes of NMR-style one-dimensional pulses 14 , self-refocusing to a given order, can be found by analyzing the single-spin dynamics with the system Hamiltonian (3) and the control Hamiltonian Specifically, we encoded the trial pulse shapes in terms of their Fourier coefficients, where τ is the pulse duration, Ω ≡ 2π/τ is the corresponding angular frequency, and φ 0 is the requested rotation angle of the pulse. Note that the form (31) guarantees the symmetry of the pulse, V (τ − t) = V (t). In addition, in order to reduce the spectral width of the control fields, we also constrained a certain number of derivatives of the function (31) to vanish at t = 0 and We implemented the computational algorithm described in the previous section using the standard fourthorder Runge-Kutta algorithm for solving coupled differential equations (22), (25), and the GSL library 29 for matrix operations. The coefficient optimization was done using a combination of simulated annealing and the steepest descent method. The target function for single-pulse optimization included the sum of the magnitudes squared of the matrix elements of the matrices R k ≡ R k (τ ), k = 1, . . . , K, as well as the weighed sum of the squares of the coefficients A m , The second sum serves to provide some suppression of the higher Fourier harmonics of the pulse. In our simulations, the minimization was considered as having converged only after the first term evaluated to zero with numerical precision (typically, eight digits or better). For such a minimum to exist, the coefficient ǫ in Eq. (32) should be sufficiently small (we used ǫ = 10 −4 ). For given pulse order K and the given number of additional constraints L, there is a minimum number of harmonics M min (K, L) necessary for convergence. However, we found that the shapes obtained with M = M min (K, L) tend be over constrained and simply do not look nice. Our solution was to add one or two additional Fourier harmonics by increasing M .

B. Pulse shapes
Previously, in Ref. 18, we gave the coefficients of the first-order self-refocusing (K = 1) inversion (φ 0 = π) pulse shapes S L , as well as the second-order (K = 2) inversion shapes Q L , L = 1, 2. Here L is the parameter that indicates the number of constraints at the ends of the interval: the value of the function and its derivatives up to (2L − 1)st vanish at the ends of the interval [note that all odd derivatives are suppressed automatically due to the symmetry of the function, see Eq. (31)].
In this work, we extend the list of constructed pulses to rotation angles φ 0 = 10 • , 20 • , . . . , 180 • . In Table I, we list the coefficients for the pulses used in the simulations. The coefficients for all of the constructed pulses are available upon request.   (31). Shapes SL(φ0) and QL(φ0) are the pulse shapes for rotation angle φ0, respectively first (K = 1) and second (K = 2) order for the Hamiltonian (3). These shapes have 2L derivatives vanishing at the ends of the interval.

C. Pulse shape analysis
The pulse shapes Q L (φ 0 ), S L (φ 0 ) are constructed as first-or second-order self-refocusing pulses for the chemical shift Hamiltonian, Eq. (3). We would like, however, to have a universally good pulse shape that would work in most settings. To analyze the performance of the constructed pulses in most general circumstances, we construct the Magnus expansion of the evolution operator for the most general system Hamiltonian, where A µ are the operators responsible for coupling with the outside worlds, and A 0 is the external Hamiltonian. The analysis of the inversion pulses with φ 0 = π appeared previously in Ref. 20; here we extend it to φ 0 = π.

Driven evolution
The control Hamiltonian (30) alone [to zeroth order in H S ] produces the following unitary evolution operator [cf. Eq. (12)] When acting on the spin operators, this is just a rotation, . Consequently, the system Hamiltonian in the interaction representation has the form

Leading-order average Hamiltonian
The zeroth order average Hamiltonian (17) is just the average of Eq. (35) over the pulse duration. We assume is the overall notation angle. It is convenient to introduce the symmetrized rotation angle, such that ϕ(τ p − t) = −ϕ(τ p ). Then, the average of the sine over the pulse duration vanishes, sin ϕ = 0. This implies that the averages of the cosine and sine of the original rotation angle are where If we denote then the zeroth order average Hamiltonian for a onedimensional pulse becomes For the special case of the chemical-shift Hamiltonian (3), we have A 0 = A x = A y = 0, A z = ∆/2, and Eq. (41) givesH As seen, the spectral weight is almost entirely confined to ω < 5ω0.
Clearly, the 1st-order self-refocusing condition corresponds to υ = 0. For such pulses the full zeroth-order average Hamiltonian is given just by the two first terms in Eq. (42).

1st-order average Hamiltonian
The 1st-order average Hamiltonian (17) is given by a double integral of the commutator of the system Hamiltonian in the interaction representation (35) evaluated at two different times. We note that every term in Eq. (35) can be classified as either time-independent [proportional to e(t) ≡ 1], proportional to c(t) ≡ cos φ(t), or to s(t) ≡ sin φ(t). Therefore, most generally, the secondorder terms in the evolution operator can contain the following nine integrals, where we used the notation and 1 ≡ e(t) or 1 ′ ≡ e(t ′ ) indicate an identity factor at the corresponding position of the average. However, because of the commutator structure in Eq. (17), only the following antisymmetric combinations appear in the expression for the corresponding term in the average Hamiltonian theory, H (1) , where In fact, the coefficients ζ C and ζ S can be reduced further, Thus, to second order, the average Hamiltonian of a symmetric angle-φ 0 one-dimensional pulse is determined by only three dimensionless coefficients, υ, α, and ζ, see Eqs. (40), (46), and (48). These coefficients contain all the relevant information about the shape of the pulse. An explicit calculation of the 1st-order average Hamiltonian gives For the Hamiltonian (3), the terms with ζ disappear, and we have, simply Thus, the second-order self-refocusing pulses have both υ = 0 and α = 0.
The actual parameters for the pulses with φ 0 = π/2, π, and 2π are listed in TAB. II.  II: Parameters of several common pulse shapes. The first line represents the "hard" δ-function pulse applied at the center of the interval of duration τ , G001 denotes the Gaussian pulse with the width 0.01τp, while Sn and Qn denote the 1st and 2nd-order self-refocusing pulses from Tab. I.

IV. OPEN SYSTEMS
In this work we concentrate on the performance of high-order pulses and pulse sequences in closed quantum systems. However, it turns out that such sequences also remain efficient in open systems, in the presence of lowfrequency bath modes 20,30 .
The analysis is done in general form with the help of an assumption that the bath couplings have the same form as the existing terms in the system Hamiltonian (7), which are assumed to be suppressed to order K = 1 or K = 2. The bath modes are assumed to be lowfrequency; in addition to the expansion in powers of the corresponding couplings, one needs a low frequency expansion in powers of the adiabaticity parameter τ c /τ 0 , where τ c is the decoupling cycle duration and τ 0 is the bath correlation time.
With K = 1 decoupling, the effect in the open system is a suppression of direct decay (T 1 ) processes, as well as the reduction of the dephasing rate (T 2 ) by the factor of order of the adiabaticity parameter τ c /τ 0 . The former result can be understood by analyzing the spec-tral properties of the driven system 31,32 , while the latter can be viewed as due to a reduction of the time step for phase diffusion. With second-order decoupling, K = 2, the decoherence rate is additionally suppressed, and with time-reversal invariant bath coupling all orders of the expansion in powers of adiabaticity parameter may vanish, in which case the leading-order dephasing term becomes exponentially small and dephasing would likely be determined by terms of higher order in bath coupling. Along with the decoherence rates characterizing the exponential decay of quantum correlations with time, the corresponding prefactor, which determines the "visibility" (or "initial decoherence" 33 ), was also analyzed 30 . While for generic refocusing sequences with K ≥ 1 the initial decoherence is quadratic in τ c and does not scale with the thermal bath correlation time τ 0 , for symmetric pulse sequences it is reduced by an additional power of the adiabaticity parameter (τ /τ 0 ). These results were originally derived for a generic featureless bath, but they also hold in a vicinity of a sharp resonance as long as the effective (i.e., renormalized as in the average Hamiltonian) coupling to the corresponding mode is small compared to its width 20 .

V. APPLICATION EXAMPLES.
A. Decoupling sequences for a chain of qubits Decoupling sequences are designed to prevent quantum evolution from happening. Thus, we want to construct a sequence such that the resulting evolution operator over the period τ c is identity, U (τ c ) = 1 1. We illustrate the scalability of dynamical decoupling to large system sizes by considering linear chains of qubits with either Ising or XXZ n.n. random-valued couplings [only J z n,n+1 or both J z n,n+1 and J x n,n+1 = J y n,n+1 in Eq. (7)], plus the local fields either along z axis or in arbitrary direction [∆ z n = 0 or ∆ µ n = 0 for µ = x, y, z in Eq. (7)]. With such a system Hamiltonian, zeroth-order average Hamiltonian (16) contains only individual qubits or pairs of neighboring qubits, the largest clusters contributing to the 1st-order average Hamiltonian (17) originate from two bonds sharing a site (three qubits), and in general H (n) contains terms spanning contiguous clusters of up to n + 1 bonds, that is, n + 2 qubits. Thus, to design a K-th order decoupling sequence, one needs to consider individual clusters of up to K + 1 qubits.
With nearest-neighbor and local couplings only, the decoupling can be implemented by simultaneously applying pulses on either odd or even sublattice. We note that in our setup there is no gap between subsequent pulses, the pulses follow back to back with the repetition period τ . The system is "focused" at the end of each cycle consisting of several pulses of length τ . Such a scheme with a common "clock" time τ is convenient, e.g., for parallel execution of quantum gates in different parts of the system. For each qubit, various pulses (or intervals of no signal) can be executed in sequence.
In this work we consider the following two sequences from Ref. 18, 4 = X 1 Y 2X1Ȳ2 and its symmetrized version 8 = X 1 Y 2X1Ȳ2Ȳ2X1 Y 2 X 1 , which provide universal refocusing of the couplings between the sublattices, and also suppress the on-site chemical shifts ∆ z n . Here, X 1 is a π x pulse simultaneously applied on all odd sites,Ȳ 2 is a (−π) y = π −y pulse applied on all even sites, etc. These sequences are "best" sequences at given length for all pulse shapes found by exhaustive search (high-order sequences 34,35 equivalent for hard pulses do not necessarily have equal orders here). The fact that such a bruteforce optimization approach works is entirely due to the efficiency of the numerical method.
In addition, we constructed two longer sequences, , and its symmetrized version 32, constructed by running the sequence 16 first directly and then in reverse order. Here 0 denotes zero pulse, an empty interval of duration τ . These two sequences provide universal decoupling both for any couplings between the sublattices and for arbitrary on-site fields (∆ µ n = 0). In addition to the system Hamiltonian, the effectiveness of a sequence application depends on the quality of the pulses. In Table III, we list orders of the sequences when applied with different pulse shapes, computed using the numerical time-dependent perturbation theory as described in sec. II D. The term R k (τ c ) was considered to be zero if its norm vanished with numerical precision, typically 10 −8 or better, compared to typical values of order one for orders where R k (τ c ) = 0. The orders K do not depend on the chain length; we verified this statement on chains up to n = 7 qubits. Also, the computed orders are the same for all self-refocusing pulse shapes of particular order; we believe that the results will remain valid for other symmetric pulse shapes of the same order as indicated in the 1st column of Table III.

B. Error scaling
We illustrate the predicted power laws in Fig. 6, where the average infidelity (A4) is computed for different ratios of t/τ , where t is the fixed evolution time and the pulse duration τ was reduced to accommodate a different number of decoupling cycles. The simulation is done for chains of n = 4 qubits with randomly chosen but fixed parameters corresponding to different chain models as indicated. The steepest lines correspond to largest order K of the sequence decoupling order. For symmetric sequence 8 with Ising chain, K = 2 for Gaussian pulses, Fig. 6(a), K = 4 for 1st-order pulses, Fig. 6(b), and K = 6 for 2nd-order pulses, Fig. 6(c). The corresponding infidelities for fixed evolution time scale as ∝ (J z τ ) 4 , ∝ (J z τ ) 8 , and ∝ (J z τ ) 12 . Larger values of K can improve accuracy by orders of magnitude, or, at fixed required fidelity, substantially reduce the number of decoupling cycles.   (14) is H (K) , so that for small enough τ the mismatch in the unitary evolution operator after n decoupling cycles (evolution time t = nτc) scales as U − 1 1 ∝ tτ K c , and the corresponding infidelity 1 − F ∝ t 2 τ 2K cycle . Sequence 8 a sequence of 8 pulses applied intermittently on odd or even sublattices, see text for actual definitions.
We saw that with order-K decoupling in multi-qubit systems with local couplings, the decoupling error operators can be represented as connected clusters of up to K + 1 bonds. For a linear chain, these involve up to K + 2 qubits, and the number of such operators scales linearly with the total number n of qubits, as long as n > K + 2. In an n-qubit system, each of such operators can be written as an outer product of the cluster contribution, and the identity operators for the remaining qubits. As a result, the square of the Frobenius norm of the error operators scales linearly with the size of the Hilbert space, that is, exponentially with the number of qubits. However, this exponential scaling is suppressed when we compute the infidelity [see Eq. (A4)], so that the infidelity scales only linearly with the number of clusters, that is, linearly with the number of qubits. The same scaling with the system size is expected in higher dimensional arrangements of qubits (planar, 3D).
We illustrate the scaling of decoupling errors with the qubit number n in Fig. 7. The plots show the scaling of the average infidelity at the end of the interval in Fig. 6 and other data with the chain length n.

C. Composite pulses
Composite pulses are, in fact, pulse sequences designed to replace a single pulse and specially designed to compensate for some particular systematic errors, including off-resonance application, pulse amplitude, and pulse phase errors 36,37,38,39,40,41,42,43,44 .   The off-resonance errors appear when the carrier frequency of the applied pulse is off the transition frequency between the |0 and |1 state of a qubit. In the rotating reference frame this is equivalent to a non-zero chemical shift Hamiltonian (3), with ∆ equal to the corresponding frequency bias. We note that our 1st and 2nd-order selfrefocusing pulse shapes already offer a degree of stability against such errors.
For this reason we concentrate on the pulse amplitude errors, where the correct pulse shape is applied with the wrong amplitude, producing an incorrect rotation anglẽ φ 0 = φ 0 . Note that no one-dimensional pulse shaping can compensate for this kind of errors, since the modified rotation angle is simply proportional to the pulse amplitude, On the other hand, one can expect that the pulse amplitude offset f remains the same for all the pulses applied at a particular frequency. This uniformity is utilized in several composite pulses designed so that the net rotation would be insensitive to such uniform errors.

SCROFULOUS
The three-pulse sequence SCROFULOUS 38 is based on the sequence originally proposed by Tycko 45,46 . Particularly, an improved π pulse is obtained by applying three π pulses, at 60 • , 300 • , and again at 60 • , or just π 60 π 300 π 60 . In the case of ideal δ-pulses, the resulting pulse compensates for pulse amplitude errors to linear order. With finite-width shaped pulses, an additional error is generated due to the presence of the system Hamiltonian. In particular, for the chemical-shift Hamiltonian (3) Combining the corresponding expressions appropriately rotated in the x-y plane, and expanding the result to quadratic power in the relative amplitude offset f , we obtain for the composite pulse π 60 π 300 π 60 , whereυ is the parameter υ [Eq. (40)] but for the pulse with rescaled amplitude, and the further terms are of order fυτ ∆,υ 2 τ 2 ∆ 2 ,αf τ 2 ∆ 2 . Clearly, with Gaussian or other pulse shape such that υ = 0, the error is linear in τ ∆ and quadratic in the amplitude shift f (although generally there will also be a cross-term ∝ f τ ∆). This situation is illustrated in Fig. 8(a), where the average infidelity is plotted for the SCROFULOUS sequence with pulses G 010 on the plane ∆τ -f . The region for 1 − F = 10 −5 is a narrow vertical line which corresponds to great sensitivity to frequency shift. With 1st-order self-refocusing pulses such that υ = 0,υ(f ) ∝ f , the error is dominated by the term ∝ f τ ∆. The corresponding region corresponds to a diamond-like shape in the center of Fig. 8(b). We have also generated self-refocusing pulse shapes such that both υ = 0 and υ ′ = 0. Then, by symmetry, υ ′′ = 0, and genericallỹ υ ∝
0.136, compared with |f | < 0.090 for the SCROFULOUS sequence. Thus, even though order of the BB 1 family of composite pulses is higher (and thus, for small |f | their performance is much better asymptotically), at this level of infidelity their performance is comparable. We now turn to off-resonance correction terms which differ between implementations of the BB 1 sequence. In particular, with generic pulses such that υ = 0, already at f = 0, all of these sequences acquire linear corrections scaling with τ ∆. For example, the expansion of the sequence BB (W ) 1 at f = 0 can be written as where α ≡ 2α 1 −α 2 , υ 2 ≡ υ 2 1 +υ 2 2 , and the parameters υ 1 , α 1 and υ 2 , α 2 correspond to the π and 2π pulses respectively. Notice that the second-order coefficients α i enter only in the combination 2α 1 − α 2 . Not surprisingly, if we replace the 2π-pulse with two π pulses, the coefficients α cancel out, With specially designed pulses such that both υ = 0 and υ ′ = 0, due to pulse symmetry, also υ ′′ = 0, so that υ(f ) = υ ′′′ f 3 + O(f 4 ). Then, for sequences other than BB , with 1st-order pulses the error is dominated by the terms quadratic in τ ∆ due to coefficients α i [cf. Eq. (53)]. As a result, the high-fidelity regions in Figs. 9(c), 10(b), and 11(a) have much more rounded shape. With 2nd-order pulses with amplitude correction, such that α i = 0 but α ′ i = 0, the leading-order error scales as O(f τ 2 ∆ 2 ), which extends the high-fidelity regions out to larger values of τ ∆ in a characteristic "smile" pattern. For the sequence BB (W ′ ) 1 these terms cancel out, and the leading-order error term comes from the non-zero υ ′′′ , the errors scale as ∝ f 3 τ ∆, This error has the same order in f as that in Eq. (52), and it can compensate or increase the contribution linear in σ z . The result is a somewhat skewed in the center high-fidelity region widely stretched horizontally.

D. Stability of decoupling against amplitude errors
We now return to the problem of decoupling for a chain of qubits, but now consider the effect of the amplitude errors. This needs a separate study since single-qubit errors with composite pulses have structure different from those due to, say, finite pulse width.
In Fig. 12(a) we present the results of simulations for a particular 4-qubit Ising chain with on-site chemical potential shifts ∆ z i over time interval t = 128τ (exactly the same parameters as in Figs. 6, 7). Specifically, we plot the infidelity (1 − F ) in units of 10 −4 , as a function of the fractional pulse amplitude 1 + f . In the first half of the symmetric sequence 8 = X 1 Y 2X1Ȳ2Ȳ2X1 Y 2 X 1 , we used the original BB (W ) 1 pulse X → π 0 π φ (2π) 3φ π φ for π x , rotated appropriately to implement Y → π π/2 π φ+π/2 (2π) 3φ+π/2 π φ+π/2 as well as X ≡ X π , while in the second part of the sequence we used the same decomposition but backwards, e.g., With no amplitude mismatch, thus constructed decoupling sequence has order K = 1 with 1st-order pulses, and an almost-2nd order "1*" with 2nd-order pulses, with the norm of the 2nd-order term reduced by three orders of magnitude compared to 1st-order pulses. These should be compared with K = 1 and K = 3 respectively for the regular 8-pulse sequence [Tab. III].
It is seen from Fig. 12(a) that regular 1st-and 2ndorder pulses S 1 and Q 1 perform in the BB 1 -corrected composite sequence almost as well as the self-refocusing pulses with additional amplitude protection (1st-order L 1 and 2nd-order L 3 ). In addition, the 2nd-order pulses with amplitude protection (L 3 ) work almost as well in the regular 8-pulse sequence. All of these allow to achieve the infidelity level of 10 −4 at amplitude mismatch of 3% or higher (up to 6% with pulses L 3 in the sequence 8BB 1 ). With the regular 8-pulse sequence, the region of high fidelity shrinks substantially for pulses Q 1 and L 1 , and it all but disappears for the pulse S 1 [the coefficient α for the pulse L 1 happens to be a few times smaller than that for S 1 ].
five orders of magnitude with pulses Q 1 [Fig. 6(c)], and by some two orders of magnitude for pulses S 1 [ Fig. 6(b)]. As for the Ising chain, the effect of the amplitude errors is weaker with specially-designed pulses L 1 and L 3 such that the 1st-order coefficientυ(f ) scales as a higher power of f . With the pulse L 1 , the 2nd-order coefficient α is non-zero but small; it is seen from Fig. 12(b) that its effect is to introduce a linear term 1 − F ∝ f which tends to skew the infidelity minimum away from f = 0. We should also note that with the Gaussian pulses G 010 (not shown), even the on-resonance infidelity is out of range of the plots Fig. 12.  FIG. 11: (color online) As in Fig. 9 but for the sequence BB W ′ 1 [π0π φ π 3φ π 3φ π φ ] using only π-pulses. (a) 1st-order pulses with amplitude correction; (b) 2nd-order pulses with amplitude correction. The absence of the error terms linear in α [see Eqs. (54), (55)] produces a much wider high-fidelity region already with 1st-order pulses.

VI. CONCLUSIONS
We presented a comprehensive study targeting pulse and sequence design and analysis based on a consistent high-order average Hamiltonian expansion. The numerical technique for expanding the evolution operator was originally introduced by us in Ref. 18, and a complimentary analytical technique was developed for π-pulses by one of the authors in Ref. 20.
The overall approach is to start with a closed system described by a finite-dimensional Hamiltonian H S and design a sequence of shaped pulses such that the evolution operator would be accurate to a given order K in powers of H S . The key to this approach are the NMR-style 1st-and 2nd-order self-refocusing one-dimensional pulses constructed for a single-qubit chemical-shift Hamiltonian (3). In this work we designed a number of such shapes for different rotation angles φ 0 , and presented a careful analytical analysis of the first two leading orders of the average Hamiltonian theory for driven qubit evolution with the most general system Hamiltonian H S . While any symmetric one-dimensional pulse shape is characterized by only three parameters, two of these can be set to zero by pulse shaping. The remaining parameter is also non-zero for an ideal "hard" δ-function pulse. This leads to an important conclusion that the constructed pulses can be used as drop-in re-   (a) Ising chain in the presence of individual chemical shifts ∆ z i . 8-pulse sequence with BB1 composite pulses offers the best accuracy, which remains essentially the same whether the sequence is used with regular 1st-or 2nd-order pulses or with the pulses stabilized against amplitude errors (1st-order pulses L1 and 2nd-order pulses L3). However, the pulse L3 works well enough even with regular 8-pulse sequence. While the details of the amplitude scaling differ, at the level of 1 − F = 10 −4 , the 1st-order amplitude-protected pulses L1 and regular 2nd-order pulses Q1 have comparable accuracy. The use of 1st-order pulses show relatively poor performance even on resonance. (b) XXZ chain in the presence of individual chemical shifts ∆ z i . With XXZ coupling, The BB1 composite pulse is no longer accurate, as the errors appear already in the linear order in τ J ⊥ (not shown). With the regular 8-pulse sequence 8, the best accuracy is obtained for the pulses with amplitude correction. placement for hard pulses; with proper pulse placement the results should be identical to first two orders. The structure of errors appearing in higher orders of the evolution operator can be understood by analyzing the numerical time-dependent perturbation series for the evolution operator of a closed system.
An important advantage of this approach is that the expansion order offers a natural classification of the error operators. As a result, (i) the convergence regions have regular shapes as a function of parameters [see Figs. 10(c,d) and 11(b,c)]. Furthermore, with local two-(or few-)qubit couplings dominant, (ii) the error operators can be placed on connected clusters of up to k + 1 qubits for terms of order k, which allows one to understand their structure in terms, e.g., the direct products of up to (k + 1) Pauli matrices. Once their structure understood, the convergence can be readily improved by suppressing the error operators, as in our analysis of pulseamplitude errors. We emphasize, that such an analysis can be performed even for very large qubit systems. Thus, (iii) this approach is characterized by scalability with the system size, as we illustrated by analyzing decoupling infidelity with the system size [ Fig. 7]. Although in this work we concentrated on the dynamics of closed systems, another important advantage is that (iv) the high-order control sequences result in lower decoherence in the presence of slow environmental modes 20,30 .
Most obvious application of highly-optimized shaped pulses of the sort presented in this work is in solid-state quantum computation, where the bandwidth available for quantum gates is typically limited. Our techniques based on analytical and numerical high-order average Hamiltonian theory offers a systematic scalable approach for constructing gates for such multi-qubit systems, without need of solving their full dynamics. However, even if the bandwidth does not appear to be at premium, simple pulse shaping (e.g., using 1st-order pulses) can still offer a substantial improvement of control accuracy.

VII. ACKNOWLEDGMENTS.
This research was supported in part by the NSF grant No. 0622242 (LP). LANL is supported by US DOE under Contract No. W-7405-ENG-36. NHMFL is supported by the DOE, the NSF, and the state of Florida.

APPENDIX A: AVERAGE FIDELITY
Here we discuss the calculation of the fidelity averaged over the initial state, in the case of unitary evolution with known evolution matrix U , while the desired evolution matrix is U 0 . Let us write the density matrix of the initial pure state as ρ 0 = ψψ † , where ψ is an N -component complex vector. Then the actual density matrix is ρ = U ψψ † U † , while the desired density matrix is ρ ideal = U 0 ψψ † U † 0 . The fidelity with the given initial state The only condition on the components ψ i of the wavefunction is the normalization, 1 = ψ † ψ = i |ψ i | 2 . Generally, this means that the average of the product in Eq. (A1) can only depend on the identity tensor δ ij . By symmetry, ψ i ψ * j ψ k ψ * l = A(δ ij δ kl + δ il δ kj ), where the unknown coefficient A can be computed from the normalization by tracing over i = j, k = l. We obtain 1 = A(N 2 + N ), so that the average fidelity where V ≡ U † 0 U . Numerically, with V close to identity matrix, the loss of precision can be avoided by expressing the infidelity 1 − F in terms of the modified mismatch, Namely, since δ 2 = 2N − 2| tr V |, the average infidelity can be written as