Filter Function Formalism and Software Package to Compute Quantum Processes of Gate Sequences for Classical Non-Markovian Noise

Correlated, non-Markovian noise is present in many solid-state systems employed as hosts for quantum information technologies, significantly complicating the realistic theoretical description of these systems. In this regime, the effects of noise on sequences of quantum gates cannot be described by concatenating isolated quantum operations if the environmental correlation times are on the scale of the typical gate durations. The filter function formalism has been successful in characterizing the decay of coherence under the influence of such classical, non-Markovian environments and here we show it can be applied to describe unital evolution within the quantum operations formalism. We find exact results for the quantum process and a simple composition rule for a sequence of operations. This enables the detailed study of effects of noise correlations on algorithms and periodically driven systems. Moreover, we point out the method's suitability for numerical applications and present the open-source Python software package filter_functions. Amongst other things, it facilitates computing the noise-averaged transfer matrix representation of a unital quantum operation in the presence of universal classical noise for arbitrary control sequences. We apply the presented methods to selected examples.


I. INTRODUCTION
In the circuit model of quantum computing, computations are driven by applying time-local quantum gates. Any algorithm can be compiled using sequences of oneand two-qubit gates [1]. Ideal, error-free gates are represented by unitary transformations, so that simulating the action of an algorithm on an initial state of a quantum computer amounts to simple matrix multiplication. Real implementations are subject to noise that causes decoherence resulting in gate errors. If the noise is uncorrelated between gates, its effect can be described by quantum operations acting as linear maps on density matrices, even when several gates are concatenated. A closely related approach is the use of a Master equation in Lindblad form [2], which governs the dynamics of density matrices under the influence of Markovian noise with a flat power spectral density.
Yet many physical systems used as hosts for qubits do not satisfy the condition of uncorrelated noise. One example frequently encountered in solid state systems is that of 1/f noise, which in principle contains arbitrarily long correlation times. It emerges for instance as flux noise in superconducting qubits and electrical noise in quantum dot qubits [3][4][5][6]. Whereas simple approaches exist to treat for example quasistatic noise, which corresponds to perfectly correlated noise (i.e. a spectrum with weight only at zero frequency), they cannot be applied to 1/f noise because of the wide distribution of correlation times it contains. Thus, there is a gap in the mathematical descriptions of gate operations for noises with arbitrary power spectra that exist between the extremal cases of perfectly flat (white) and sharply peaked * tobias.hangleiter@rwth-aachen.de † pascal.cerfontaine@rwth-aachen.de (quasistatic) spectra. To capture experimentally relevant effects important to understand the capabilities of quantum computing systems, a universally applicable formalism is hence desirable. For example, one may expect the fidelity requirements for quantum error correction to be more stringent for correlated noise as errors of different gates can interfere constructively [7]. On the other hand, it might also be possible to use correlation effects to one's benefit, attenuating decoherence by cleverly constructing the gate sequences in algorithms. As experimental platforms begin to approach fidelity limits set by employing primitive pulse schemes [5,8,9] and detailed knowledge about noise sources and spectra in solid-state systems becomes available [10][11][12], control pulse optimizations tailored towards specific systems will be required to further push fidelities beyond the error correction threshold [13,14]. This calls for flexible and generically applicable tools as a basis for the numerical optimization of pulses as well as the detailed analysis of the quantum processes they effect. In order to obtain a useful description also for gate operations that decouple from leading orders of noise, such as dynamically corrected gates (DCGs) [15], beyond leading order or exact results are required.
In an accompanying publication [16] we presented a formalism based on filter functions and the Magnus expansion (ME) that addresses these needs and limitations of the canonical master equation approach for correlated noise. Specifically, we showed how process descriptions can be obtained perturbatively for arbitrary classical noise spectra and derived a concatenation rule to obtain the filter function of a sequence of gates from those of the individual gates. This paper extends these results.
Filter functions (FFs) were originally introduced to describe the decay of phase coherence under dynamical decoupling (DD) sequences [17][18][19][20] consisting of wait times and perfect π-pulses. The formalism facilitated recog-nizing these sequences as band-pass filters that allow for probing the environmental noise characteristics of a quantum system through noise spectroscopy [12,[21][22][23] or optimizing sequences to suppress specific noise bands [24][25][26][27]. It can also be extended to fidelities of gate operations for single [28,29] or multiple [30,31] qubits using the ME [32,33] as well as more general DD protocols [34]. The works by Green et al. [29] and Clausen et al. [35] also introduced the notion of the control matrix as a quantity closely related to the canonical filter function that is convenient for calculations. In this context, the formalism's capability to predict fidelities of gate implementations has been identified and experimentally tested [26,29,36,37]. Recently, it has also proved useful in assessing the performance requirements for classical control electronics [38].
While analytic approaches allow for the calculation of filter functions of arbitrary quantum control protocols in principle, it is in practice often a tedious task to determine analytic solutions to the integrals involved if the complexity of the applied wave forms goes beyond simple square pulses or extend to multiple qubits. Moreover, one does not always have a closed-form expression of the control at hand, such as is the case for numerically optimized control pulses. This calls for a numerical approach which, while giving up some of the insights an analytical form offers, is universally applicable and eliminates the need for laborious analytic calculations.
Here, we build and extend upon our accompanying work of Ref. 16 and that of Ref. 29 to show that the formalism can be recast within the framework of stochastic Liouville equations by means of the cumulant expansion [39,40] which, for Gaussian noise, entails exact results for the quantum process of an arbitrary control operation using only first and second order terms of the ME [32]. Moreover, due to the fact that the ME retains the algebraic structure of the expanded quantity [33] we are able to separate decoherent and coherent contributions to the process. We give explicit methods to evaluate these terms for piecewise-constant control pulses. Moreover, we show that the formalism naturally lends itself as a tool for numerical calculations and present the filter_functions Python software package that enables calculating the filter function of arbitrary, piecewise constant defined pulses [41]. On top of providing methods to handle individual quantum gates, the package also implements the concatenation operation as well as parallelized execution of pulses on different groups of qubits, allowing for a highly modular and hence computationally powerful treatment of quantum algorithms in the presence of correlated noise. Given an arbitrary, classical noise spectral density, it can be used to calculate a matrix representation of the error process. From this matrix one can extract average gate fidelities, transition probabilities, and leakage rates as we derive below. To simplify adaptation the software's API is strongly inspired by and compatible with QuTiP [42] as well as qopt [43]. This allows users to use these packages in conjunction. Assess-ing the computational performance, we show that our method outperforms Monte Carlo simulations for single gates. New analytic results applicable to periodic Hamiltonians and employing the concatenation property make this advantage even more pronounced for sequences of gates. To highlight the main software features, we show example applications below.
We provide this package in the expectation that it will be a useful tool for the community. Besides recasting and expanding on our earlier introduction of the formalism in Ref. 16, the present paper is intended to provide an overview of the software and its capabilities. It is structured as follows: In Section II we derive an exact expression for unital quantum operations in the presence of non-Markovian Gaussian noise and lay out how it may be evaluated using the filter function formalism. We review the concatenation of quantum operations shown in Ref. 16 and furthermore adapt the method by Green et al. [29] to calculate the filter function of an arbitrary control sequence numerically. We will specifically focus on computational aspects of the formalism and lay out how to compute various quantities of interest. Moreover, we classify its computational complexity for calculating average gate fidelities and remark on simplifications that allow for drastic improvements in performance in certain applications. In Section IV, we introduce the software package by outlining the programmatic structure and giving a brief overview over the API. Lastly, in Section V, we show the application of the software by means of four examples that highlight various features of the formalism and its implementation. Therein, we first demonstrate that the formalism can predict average gate fidelities for complex two-qubit quantum gates in agreement with computationally much more costly Monte Carlo calculations. Next, we show how it can be applied to periodically driven systems to efficiently analyze Rabi oscillations. We finally establish the formalism's ability to predict deviations from the simple concatenation of unitary gates for sequences and algorithms in the presence of correlated noise by simulating a randomized benchmarking experiment as well as assembling a quantum Fourier transform (QFT) circuit from numerically optimized gates. We conclude by briefly remarking on possible future application and extension of our method in Section VI.
Throughout the paper we will denote operators by Roman font, e.g. U , and quantum operations and their representations as transfer matrices by calligraphic font, e.g. U, which we also use for the control matrixB to emphasize its innate connection to a transfer matrix. For consistency, a unitary quantum operation will share the same character as the corresponding unitary operator. An operator in the interaction picture will furthermore be designated by an overset tilde, e.g.H = U † HU with U the unitary operator defining the co-moving frame. Definitions of new quantities on the left and right side of an equality are denoted by := and =:, respectively. We use a central dot ( · ) as a placeholder in some definitions of abstract operators such as the Liouvillian, denoted by L := [H, · ], which is to be understood as the commutator of the corresponding Hamiltonian H and the operator that L acts on. The identity matrix is denoted by 1 and its dimension always inferred from context. Furthermore, we will use Greek letters for indices that correspond to noise operators in order to distinguish them clearly from those that correspond to basis or matrix elements. Lastly, we work in units where = 1.

II. FILTER FUNCTION FORMALISM FOR UNITAL QUANTUM OPERATIONS
We begin the theoretical part of this article by showing how a superoperator matrix representation of the error process, the "error transfer matrix", of a unital quantum operation can be computed from the control matrix of the pulse implementing the operation. The control matrix relates the operators through which noise couples into the system to a set of basis operators in the interaction picture and we detail how it can be calculated in a relatively efficient manner for two different situations. First, we consider a sequence of gates whose control matrices have been precomputed. Second, we lay out how the control matrix can be obtained from scratch under the assumption of piecewise constant control, which is often convenient for approximating continuous pulse shapes. Other wave forms can be dealt with analogously by solving the corresponding integrals. We then move on to show how several quantities of interest can be extracted and present optimized strategies for computing the central objects of the formalism.
A. Transfer matrix representation of quantum operations

Brief review of quantum operations and superoperators
The quantum operations formalism provides a general framework for the description of open quantum systems [44,45]. It forms the mathematical basis for quantum process tomography (QPT) [46,47] as well as gate set tomography (GST) [48,49] and has also been extensively employed in the context of randomized benchmarking (RB) [50,51]. Several different representations of quantum operations exist. While all of them are equivalent one typically chooses the most convenient for the problem at hand. For an overview of the most commonly used representations see Ref. 49 and for matrix representations in particular Ref. 52 and the references therein. In this work we employ the Liouville representation, to the best of our knowledge first formalized by Fano [53], to profit from its simple properties under composition. It is also known as the transfer matrix representation and we will use the terms interchangeably below. We now briefly review the concept and refer the reader to the literature for further details. Concretely, the Liouville representation of an operation E : ρ → E(ρ) acting on density operators in a Hilbert space H of dimension d is given by with an operator basis C = {C 0 , C 1 , . . . , C d 2 −1 } for the space of linear operators over H , L(H ), orthonormal with respect to the Hilbert-Schmidt product A, B := tr A † B . In the case that the operator basis corresponds to the Pauli matrices Eq. (1) is known as the Pauli transfer matrix (PTM). The operation E is thus associated with a d 2 × d 2 matrix in Liouville space L that describes its action as the degree to which the j-th basis element is mapped onto the i-th. On L one can identify a set of basis kets isomorphic to the operators C i (and correspondingly bras i| to the adjoint C † i ) as well as the inner product i|j = C i , C j . As the vectors |i form an orthonormal basis, any operator on H can be written as a vector on L , |A = i |i i|A , whereas a superoperator on H becomes a matrix on L , see Eq. (1). It can then be shown that density operators represented by vectors are propagated by transfer matrices so that the action of a quantum operation E on a density operator ρ is given by |E(ρ) = E |ρ = ij |i i|E|j j|ρ . Thus, the composition of two operations E 1 and E 2 corresponds to matrix multiplication in Liouville space, a property which makes the representation particularly attractive for sequences of operations. Although from a numerical perspective the computational complexity scales unfavorably with the system dimension d (c.f. Section III D), we will employ the Liouville representation for its transparent interpretation and concise behavior under composition in the following analytical considerations. Lastly, we note that for C 0 ∝ 1, trace-preservation and unitality are encoded in the relations E 0j = δ 0j and E j0 = δ j0 , respectively.

Liouville representation of the error channel
We will now derive an expression for the quantum process of a quantum gate in the presence of arbitrary classical noise. As a single realization of a classical noise generates strictly unitary dynamics, we will be interested in the expectation value of the dynamics over many such realizations, which will lead to a quantum process including decoherence. If the noise is additionally Gaussian, these results are exact and therefore apply without restrictions to arbitrarily large noise strength as well as to gates that partially decouple from noise. For such DCGs or DD sequences [15,20] higher order terms can become dominant. In the case that the environment is not strictly Gaussian, our approach becomes perturbative and we recover the results presented in Ref. 16. As most of our discussion later on in this article will focus on the leading order approximation, readers not interested in the full generality may refer to that publication for a less general but perhaps more accessible derivation and skip ahead to Section II B.
The difference is that in Ref. 16, the Magnus expansion is applied to the solution of the Schrödinger equation, whereas the approach presented here is based on the theory of stochastic Liouville equations and the cumulant expansion [39,40]. In the filter function context, the cumulant expansion has been used to express the decay of the off-diagonal terms of a single-qubit density matrix in Ref. 20. More recently, Paz-Silva and Viola [34] employed it in conjunction with the ME to obtain the matrix elements of the perturbed density operator after a time T of noisy evolution. In Ref. 54, the authors made use of the cumulant expansion and stochastic Liouville equations for the purpose of gate optimization. Here, we combine different aspects of these works and make the connection to the quantum operations formalism by determining the noise-averaged error propagator in the Liouville representation. This form completely characterizes the error process and hence allows for detailed insight into the decoherence mechanisms of the operation.
Concretely, we consider a system described by the stochastic Hamiltonian H c (t) is implemented by the experimentalist to generate the desired control operation during the time t ∈ [0, τ ] and H n (t) describes classical fluctuating noise environments b α (t) ∈ R that couple to the quantum system via the Hermitian noise operators B α (t) ∈ L(H ). These may carry a general, deterministic time dependence and without loss of generality, we can require them to be traceless since any contributions proportional to the identity do not contribute to noisy evolution in any case [55]. The b α (t) are random variables drawn from (not necessarily Gaussian) distributions with zero mean that are assumed to be independent and identically distributed (i.i.d.) both with respect to repetitions of the experiment. Note that this concept of independence does not preclude correlations between different noise sources α = β nor between one noise source at different times t = t , but only serves to obtain a well-defined ensemble average. Lastly, to be able to later on relate the correlation functions of the b α (t) to their spectral density, we require the noise fields to be wide-sense stationary, meaning that their correlation function depends only on the time difference. For noise operators without explicit time dependence, Eq. (3) constitutes a universal decomposition as can be seen by choosing the B α from an orthonormal basis for L(H ). To motivate the time-dependent form of Eq. (3), assume the true Hamiltonian is a function of a set of noisy parameters λ (t) = λ(t) + δλ(t) where δλ(t) = vec({b α (t)} α ) are the stochastic variables. Ex-panding the Hamiltonian in an orthonormal operator basis yields H( λ (t)) = α f α ( λ(t), δλ(t))B α . In general, however, the expansion coefficients f α will be arbitrary functions of both the deterministic parameters λ(t) and the stochastic noises δλ(t), which prohibits a factorized form like Eq. (3). We can address this problem by first expanding H around λ(t) for small fluctuations δλ(t). Then, the Hamiltonian approximately becomes H( λ (t)) ≈ H( λ(t)) + δλ(t) · ∇ λ H( λ(t)), where we can define the control Hamiltonian as H c (t) := H( λ(t)). Expanding the second term in the operator basis now results in the form (3) for the noise Hamiltonian as it is linear in δλ(t) and the deterministic time dependence is contained in ∇ λ H( λ(t)) alone.
This permits us to model complex relations between physical noise sources and the noise operators that capture the coupling to the quantum system, arising for example through control hardware or effective Hamiltonians obtained from e.g. Schrieffer-Wolff transformations. While the linearization is in most cases an approximation, it does not impose significant constraints since the noise is typically weak compared to the control [56]. As an example, we could capture a dependence of the device sensitivity on external controls (see also Ref. 30). In a widely used setting electrons confined in solid-state quantum dots are manipulated using the exchange interaction J that depends non-linearly on the potential difference between two dots. Since the dominant physical noise source affecting this control is charge noise, one could include the effect on J( ) to first order with s (t) = ∂J( (t))/∂ (t) so that H n (t) = b (t)B (t) = b (t)s (t)B for some operator B which represents the exchange coupling.
We proceed in our derivation by noting that the control Hamiltonian H c gives rise to the noise-free Liouville-von Neumann equation where the action of U c on a state ρ is to be understood as U c : ρ → U c ρU † c with U c the usual time evolution operator satisfying the corresponding Schrödinger equation. This allows us to write the superpropagator for the total Liouvillian L = L c + L n as U(t) = U c (t)Ũ(t) where the unitary error superpropagatorŨ(t) contains the effect of a specific noise realization in Eq. (3). Next, we transform the noise Liouvillian L n to the interaction picture with respect to the control Liouvillian L c so thatŨ(t) satisfies the modified Liouville equation Equation (5) may be formally solved using the Magnus expansion [32] so that at time t = τ with L eff (τ ) = ∞ n=1 L eff,n (τ ). A sufficient criterion for the convergence of the expansion is given by Moan et al. [57] as τ 0 dt L n (t) < π where · = · , · is the Frobenius (Hilbert-Schmidt) norm. The first and second terms of the ME are given by [32,33] The n-th term of the expansion contains n factors of the noise variables b α (t) and scales with n factors of the control duration τ , suggesting that higher-order terms can be neglected if their product is small. In the Bloch sphere picture this corresponds to requiring that the angle by which the Bloch vector is rotated away from its intended trajectory due to the noise be small. Below, we will use the parameter ξ to denote the magnitude of this deviation. It is properly defined in Appendix D 1 where also bounds for the convergence of the ME are discussed.
Here, we only state that L eff,n ∼ ξ n (see also Ref. 29).
We have suggestively written the ME in terms of an effective Liouvillian L eff = [H eff , · ] to interpret it as the generator of a time-averaged evolution of a single noise realization up to time τ . In order to obtain the ensembleaveraged evolution of many realizations of the stochastic Hamiltonian in Eq. (3), we apply the cumulant expansion toŨ (see also Refs. 58 and 59), with · denoting the ensemble average [60] and the cumulant function [39] The notation · c denotes the cumulant average which prescribes a certain averaging operation. The first cumulant of a set of random variables {X i (t)} i is simply the expectation value, X i (t) c = X i (t) , whereas the second cumulant corresponds to the covariance, Remarkably, third and higher-order cumulants vanish for Gaussian processes [40,61], making Eq. (11) exact by truncating the sums already at k = 2 and n = 2. In this case, the convergence radius of the ME becomes infinite. The terms with k = n = 2 do not contribute as they involve fourth-order cumulants. Since furthermore we assume that the noise fields b α (t) have zero mean, also the terms with k = n = 1 vanish and X i (t)X j (t) c = X i (t)X j (t) . We can hence write the cumulant function succinctly as Equations (9) and (12) allow us to exactly compute the full quantum process Ũ : ρ → Ũ (ρ) for Gaussian noise with arbitrary spectral density and power. For non-Gaussian noise these expressions are approximate up to O ξ 2 and higher order terms include both higher orders of the ME and the cumulant expansion. Inspecting Eq. (12), we observe that the first term is anti-Hermitian as it is a pure Magnus term (remember that the ME preserves algebraic structure to every order) and thus generates unitary, coherent time evolution. Conversely, the second term is Hermitian and thus generates decoherence [62]. The former is more difficult to compute than the latter because the second order of the ME, Eq. (8b), contains nested time integrals. Arguments can be made [16], however, that for single gates in an experimental context the coherent errors captured by this term can be calibrated out to a large degree [63,64]. Moreover, many of the central quantities of interest that can be extracted from the quantum process, among which are gate fidelities and certain measurement probabilities, are functions of only the diagonal elements of K. By virtue of the antisymmetry of the second order terms, they do not contribute to these quantities to leading order as we show in Section II D. While we will also lay out how to compute the second order, our discussion will therefore focus on contributions from the incoherent term below. As it turns out, this term can be computed using a filter function formalism based on that by Green et al. [29]. To see this, we insert the explicit forms of the ME given in Eq.
is the cross-correlation function of noise sources α and β which we will later relate to the spectral density. For now, we stay in the time domain and introduce an orthonormal and Hermitian operator basis for the Hilbert space H to define the Liouville representation, where we choose C 0 = 1 d 1/2 for convenience so that the remaining elements are traceless. In order to separate the commutators from the time-dependence and hence the integral in Eq. (13), we expand the noise operators in this basis so that The expansion coefficientsB αk (t) ∈ R are given by the inner product of a noise operator in the interaction picture on the one hand and a basis element on the other: In line with Green et al. [29], we call these coefficients the control matrix (see also Refs. 65 and 35). In the transfer matrix (superoperator) picture we can take up the following interpretation for the control matrix by virtue of the cyclicity of the trace: it describes a mapping of a state, represented by the basis element C k and subject to the control operation U c (t) : C k → U c (t)C k U † c (t), onto the noise operator B α (t). That is, we can write the αth row of the control matrix as B α (t)| = B α (t)| U c (t). In this connection lies the power of the FF formalism as will become clear shortly; we can first determine the ideal evolution without noise and subsequently evaluate the error process by linking the unitary control operation to the noise operators.
Having expanded the noise operators in the basis C, we can already anticipate that upon substituting them, Eq. (13) will separate into a time-dependent part involving on one hand the control matrix and cross-correlation functions and on the other a time-independent part involving commutators of basis elements. This will simplify our calculations in the following. To see this, we recall the definition of the Liouville representation in Eq. (1) and apply it to the cumulant function so that K ij = tr(C i K[C j ]), where the notation K[C j ] means substituting C j for the placeholder · in the commutators in Eq. (13) and we suppressed the time argument for legibility. Finally, we insert the expanded noise operators given by Eq. (15) and obtain the Liouville representation of the cumulant function, Here, we captured the ordering of the noise operators due to the commutators in Eq. (13) in the coefficients f ijkl and g ijkl . These are trivial functions of the fourth order trace tensor [66] T given by Furthermore, we introduced the frequency (Lamb) shifts ∆ and decay amplitudes Γ which contain all information on the noise and qubit dynamics as captured by the control matrixB(t): The frequency shifts ∆ correspond to the first term in Eq. (12), hence incurring coherent errors, i.e. generalized axis and overrotation errors. They reflect a perturbative correction to the quantum evolution due to a change of the Hamiltonian at two points in time, and thus time ordering matters. Conversely, the decay amplitudes Γ correspond to the second term and capture the decoherence. These terms are due to an incoherent average that only takes classical correlations into account, so that time ordering does not play a role. Note that Eq. (17) together with Eq. (9) constitutes an exact version (in the Liouville representation) of Eq. (4) from Ref. 16 for Gaussian noise. The approximation of Ref. 16 is obtained by expanding the exponential to linear order and neglecting the second order terms ∆.
For a single qubit and C the Pauli basis one can make use of the simple commutation relations so that the cumulant function takes the form (see Appendix A) for i, j > 0 and any α, β. As mentioned in Section II A the cases j = 0 and i = 0 encode trace-preservation and unitality, respectively, and as such K 0j = K i0 = 0 since our model is both trace-preserving and unital.

B. Calculating the decay amplitudes
In order to evaluate the cumulant function K(τ ) given by Eq. (17) and thus the transfer matrix Ũ (τ ) from Eq. (9) for a given control operation, we solely require the decay amplitudes Γ kl and frequency shifts ∆ kl since the trace tensor T ijkl depends only on the choice of basis and is therefore trivial (although quite costly for large dimensions, c.f. Section III C) to calculate. In this section, we describe simple methods for calculating Γ kl using an extension of the filter function formalism developed by Green et al. [29] that we introduced in Ref. 16. The central quantity of interest will be the control matrix that we already introduced above. It relates the interaction picture noise operators to the operator basis and we will compute it in Fourier space in order to identify the crosscorrelation functions with the noise spectral density in Eq. (21). We distinguish between a sequence of quantum gates, as already presented in our related work [16], and a single gate. In the first case the control matrix of the entire sequence can be calculated from those of the individual gates, greatly simplifying the calculation if the latter have been precomputed. This approach gives rise to correlation terms in the expression for Γ kl that capture the effects of sequencing gates. In the second case, as was shown by Green et al. [29], one can calculate the control matrix for arbitrary single pulses under the assumption of piecewise constant control and we lay out how to adapt the approach for numerical applications.
We start by noting that, because we assumed the noise fields b α (t) to be wide-sense stationary, that is to say the cross-correlation functions evaluated at two different points in time t 1 and t 2 depend only on their difference t 1 − t 2 , we can define the two-sided noise power spectral density S αβ (ω) as the Fourier transform of the cross-correlation functions Note that the spectrum only characterizes the noise fully in the case of Gaussian noise. For non-Gaussian components in the noise, additional polyspectra have in principle to be considered for higher-order correlation functions [67]. However, since we only discuss second-order contributions which involve two-point correlation functions here, we only need to take S αβ (ω) into account. Inserting the definition of the spectral density into Eq. (21), one finds In the above equation, the fourth order tensor is the generalized filter function that captures the susceptibility of the decay amplitudes to noise at frequency ω. For α = β, k = l, and by summing over the basis elements, and this tensor reduces to the canonical fidelity filter function [28] from which the entanglement fidelity can be obtained, see Section II D 1. Thus, if the frequencydomain control matrixB αk (ω) for noise source α and basis element k is known, the transfer matrix can be evaluated by integrating Eq. (24). Moreover, one can study the contributions of each pair of noise sources (α, β) both separately or, at virtually no additional cost and to leading order, collectively by summing over them, We now discuss how to calculate the control matrix B(ω) in frequency space for a given control operation. We focus first on sequences of quantum gates, assuming that the control matricesB (g) (ω) for each gate g have been calculated before.

Control matrix of a gate sequence
For a sequence of gates with precomputed interaction picture noise operators the approach developed by Green et al. [29] based on piecewise constant control can be adapted to yield an analytical expression for those of the composite gate sequence that is computationally efficient to evaluate [16]. Here we review these results to give a complete picture of the formalism. While our results are general and apply to any superoperator representation, we employ the Liouville representation here for its simple composition operation: matrix multiplication. Computationally, this is not the most efficient choice since transfer matrices have dimension d 2 × d 2 and thus their matrix multiplication scales unfavorably compared to, for example, left-right conjugation by unitaries (c.f. Section III D). However, because the structure of the control matrixB is similar to that of a transfer matrix (remember that it corresponds to a basis expansion of the interaction picture noise operators), we will obtain a particularly concise expression for the sequence in the following. For a Figure 1. Illustration of a sequence of G gates. Individual gates with propagators Pg start at time tg−1 and complete at time tg. The total action from t0 to tg is given by Qg.
perhaps more intuitive description employing exclusively conjugation by unitaries, we refer the reader to our accompanying publication Ref. 16. Fig. 1 is considered. The cumulative propagator of the sequence up to time t g is then given by Q g = 0 g =g P g with P 0 = 1 and its Liouville representation denoted by Q (g) . Furthermore, the control matrix of the g-th pulse at the time t − t g−1 relative to the start of segment g is We can now exploit the fact that in the transfer matrix picture quantum operations compose by matrix multiplication to write the total control matrix at time The Fourier transform of Eq. (28) can then be obtained by evaluating the transform of each gate separately, with ∆t g = t g − t g−1 the duration of gate g. Hence, calculating the control matrix of the full sequence requires only the knowledge of the temporal positions, encoded in the phase factors e iωtg−1 , and the total intended action Q (g−1) of the individual pulses if their control matrices have been precomputed. The sequence structure can thus be exploited to one's benefit. If the same gates appear multiple times during the sequence one can reuse control matrices for equal pulses to facilitate calculating filter functions for complex sequences with modest computational effort. Most importantly, Eq. (29) is independent of the inner structure of the individual pulses and therefore takes the same time to evaluate whether they are highly complex or very simple. In Section III D, we will analyze the computational efficiency of capitalizing on this feature in more detail.
As we have seen, the total control matrix of a composite pulse sequence is given by a sum over the individual control matrices. SinceB(ω) enters Eq. (24) twice, this leads to correlation terms between two gates at different positions in the sequence when computing the total decay amplitudes Γ αβ,kl . Inserting Eq. (29) into Eq. (24) gives where we defined the pulse correlation filter function F (gg ) αβ,kl (ω) that captures the temporal correlations between pulses at different positions g and g in the sequence. Unlike regular filter functions, these can be negative for g = g and therefore reduce the overall noise susceptibility of a sequence given by F (ω) = gg F (gg ) (ω). We have thus gained a concise description of the noisecancelling properties of gate sequences: in this picture, they arise purely from the concatenation of different pulses, quantifying, for instance, the effectiveness of dynamical decoupling (DD) sequences [16].

Control matrix of a single gate
Previous efforts have derived the control matrix analytically for selected pulses such as dynamical decoupling (DD) sequences [20], special dynamically corrected gates (DCGs) [30], as well as developed a general analytic framework [28,29]. However, analytical solutions might not always be accessible, e.g. for numerically optimized pulses, and are generally laborious to obtain. Therefore, we now detail a method to obtain the control matrix numerically under the assumption of piecewise constant control. Our method is similar in spirit to that of Green et al. [28] for single qubits with d = 2, but whereas those authors computed analytical solutions to the relevant integrals during each time step, here we use matrix diago-nalization to obtain the propagator of a control operation to make the approach amenable to numerical implementation. This allows carrying out the Fourier transform of the control matrix Eq. (16) analytically by writing the control propagators in terms of their eigenvalues in diagonal form.
We divide the total duration of the control operation, τ , into G intervals (t g−1 , t g ] of duration ∆t g with g ∈ {0, . . . , G} and t 0 = 0, t G = τ . We then approximate the control Hamiltonian as constant within each interval so that within the g-th and similarly the deterministic time dependence of the noise operators as } contains the time evolution of the eigevalues. Using this result together with Q g−1 , the cumulative propagator up to time t g−1 , we can acquire the total time evolution operator at We then substitute this relation into the definition of the control matrix, Eq. (16), and obtaiñ where Ω Carrying out the Fourier transform of Eq. (34) to get the frequency-domain control matrix of the pulse generated by the Hamiltonian from Eq. (32) is now straightforward since the integrals involved are over simple exponential functions. We obtaiñ ij ) and the Hadamard product (A•B) ij := A ij ·B ij . Equation (35) is readily evaluated on a computer and thus enables the calculation of filter functions of arbitrary control sequences, either on its own or in conjunction with Eq. (29). A similar expression is obtained for representations other than the Liouville representation.

C. Calculating the frequency shifts
The frequency shifts ∆ αβ,kl in Eq. (17) correspond to the second order of the ME and thus involve a double integral with a nested time dependence. This makes their evaluation more involved than that of the decay amplitudes Γ αβ,kl and, in contrast to the previous section, we cannot identify a concatenation rule or single out correlation terms as in Eq. (31). However, we can still apply the approximation of piecewise constant control and follow a similar approach as in Section II B 2 to compute ∆ in Fourier space. Since these terms correspond to a coherent gate error that can in principle be calibrated out in experiments we will not go into much detail here.
We follow the arguments made above for the decay amplitudes and express the cross-correlation functions b α (t)b β (t ) by their Fourier transform, the spectral density S αβ (ω), using Eq. (23). Inserting this equation into the definition of the frequency shifts in the time domain, Eq. (20), yields We again assume piecewise constant time segments so that the inner time integral can be split up into a sum of integrals over complete constant segments (t g −1 , t g ] as well as a single integral that contains the last, incomplete segment up to time t. That is, taking the time t of the outer integral to be within the interval (t g−1 , t g ] we perform the replacement We have thus divided our task into two: The first term allows, as before in Sections II B 1 and II B 2, to identify the Fourier transform of the control matrix during time steps g and g for both the inner and the outer integral according to Eq. (35). The second term remains a nested double integral, but now the integrand contains only products of complex exponentials because we assume the control to be constant within the limits of integration. As a next step, we also replace the outer time integral by a sum of integrals over single segments, Before continuing, we ease notation and defineB(ω) =: g G (g) (ω) with G (g) (ω) obtained from Eq. (35) and furthermore adopt the Einstein summation convention for the remainder of this section, meaning multiple subscript indices that appear on only one side of an equality are summed over implicitly. We now proceed like in Section II B 2 and make use of the piecewise constant approximation to diagonalize the control Hamiltonian during each segment. For the nested integrals, we obtainB αk (t) from Eq. (34), whereas the remaining integrals factorize and we can identify the Fourier transformed quantity G (g) (ω). Equation (38) then becomes ij as defined above in Section II B 2 and Explicit results for the integration in Eq. (40) are given in Appendix A 2. To calculate the frequency shifts ∆, we can thus reuse the quantity G (g) (ω) also required for the decay amplitudes Γ. The only additional computation, apart from contraction, involves the G integrations I (g) ijmn (ω). Importantly, Eq. (39) has the same structure as the corresponding Eq. (24) for Γ in that the individual entries of ∆ are given by an integral over the spectral density of the noise multiplied with a -in this case second order -filter function that describes the susceptibility to noise at frequency ω:

D. Computing derived quantities
By means of Eqs. (29), (35) and (39), one can obtain the cumulant function K(τ ) from Eq. (17) and hence the error process Ũ (τ ) from Eq. (9) for an arbitrary sequence of gates. From this, several quantities of interest for the characterization of a given control operation can be extracted. We explicitly review the average gate and state fidelities as well as expressions to quantify leakage here, but emphasize that this is not exhaustive. Because for many applications the noise is weak and hence the parameter ξ 1, we will in the following expand the exponential in Eq. (7) to leading order in ξ in the following. That is, we approximate (remember that For Gaussian noise, higher order corrections can be obtained either by explicitly calculating higher powers of K or by numerically evaluating the exponential of the cumulant function. The former method often leads to simpler expressions than Eq. (17) for which the trace tensor T ijkl need not be computed directly. In the weak-noise regime, one can also define specific filter functions for each derived quantity that are given in terms of linear combinations of the generalized filter functions F αβ,kl (ω).
The ensemble expectation value of the quantity can then be obtained directly from the overlap with the spectral density, dω/2π F (ω)S(ω). Finally, we will drop the averaging brackets and the argument of the error transfer matrix Ũ (τ ) for brevity in the following.

Average gate and entanglement fidelity
The average gate fidelity is a commonly quoted figure of merit used to characterize physical gate implementations [5,8,[68][69][70]. It represents the fidelity between an implementation U and the ideal gate Q averaged over the uniform Haar measure. Since F avg (U, Q) = F avg (Q † • U, 1) = F avg (Ũ), the average gate fidelity can be obtained from the error channelŨ as [71,72] where d is the system dimension and F e (Ũ) = trŨ/d 2 is the entanglement fidelity. In the low-noise regime where Eq. (42) holds, we can write the entanglement fidelity in terms of the cumulant function K αβ approximately as Here, we defined I αβ , the infidelity due to a pair of noise sources (α, β). As we show in Appendix A 3, we can simplify the trace of the cumulant function so that the infidelity reads Equation (47) (24), we recover the relation (setting α = β for simplicity) with the fidelity filter function F α (ω) given by Eq. (26).
Notably, only the decay amplitudes Γ contribute to the fidelity to leading order since the frequency shifts ∆ are antisymmetric and therefore vanish under the trace.

State fidelity and measurements
In the context of quantum information processing we are often interested in the probability of measuring the expected state during readout. We can extract this projective readout probability from the transfer matrix in Eq. (7) by inspecting the transition probability, or state fidelity, between a pure state ρ = |ψ ψ| and an arbitrary state σ that evolves according to the quantum operation E : σ → E(σ). Using the double braket notation introduced at the beginning of Section II we then define the state fidelity as where we have expressed the density matrices by vectors on the Liouville space L and E as a transfer matrix. We can thus calculate arbitrary pure state fidelities by simple matrix-vector multiplications of the transfer matrices E = QŨ and the vectorized density matrices |ρ and |σ .
In Section V C we employ this measure to simulate a RB experiment where return probabilities are of interest so that F(|ψ , E(ρ)) = ρ| QŨ |ρ . General measurements can be incorporated in the superoperator formalism we have employed here in a straightforward manner using the positive operatorvalued measure (POVM) formalism [49,73]. POVMs constitute a set of Hermitian, positive semidefinite operators {E i } i (in contrast to the projective measurement {|ψ ψ| , 1 − |ψ ψ|}) that fulfill the completeness relation i E i = 1 and in the double braket notation may be represented as the row vectors { E i | } i in Liouville space. Consequently, the measurement probability for outcome E i is given by E i |E(σ) = E i |E|σ if the system was prepared in the state σ and evolved according to E.

Leakage
In many physical implementations qubits are not encoded in real two-level systems but in two levels of a larger Hilbert space (e.g. transmon [74] or singlet-triplet [75] spin qubits) such that population can leak between this computational subspace and other energy levels. Thus, it is often of interest to quantify leakage when assessing gate performance. Recently, Wood and Gambetta [76] have suggested two separate measures for quantifying leakage out of the computational subspace on the one hand and seepage into the subspace on the other. With the filter function formalism and the transfer matrix of the error process given by Eqs. (7) and (17), we can easily extract these quantities.
Using the definitions from Ref. 76 and the double braket notation we can write the leakage rate generated by a quantum operation E as and the seepage rate as Here, Π c, are projectors onto the computational and leakage subspaces, respectively, and d c, the corresponding dimensions. For unital channels the leakage and seepage rates are not independent but satisfy d c L c = d L [76] so that we only need to consider one of the above expressions here (c.f. Section II A 2). Equations (50a) and (50b) can be used to determine both coherent and incoherent leakage separately by substituting Q orŨ, respectively, for E. While the former is due to systematic errors of the applied pulse and could thus be corrected for by calibration, the latter is induced by noise only. Alternatively, the leakage from both contributions can also be determined collectively by substituting U for E.

III. PERFORMANCE ANALYSIS AND EFFICIENCY IMPROVEMENTS
In this section we focus on computational aspects of the formalism, remarking first on several mathematical simplifications that make the calculation of control matrices and decay amplitudes more economical. Following this, we investigate the computational complexity of the method in comparison with Monte Carlo techniques and show that our software implementation surpasses the latter's performance in relevant parameter regimes.

A. Periodic Hamiltonians
If the control Hamiltonian is periodic, that is H c (t) = H c (t + T ), we can reduce the computational effort of calculating the control matrix by potentially orders of magnitude (see Section V B for an application in Rabi driving). We start by making the following observations: First, the frequency domain control matrix of every period of the control is the same so thatB (g) (ω) =B (1) (ω).
Moreover, e iω∆tg = e iωT for all g so that e iωtg−1 = e iωT (g−1) and by the composition property of transfer matrices Q (g−1) = Q (1) g−1 where the superscript without parentheses denotes matrix power. We can then simplify Eq. (29) to read which is typically the case for the vast majority of values of ω, the previous expression can be rewritten as by evaluating the sum as a finite Neumann series. Equation (52) offers a significant performance benefit over regular concatenation in the case of many periods G as we will show in Section III D. Beyond numerical advantages, it also provides an analytic method for studying filter functions of periodic driving Hamiltonians.

B. Extending Hilbert spaces
Examining Eq. (16), we can see that the columns of the control matrix and therefore also the filter function are invariant (up to normalization) under an extension of the Hilbert space. This allows parallelizing pulses with precomputed control matrices in a very resource-efficient manner if one chooses a suitable operator basis. Note that the same also applies to other representations of quantum operations.
Suppose we extend the Hilbert space H 1 of a gate for which we have already computed the control matrix by a second Hilbert space H 2 so that H 12 = H 1 ⊗ H 2 . If we can find an operator basis whose elements separate into tensor products themselves, i.e. C 12 = C 1 ⊗ C 2 as for the Pauli basis (c.f. Section III C), the control matrix of the composite gate defined on H 12 has the same nontrivial columns as that of the original gate on H 1 up to a different normalization factor. The remaining columns are simply zero. This is because the trace over a tensor product factors into traces over the individual subsystems so thatB since we assumed that the noise operators B α are traceless (c.f. Section II A 2).
Generalizing this result to multiple originally disjoint Hilbert spaces we write the composite space as H = i H i and the corresponding basis as C = i C i . The control matrix of the composite pulse on H is then a combination of the columns of the control matrices on H i for noise operators B α that are non-trivial, i.e. not the identity, only on their original space. For noise operators defined on more than one subspace, e.g. B ij = B i ⊗ B j , B i ∈ H i , B j ∈ H j , this does not hold anymore and the corresponding row in the composite control matrix needs to be computed from scratch.
One can thus reuse precomputed control matrices beyond the concatenation laid out above when studying multi-qubit pulses or algorithms. For concreteness, consider a set of one-and two-qubit pulses whose control matrices have been precomputed. We can then remap those control matrices to any other qubit in a larger register if the entire Hilbert space is defined by the tensor product of the single-qubit Hilbert spaces, and even map the control matrices of two different pulses to the same time slot on different qubits. Thus, we do not need to perform the possibly costly computation of the control matrices again but instead only need to remap the columns ofB to the equivalent basis elements in the basis of the complete Hilbert space, making the assembly of algorithms that consist of a limited set of gates which are used at several points in the algorithm more efficient. In Section V D we simulate a four-qubit QFT algorithm making use of the shortcuts described here.

C. Operator bases
Up to this point, we have not explicitly specified the basis that defines the Liouville representation. The only conditions imposed by Eq. (14) are orthonormality with respect to the Hilbert-Schmidt product and that the basis elements are Hermitian. Yet, the choice of operator basis can have a large impact on the time it takes to compute the control matrix as discussed in the previous section. We therefore give a short overview over two possible choices in the following. As we are mostly interested in the computational properties, we represent linear operators in L(H ) as matrices on C d×d .
The n-qubit Pauli basis fulfills the requirements set by Eq. (14) and furthermore allows for the simplifications described before. In our normalization convention where C i , C i = 1 it can be written as with the Pauli matrices σ x , σ y and σ z . While it is obvious that it is separable, meaning it factors into tensor products of the single-qubit Pauli matrices, the dimension of the Pauli basis is restricted to powers of two, i.e. d = 2 n . An operator basis without this restriction is the generalized Gell-Mann (GGM) basis [77,78]. In the following we will discuss optimizations pertaining to this basis that are also implemented in the software (see Section IV). The GGM matrices are a generalization of the Gell-Mann matrices known from particle physics to arbitrary dimensions. In our normalization convention, the basis (excluding the identity element) is given by [79] with and an orthonormal vector basis {|j } d j=1 of the Hilbert space. Expanding an arbitrary matrix A ∈ C d×d in the basis of Eq. (54) is then simply a matter of adding up the corresponding matrix elements of A according to Eqs. (54a) to (54c). For instance, the expansion coefficient for the first symmetric basis element is given by u 12 = (A 12 |1 2| + A 21 |2 1|) √ 2 . The explicit construction prescription of the GGM basis thus allows calculating inner products of the form Λ j |A at constant cost instead of the quadratic cost of the trace of a matrix product, speeding up the computation of the transfer matrix from Eq. (1) (in which case A = E(Λ k )). In numerical experiments, calculating the transfer matrix of a unitary U with dimension d and precomputed matrix products A k = U Λ k U † scaled as ∼ d 4. 16 . This agrees with the expected scaling of ∼ d 4 (a transfer matrix has d 2 × d 2 elements) and presents a significant improvement over the explicit calculation with trace overlaps tr(Λ j A k ) that we observed to scale as ∼ d 5.93 (we expected ∼ d 6 ).
Further inspection of the GGM basis additionally reveals an increasing sparsity for large d (the filling factor scales roughly with d −2 ), so that it is well suited for computing the trace tensor Eq. (18). Since this tensor has d 8 elements, the amount of memory required for a dense representation becomes unreasonably large quite quickly. To overcome this constraint, we can use a GGM basis instead of a dense basis like the Pauli basis (which has a filling factor of 1/2 ). In this case, the resulting tensor is also sparse because the overlap between different basis elements is small. This not only enables storing the tensor in memory but also makes the calculation much faster since one can employ algorithms optimized for operations on sparse arrays (see Section IV).
As an illustration, consider a system of four qubits so that the Hilbert space has dimension d = 2 4 . An operator basis for this space has d 2 = 2 8 elements and consequently the tensor T ijkl has (2 8 ) 4 = 2 32 entries. Using 128 bit complex floats to represent the entries the tensor would take up ≈ 68 GB of memory in a dense format. Conversely, for a GGM basis stored in a sparse data structure, the resulting trace tensor only takes up ≈ 100 MB of memory. Furthermore, calculating T takes ≈ 2.89 s on an Intel ® Core™ i9-9900K eight-core processor since a GGM has a low filling factor. By contrast, the same calculation with a Pauli basis takes ≈ 217 s. This is due to the larger filling factor on the one hand and because sparse matrix multiplication algorithms perform poorly with dense matrices on the other.

D. Computational complexity
In order to assess the performance of filter functions (FF) for computing fidelities compared to Monte Carlo (MC) methods, we determine each method's asymptotic scaling behavior as a function of the system dimension d. For the filter functions, we calculate the fidelities using Eqs. (26) and (48) in our software implementation, described in more detail in Section IV and hence neglect contributions of O ξ 4 from the frequency shifts ∆. Additionally, we distinguish between three different approaches for calculating the control matrix; first, for a single pulse following Eq. (35), second for an arbitrary sequence of pulses following Eq. (29), and third for a periodic sequence of pulses following Eq. (52). For the single pulse, we run benchmarks using exemplary values for the various parameters on a machine with an AMD FX™-6300 processor with six logical cores and 24 GB of memory. We also discuss the filter function method using left-right conjugation by unitaries instead of the Liouville representation. The latter has higher memory requirements and is expected to perform poorly for large system dimensions d since one deals with d 2 × d 2 transfer matrices on a Liouville space L instead of d × d unitaries on a Hilbert space H . In the software package, the calculations are currently implemented in Liouville space and calculation by conjugation is only partially supported through the low-level API. However, both representations perform similarly for small dimensions as we show below. Note that for a fair performance comparison the different nature of errors needs to be kept in mind. Monte Carlo becomes less costly if larger statistical errors can be tolerated, whereas the filter function formalism is typically limited by higher order errors. For reference, the following considerations are summarized in Table I for each approach and a representative set of parameters.
To calculate the fidelity using MC, we generate n MC different noise traces that slice every time step ∆t of the pulse into n seg = f UV ∆t segments to appropriately sample the spectral density with f UV being the ultraviolet cutoff frequency. In total, there are n ∆t n MC n seg . For simplicity, we use a white noise spectrum for which S(ω) = const. but note that sampling arbitrary spectra induces additional overhead for MC, depending on which method is used to generate the noise traces. Typical time-domain methods include the simulation of the underlying physical process (like two-state fluctuators) or the application of an inverse Fourier transform to white noise multiplied by a frequency-domain transfer function. By contrast, the computational cost of the filter function formalism as realized by Eq. (35) is independent of the form of the spectrum. For this approach we find the leading terms to scale as ∼ n ∆t n ω n α d 4 +n ∆t d b+2 with n α the number of noise operators and n ω the number of frequency samples. Here, the first term is due to the trace in Eq. (35) which boils down to the trace of a matrix product, ij A ij B ji , that scales with d 2 and is performed once for each of the d 2 basis elements, n α noise operators, n ∆t time steps, and n ω frequency points. The second term is due to the transformation C k →C (g) k which requires multiplication of d × d matrices for every time step and basis element. As n α n ω < n MC n seg for realistic parameters because the ultraviolet cutoff frequency needs to be chosen sufficiently high and the relative error of the method decreases with 1 √ n MC , we expect that in the case of a single pulse the filter function formalism in Liouville representation should outperform Monte Carlo calculations for reasonably small dimensions d. Using left-right conjugation, this advantage should hold also for large d. In this case the Hadamard product (∼ d 2 ) as well as matrix multiplication (∼ d b ) are carried out for each frequency, noise operator, and time step to calculate the interaction picture noise operatorsB α (ω). We thus find this method to scale with ∼ n ∆t n ω n α (d 2 + d b ). Figure 2 shows exemplary wall times for both methods and d ∈ [2, 120] that confirm our expectation. Only for about d ≈ 100 the overhead from the extra time steps and trajectories over which is averaged is compensated for MC. For smaller dimensions the calculation using FFs is faster by almost two orders of magnitude (see the inset showing the same data in a log-log plot). The lines show fits to t = ad b . The data is not quite in the asymptotic regime due to limited memory so that even for large dimension terms of lower power in d contribute significantly to the run time. Even though this causes the fits to underestimate the exponent b, the general trend agrees with our expectation. Note that the crossover does not always occur at the same dimension d. On a different system with an Intel ® Core™ i9-9900K eight-core processor the FF method outperformed MC even for d = 120 beyond which available memory limited the simulation.
Quantifying the performance gain from using the control matrices' concatenation property to calculate fidelities of gate sequences is more difficult since it strongly depends on the number of gates occurring multiple times in the sequence (enabling reuse of precomputed control matrices) as well as the complexity of the individual gates. The evaluation using the concatenation rule Eq. (29) performs asymptotically worse than the evaluation for a complete pulse according to Eq. (35) because of higher powers of d dominating the calculation in the former case. Performing the G matrix multiplica-tionsB (g) (ω)Q (g−1) from Eq. (29) is of order ∼ Gn ω n α d 4 , with G the number of pulses in the sequence. Furthermore, calculating the transfer matrix of the total propagators Q g−1 involves multiplication of d × d matrices for all d 4 combinations of basis elements amounting to ∼ Gd b+4 . In case the Liouville representation of the individual pulses' total propagators P g , P (g) , have been precomputed, the latter computation can be made more efficient since one can just propagate the transfer matrices P (g) to obtain the cumulative transfer matrices for the sequence, Q (g) = 0 g =g Q (g ) , at cost ∼ Gd 2b . The restriction to small dimensions does not apply for conjugation by unitaries as in this case the matrix multiplications involve d × d matrices and we do not have to compute the Liouville representation. We thus obtain a more favorable asymptotic scaling of ∼ Gn ω n α d b .
Utilizing the concatenation property in the Liouville representation thus corresponds to effectively reducing the number of times the calculations scaling with ∼ n ω n α d 4 have to be carried out but incurs additional cal- Carlo trajectories over which is averaged, and nω the number of frequency samples. The calculation using filter functions clearly outperforms MC for small system sizes. For dimensions larger than d ≈ 100 (roughly equivalent to 7 qubits) Monte Carlo (blue squares) performs better than the FF calculation with transfer matrices (green triangles) for this set of parameters and processor due to the better scaling behavior. Using conjugation by unitaries (orange diamonds) significantly outperforms MC also for large dimensions. While the fits to t = ad b (lines) underestimate the leading order exponent due to the data not being in the asymptotic regime, they support the expected relationship of complexity between the approaches. The inset shows the same data on a linear scale, highlighting the different scaling behaviors for large d.
culations scaling with ∼ d 2b . Accordingly, it provides a performance benefit if a sequence consists of either very complex pulses, in which case single repetitions already make the calculation much more efficient, or of few pulses that occur many times. In the extremal case of G repetitions of a single gate the benefit of employing the concatenation property is most pronounced and can be improved even further utilizing the simplifications laid out in Section III A. Since matrix inversion has the same complexity as matrix multiplication and taking a matrix to the G-th power requires O(log G) matrix multiplications, Eq. (52) should scale with ∼ n ω (n α d 4 + d 2b + d 2b log G) (the first two terms are due to the final matrix multiplications and are independent of G). It hence allows for a vast speedup over Eq. (29) in that the asymptotic behavior as a function of the number of gates changes from ∼ G to ∼ log G. An example of this is presented in Section V B for the context of Rabi driving. Note that this closed form is a unique feature of the transfer matrix representation and not applicable to conjugation by unitaries.

IV. SOFTWARE IMPLEMENTATION
In this section we give an overview over the filter_functions software package implementing the main features of the formalism derived above. This includes the calculation of the decay amplitudes Γ and fidelities as well as the calculation of the control matrices for single gates and both generic and periodic sequences of gates. Moreover, control matrices may be efficiently extended to and merged on larger Hilbert spaces. Calculations using unitary conjugation instead of transfer matrices are implemented but at this point not available in the high-level API.
Our software is written in Python and available on GitHub [41] under the GPLv3 license. We also provide a current snapshot in the Supplementary Material [81]. It features a broad coverage through unit tests and extensive API documentation as well as didactic examples (see Section V). The package relies on the NumPy [82] and SciPy [83] libraries for vectorized array operations. Data visualization is handled by matplotlib [84]. For tensor multiplications with optimized contraction order we use opt_einsum [85] for which sparse [86], a library aiming to extend the SciPy sparse module to multi-dimensional arrays, serves as a backend in the calculation of the trace tensor from Eq. (18). Lastly, the package is written to interface with qopt [43,87] and QuTiP [42], frameworks for the simulation and optimization of open quantum systems, and mirrors the latter's data structure for Hamiltonians ensuring easy interoperability between the two.

A. Package overview
In the filter_functions package all operations are understood as sequences of pulses that are applied to a quantum system. These pulses are represented by instances of the PulseSequence class which holds information about the physical system (control and noise Hamiltonians) as well as the mathematical description (e.g. the basis used for the Liouville representation). As indicated above, the Hamiltonians H c (t) and H n (t) are given in a similar structure as in QuTiP. That is, a Hamiltonian is expressed as a sum of Hermitian operators with the time dependence encoded in piecewise constant coefficients so that for t ∈ (t g−1 , t g ], g ∈ {1, . . . , G} and where the a i ]. The PulseSequence class provides methods to calculate and cache the filter function according to Eq. (30). Alternatively, filter functions may also be cached manually to permit using the package with analytical solutions for the control matrix. Concatenation of pulses is implemented by the functions concatenate() and concatenate_periodic() which will attempt to use the cached attributes of the PulseSequence instances representing the pulses to efficiently calculate the filter function of the composite pulse following Eq. (29) and Eq. (52), respectively.
Operator bases fulfilling Eq. (14) are implemented by the Basis class. There are two predefined types of bases: The Pauli basis is both unitary and separable while the GGM basis is sparse for large dimensions but neither unitary nor separable. As mentioned in Section III B (see also Section V D), using a separable basis can provide significant performance benefits for calculating the filter functions of algorithms. On the other hand, a sparse basis makes the calculation of the trace tensor T ijkl and therefore also of the error transfer matrixŨ much faster (c.f. Section III C). Additionally, the user can define custom bases using the class constructor. The error transfer matrixŨ can be calculated for a given pulse and noise spectrum using the error_transfer_matrix() function [88]. Various other quantities can be computed fromŨ as outlined in Section II D. Furthermore, the package includes a plotting module that offers several functions, e.g. for the visualization of filter functions or the evolution of the Bloch vector using QuTiP.

B. Workflow
We now give a short introduction into the workflow of the filter_functions package by showing how to calculate the dephasing filter function of a simple Hahn spin echo sequence [89] as an example. The sequence consists of a single π-pulse of finite duration around the x-axis of the Bloch sphere in between two periods of free evolution. We can hence divide the control fields into three constant segments and write the control Hamiltonian as with τ the duration of the free evolution period and t π that of the π pulse. For the noise Hamiltonian we only need to define the deterministic time dependence s α (t) and operators B α since the noise strength is captured by the spectrum S(ω). Thus we have s z (t) = 1 and B z = σ z /2 for pure dephasing noise that couples linearly to the system. In the software, we first define a PulseSequence object representing the spin echo (SE) sequence. As was already mentioned, the control and noise Hamiltonians are given as a list containing lists for every control or noise operator that is considered. These sublists consist of the respective operator as a NumPy array or QuTiP Qobj and the amplitudes (a (g) i or s (g) α ) in an iterable data structure such as a list. We can hence instantiate the PulseSequence with the following code: import filter_functions as ff import qutip as qt from math import pi tau, t_pi = (1, 1e-3 where a basis is automatically chosen since we did not specify it in the constructor in the last line. Calculating the filter function of the pulse for a given frequency vector omega can then be achieved by calling where F is the dephasing filter function F zz (ω) as we only defined a single noise operator. Finally, we calculate the error transfer matrixŨ for the noise spectral density S zz (ω) = ω −2 , S = 1/omega**2 U = ff.error_transfer_matrix (ECHO, S, omega) This code uses the control matrix previously cached when the filter function was first computed. Therefore, only the integration in Eq. (24) and the calculation of the trace tensor in Eq. (18) are carried out in the last line.
An alternative approach to calculate the spin echo filter function is to employ the concatenation property. For this, we interpret the SE as a sequence consisting of three separate pulses. Each of the pulses has a single time segment during which a constant control is applied and concatenating the separate PulseSequence instances yields the PulseSequence representing a spin echo. This way analytic control matrices may be used to calculate the control matrix of the composite sequence. Pulses can be concatenated by using either the concatenate() function or the overloaded @ operator: Since we have cached the control matrices of the FID and NOT pulses, that of the ECHO object is also automatically calculated and stored. Concatenating PulseSequence objects is implemented as an arithmetic operator of the class to reflect the intrinsic composition property of the control matrices.
Further development of the software has focused on making it available in a gate optimization and simulation framework [43,87]. Besides using it to compute decoherence effects and fidelities, analytic derivatives of the filter functions have been implemented to allow for optimizing pulse parameters in the presence of non-Markovian noise within the framework of quantum optimal control [90].
Additionally, building an interface with qupulse [91,92], a software toolkit for parametrizing and sequencing control pulses and relaying them to control hardware, would implement the capability to compute filter functions of pulses assembled in qupulse, thereby allowing a user in the lab to easily inspect the noise susceptibility characteristics of the pulse they are currently applying to their device.

V. EXAMPLE APPLICATIONS
We now present example applications of the software package and the formalism. As stated before, we focus on the decay amplitudes Γ and its associated filter functions and assume that the unitary errors generated by the frequency shifts ∆ are either small (as is the case for gate fidelities) or calibrated out. All of the examples shown below are part of the software documentation as either interactive Jupyter notebooks [93] or Python scripts. In the following, we give angular frequencies and energies in units of inverse times (e.g. s −1 ) while ordinary frequencies are given in Hz and we write Ũ (τ ) =Ũ for legibility.

A. Singlet-triplet two-qubit gates
In order to benchmark fidelity predictions of our implementation as well as demonstrate its application to nontrivial pulses, we compute the first-order infidelity of the two-qubit gates presented in Ref. 94 and compare the results to the reference's Monte Carlo calculations. There, a numerically optimized gate set consisting of {X π/2 ⊗I, Y π/2 ⊗I, CNOT} for exchange-coupled singlettriplet spin qubits is introduced, taking into account different noise spectra and realistic control hardware.
For readers unfamiliar with the reference we briefly summarize the physical system and noise model entering the optimization. The authors consider four electrons confined in a linear array of four quantum dots in a semiconductor heterostructure. Each electron i ∈ {1, 2, 3, 4} experiences a different static magnetic field B i so that there is a gradient b ij = B i − B j between two adjacent dots i and j. This gives rise to spin quantization along the magnetic field axis and defines the eigenstates {|↑↓ , |↓↑ } that span the computational subspace of a single qubit so that the accessible Hilbert space of the two-qubit system is spanned by {|↑↓ , |↓↑ } ⊗2 . The magnitude of the exchange interaction J ij between two adjacent dots i and j is controlled via gate electrodes located on top of the heterostructure that can be pulsed on a nanosecond timescale with an arbitrary waveform generator (AWG). Changing the gate voltages changes the detuning ij of the electrochemical potential between dots and in turn leads to a change in exchange coupling according to the The pulses are defined by a set of discrete detuning voltages ij passed to an AWG with a sample rate of 1 GS/s and constant magnetic field gradients b ij are assumed. To reflect the fact that the qubits experience a different pulse than what is programmed into the AWG due to cable dispersion and non-ideal control hardware, the detunings are convoluted with an experimental impulse response [94]. Finally, the signal is discretized as piecewise constant by slicing each segment into five steps, yielding a time increment of ∆t = 0.2 ns.
To find optimal detuning pulses, a Levenberg-Marquardt algorithm iteratively minimizes the infidelity, leakage, and trace distance from the target unitary. For the infidelity, contributions from quasistatic magnetic field noise as well as quasistatic and white charge noise are taken into account during each iteration. Because treating colored (correlated) noise using Monte Carlo methods is computationally expensive (c.f. Section III D), the infidelity due to fast 1/f -like noise is only computed for the final gate and not used during the optimization.
Two-qubit interactions are mediated via the exchange J 23 that makes the states |↑↑↓↓ and |↓↓↑↑ accessible. They constitute levels outside of the computational subspace that ideally should only be occupied during an entangling gate operation. A non-vanishing population of these states after the operation has ended is therefore unwanted and considered leakage, the magnitude of which we could quantify following Section II D 3. However, here we limit ourselves to determine the infidelity contribution from fast, viz. non-quasistatic, charge noise entering the system through ij . That is, we consider noise sources α ∈ { 12 , 23 , 34 }. We take the non-linear dependence of the Hamiltonian on the detunings ij into account by setting s ij (t) = ∂J ij ( ij (t))/∂ ij (t) ∝ J ij ( ij (t)). Figure 3 shows the filtered (convoluted) exchange interaction J ij between each pair of dots during the pulse sequence in panel (a) and filter functions plotted as function of frequency in panel (b) for the three different detunings. For a detailed description on how the filter functions were computed in the presence of additional leakage levels refer to Appendix B. As one would expect from the fact that the intermediate (inter-qubit) exchange interaction J 23 (orange dash-dotted lines) is only turned on for short times to entangle the qubits, the filter function for 23 is smaller by roughly an order of magnitude than the intra-qubit exchange filter functions. Notably, the filter functions for 12 and 34 show clear characteristics of DCGs, that is they drop to zero as ω → 0, and decouple from quasistatic noise with an error suppression ∝ ω 2 . This is not unexpected as the optimization minimizes quasistatic noise contributions to the infidelity. In addition, one can also observe small oscillations with period 5 ns −1 in frequency space that arise as a numerical artifact of the piecewise constant discretization of the control parameters as investigations have shown. If high-frequency spectral components are expected to play a significant role, one needs to be aware of these effects and adjust the simulation parameters appropriately.
The inset of Fig. 3(b) shows the same filter functions for the DC tail on a linear scale. Most notably, F 12 and F 34 have maxima around ω = 2π/τ , i.e. exactly the frequency matching the pulse duration, and around ω = 50/τ = 1 ns −1 with τ CNOT = 50 ns. The former is the typical window in which a pulse is most susceptible to noise whereas the latter matches the absolute value of the magnetic field gradients, b 12 = −b 34 = 1 ns −1 , indicating that the peak corresponds to the qubit dynamics generated by the magnetic field gradients. Panels (c)-(e) show the cumulant functions K ij (τ ) of the detuning error channels ij on the computational subspace. K 12 displays clear characteristics of a Pauli channel with only elements on the diagonal and secondary diagonals deviating from zero significantly whereas K 34 (the target qubit) possesses a more complicated structure.
We now compute the infidelity contribution originating from fast charge noise using Eq. (43) but tracing only over the computational subspace to compare to the Monte Carlo calculations of Ref. 94 (see Appendix B for further details). Like the reference, we use a noise spectrum S ,a (f ) ∝ 1/f a with S ,a (1 MHz) = 4 × 10 −20 V 2 Hz −1 and consider white noise (a = 0) and correlated noise with a = 0.7 [10] with infrared and ultraviolet cutoffs 1/τ and 100 ns −1 , respectively. Table II compares the results in this work with the reference. The values computed here are consistent with the more elaborate Monte Carlo calculations within a few percent. Notably, the deviation is smaller for the smaller noise levels with a = 0.7, in line with the fact that we have only computed the contributions from the decay amplitudes Γ and thus the leading order perturbation. If we had additionally evaluated the frequency shifts ∆ we could have obtained the exact fidelity in the case of Gaussian noise.

B. Rabi driving
A widely used method for qubit control is Rabi driving [8,13,26,95]. If we restrict ourselves to the resonant case for simplicity, the control Hamiltonian takes on the general form H c = ω 0 σ z /2 + A sin(ω 0 t + φ)σ x . Here, ω 0 is the resonance frequency, A the drive amplitude corresponding to the Rabi frequency in the weak driving limit A/ω 0 1, Ω R ≈ A, and φ an adjustable phase giving control over the rotation axis in the xy-plane of X π/2 ⊗ I 1.7 × 10 −3 5.8 × 10 −5 1.9 × 10 −3 5.7 × 10 −5 Y π/2 ⊗ I 1.6 × 10 −3 5.7 × 10 −5 1.7 × 10 −3 5.6 × 10 −5 CNOT 1.5 × 10 −3 6.4 × 10 −5 1.6 × 10 −3 6.3 × 10 −5 the Bloch sphere. This Hamiltonian and associated decoherence mechanisms are well-studied in the weak driving regime, where the rotating wave approximation (RWA) can be applied to remove fast-oscillating terms in the rotating frame [96,97]. There is a comprehensive understanding of how spectral densities transform to this frame and which frequencies are most relevant to loss of coherence [98]. By contrast, the description of a system in the strong driving regime, where A/ω 0 ∼ 1, is more complicated since the RWA cannot be applied without making large errors. Yet, an improved understanding is desirable because strong driving allows for much shorter gate times and thus shifts the window of relevant noise frequencies towards higher energies where the total noise power is typically lower, e.g. for 1/f noise. Conversely, faster control also requires more accurate timing to prevent rotation errors. It is therefore of interest to have available tools that can provide a comprehensive picture for Rabi pulses over a wide range of driving amplitudes. By making use of the concatenation property of the filter functions, our formalism can do just that.
The problem that arises when trying to numerically investigate Rabi pulses in the weak driving regime in the lab frame is that typical control operations have a duration τ T with T = 2π/ω 0 . Since the sampling time step ∆t should additionally be chosen much smaller than a single drive period in order to sample the time evolution accurately (∆t T ), brute-force simulations are costly.
For T /∆t = 100 samples per period and assuming Rabi and drive frequencies in typical regimes for SiGe and MOS quantum dots [99,100] or trapped ions [26], Ω R = 1 µs −1 and ω 0 = 20 ns −1 , a Monte Carlo simulation of a π-rotation with approximately 3 % relative error would require 10 9 samples in total. Using the filter function formalism, we can drastically reduce the simulation time even beyond the improvement gained from concatenating precomputed filter functions of individual drive periods using Eq. (29). This can be achieved with Eq. (52), which simplifies the calculation of the control matrix for periodic Hamiltonians.
To benchmark our implementation, we use the param-   eters from above and calculate the control matrix of a NOT gate generated by a Rabi Hamiltonian with three different methods on an Intel ® Core™ i9-9900K eightcore processor. First, we use Eq. (29) in a brute force approach. Second, we utilize the concatenation property following Eq. (35). Third, we employ the simplified expression given by Eq. (52). The brute force approach takes 250 s of wall time whereas calculating the filter function using the standard concatenation is faster by two orders of magnitude, taking 1.5 s to run. Lastly, the calculation utilizing the optimized method is faster again by two orders of magnitude and is completed in 0.056 s.
As an example application, we calculate the filter functions for continuous Rabi driving in the weak and strong driving regimes. For weak driving, we use the parameters from the benchmark above for a pulse of duration τ weak ≈ 20 µs that corresponds to 20 identity rotations in total. For the strong driving regime, we use the approximate analytical solution for a flux qubit biased at its symmetry point from Ref. 101 with A = ω 0 /4 to drive the qubit for τ strong ≈ 4 ns so that we achieve the same amount of identity rotations as in the weak driving case. In the reference, strong driving in this regime is shown to give rise to non-negligible counterrotating terms that modulate the Rabi oscillations and which are welldescribed by Floquet theory applied to the Rabi driving Hamiltonian. While for the regime studied here only two additional modes appear, the results extend to the regime where A > ω 0 and up to eight different frequency components were observed. Figure 4 shows the filter functions F xx and F zz for the σ x and σ z noise operators in the weak (a) and the strong (b) driving regime. Both display sharp peaks at their Rabi frequencies and the resonance frequency for F zz and F xx , respectively. We expect these features as they correspond to perturbations of the qubit Hamiltonian that are resonant with the qubit dynamics about an axis orthogonal to them. For weak driving, F xx is constant up to the resonance frequency where it peaks sharply and then aligns with F zz . The latter has a peak at the Rabi frequency before rolling off with ω −2 and a DC level that is almost ten orders of magnitude larger than that of the transverse filter function. This behavior is consistent with the results by Yan et al. [98], who show that the noise sources dominating decoherence during driven evolution are S xx (ω 0 ) and S zz (Ω R ). Note that the piecewise constant control approximation causes the weak driving filter functions to level off towards low fre-  (a) Weak driving with A/ω0 1. The filter function Fxx for noise operator σx is approximately constant up to the resonance frequency where it peaks sharply and then aligns with the filter function Fzz for σz. Fzz peaks at the Rabi frequency before rolling off with ω −2 and a DC level that is almost ten orders of magnitude larger than the DC level of the transverse filter function Fxx. (b) Strong driving with A/ω0 ∼ 1. Again Fzz peaks at ΩR whereas Fxx has three distinct peaks at ω0 and ω0 ± ΩR. These features also appear at slightly higher frequencies in Fzz due to the strong coupling.
quencies after an initial roll-off (here at ω ∼ 1 ms −1 ). By decreasing the discretization time step ∆t, one can shift the frequency at which this effect occurs to lower frequencies and thus attribute the feature to a numerical artefact of the approximation. However, the decoupling properties depend quite sensitively on the pulse duration.
In case of strong driving, the two filter functions are closer in amplitude for lower frequencies. In addition, F xx also peaks at ω = ω 0 ± Ω R . These peaks also show up at higher frequencies in the dephasing filter function F zz , reflecting frequency mixing in the strong coupling regime. While both filter functions show characteristics of a DCG in the weak driving regime, that is they drop to zero as ω → 0, this is not the case in the strong driving regime. Instead, there they approach a constant level for small frequencies. On top of rotation errors from timing inaccuracies, we may thus expect naive strong driving gates to be more susceptible to quasistatic noise than weak driving gates. By shaping the pulse envelope of the strong driving gate the decoupling properties could be recovered.

C. Randomized Benchmarking
Standard Randomized Benchmarking (SRB) and related methods, for example interleaved RB, are popular tools to assess the quality of a qubit system and the operations used to control it [50,102,103]. The basic protocol consists of constructing K random sequences of varying length m of gates drawn from the Clifford group [104], and appending a final inversion gate so that the identity operation should be performed in total. Each of these pulse sequences is applied to an initial state |ψ in order to measure the survival probability p(|ψ ) after the sequence. In reality, the applied operations are subject to noise and experimental imprecisions. This renders them imperfect and results in a survival probability smaller than one. Assuming gate-independent errors, the average gate fidelity F avg is then obtained by fitting the measured survival probabilities for each sequence length to the zeroth-order exponential model [103] where r = 1 − F avg is the average error per single gate to be extracted from the fit, A and B are parameters capturing state preparation and measurement (SPAM) errors, and d is the dimensionality of the system. One of the main assumptions of the SRB protocol is that temporal correlations of the noise are small on timescales longer than the average gate time [103]. If this requirement is not satisfied, e.g. if 1/f noise plays a dominant role, the decay of the sequence fidelity can have non-exponential components [105][106][107] and a single exponential fit will not produce the true average gate fidelity [108,109]. The filter function formalism suggests itself to numerically probe RB experiments in such systems for two reasons. First, it enables the study of gate performance subject to noise with correlation times longer than individual gate times. This regime, where a simple description in terms of individual, isolated quantum operations fails, is accessible in the filter function formalism because universal classical noise can be included by the power spectral density S(ω). Second, the simulation of a RB experiment can be performed efficiently by using the concatenation property. Because RB sequences are compiled from a limited set of gates whose filter functions may be precomputed, one only needs to concatenate m filter functions for a single sequence of length m to gain access to the survival probability.
Since for sufficiently long RB sequences r ∈ O(1), and we would need to include the frequency shifts ∆ in a full simulation following Eq. (9) because the low-noise approximation Eq. (42) does not hold in this regime. Unfortunately, the concatenation property does not hold for ∆. Therefore, we focus on the high-fidelity regime where the exponential decay of the sequence fidelity may be approximated to linear order and only the decay amplitudes Γ need to be considered.
In order to evaluate the survival probability of a RB experiment using filter functions, we employ the state fidelity from Section II D 2 and focus on the single-qubit case with d = 2 and the (normalized) Pauli basis from Eq. (53). Because the ideal action of a RB sequence is the identity we have Q = 1. Assuming we prepare and measure in the computational basis, |ψ ∈ {|0 , |1 } so that √ 2 |ρ = |σ 0 ± |σ 3 , we simplify Eq. (49) to For the second equality we used thatŨ is trace-preserving and unital (c.f. Section II A 2) while in the last step we approximated the expression using Eqs. (22) and (42). For our simulation, we neglect SPAM errors so that A = B = 0.5, choose |ψ = |0 , and approximate Eq. (57) as for small gate errors r 1 since this is the regime which we can efficiently simulate using the concatenation property.
We simulate single-qubit SRB experiments using three different gate sets to generate the 24 elements of the Clifford group. For the first gate set we implement the group by naive "single" rotations about the symmetry axes of the cube. Each pulse corresponds to a single time segment during which one rotation is performed so that the j-th element is given by Q j = exp(−iφ j n j · σ). We compile the other two gate sets from primitive π/2 x-and y-rotations so that on average each Clifford gate consists of 3.75 primitive gates (see Ref. 110). For the specific implementation of the primitive π/2 -gates we compare "naive" rotations, i.e. with a single time segment so that Q j = exp( −iπσ j /4 ) for j ∈ {x, y}, and the "optimized" gates from Ref. 94. Pulse durations are chosen such that the average duration of all 24 Clifford gates generated from a single gate set is equal for all three gate sets. This is to ensure that the different implementations of the Clifford gates are sensitive to the same noise frequencies.
We investigate white noise and correlated noise with S(ω) ∝ ω −0.7 assuming the same noise spectrum on each Cartesian axis of the Bloch sphere and normalize the noise power for each gate set and noise type (white and correlated) so that the average Clifford infidelity r is the same throughout. We then randomly draw K = 100 sequences for 11 different lengths m ∈ [1, 101] and concatenate the m Clifford gates using Eq. (29) to compute the control matrix of the entire sequence. For the integral in Eq. (24) we choose the ultraviolet cutoff frequency two orders of magnitude above the inverse duration of the shortest pulse, f UV = 10 2 τ min . Similarly, the infrared cutoff is chosen as f IR = 10 −2 m max τ max with m max = 101 and τ max = 7τ min (since the longest gate is compiled from seven primitive gates with duration τ min ) to guarantee that all nontrivial structure of the filter functions is resolved at small frequencies [111]. Finally, we fit Eq. (59) to the infidelities computed for the different noise spectra.
The results of the simulation are shown in Fig. 5 (a) and (b) for white and correlated noise, respectively. For white noise, the survival probability agrees well with the SRB prediction for all gate types whereas for 1/f -like noise the "single" gates (green pluses) deviate considerably. Hence, fitting the zeroth-order SRB model to such data will not reveal the true average gate fidelity although errors are of order unity. We note that Refs. 105 and 37 found similar results using different methods for 1/f and perfectly correlated DC noise, respectively. The former observed SRB to estimate r within 25 % and the latter found the mean of the SRB fidelity distribution to deviate from the mode, thereby giving rise to incorrectly estimated fidelities.
On top of affirming the findings by the references, our results demonstrate that the accuracy of the predictions made by SRB theory, i.e. that the RB decay rate directly corresponds to the average error rate of the gates, not only depends on the gate implementation but also on which error channels are assumed. This can be seen from the inset of Fig. 5(b), where only dephasing noise (σ z ) contributions are shown. For this noise channel and the "naive" gates, one finds a slower RB decay than expected from the actual average gate fidelity, so that the latter would be overestimated by an RB experiment, whereas the "single" gates show the opposite behavior. Depending on the gate set and relevant error channels, non-Markovian noise may thus even lead to improved sequence fidelities due to errors interfering destructively. This behavior is captured by the pulse correlation filter functions whose contributions to the sequence fidelity lead to the deviations from the SRB prediction.
Notably, the data for the "optimized" gates agree with the prediction for every noise channel individually which implies that correlations between pulses are suppressed. This highlights the formalism's attractiveness for numerical gate optimization as the pulse correlation filter functions F (gg ) (ω) may be exploited to suppress correlation errors. To be more explicit, the correlation decay amplitudes Γ (gg ) from Eq. (31) can be used to construct cost functions for quantum optimal control algorithms like GRAPE [112,113] or CRAB [114]. By constructing linear combinations of Γ (gg ) with different pulse indices g and g , correlations between any number of pulses can be specifically targeted and suppressed using numerical pulse optimization.

D. Quantum Fourier transform
To demonstrate the flexibility of our software implementation, we calculate filter functions for a fourqubit quantum Fourier transform (QFT) [45,115] circuit. QFT plays an important role in many quantum algorithms such as Shor's algorithm [116] and quantum  59) to the data while the solid black lines correspond to a zeroth-order SRB model with A = B = 0.5 and the true average gate infidelity per Clifford r. Errorbars show the standard deviation of the SRB sequence fidelities, illustrating that for the "single" gate set noise correlations can lead to amplified destructive and constructive interference of errors. The same noise spectrum is used for all three error channels (σx, σy, σz) and the large plots show the sum of all contributions. (a) Uncorrelated white noise with the noise power adjusted for each gate type so that the average error per gate r is constant over all gate types. No notable deviation is seen between different gate types. (b) Correlated 1/f -like noise with noise power adjusted to match the average Clifford fidelity in (a). The decay of the "single" gateset differs considerably from that of the other gate sets and the SRB decay expected for the given average gate fidelity, whereas "naive" and "optimized" gates match the zeroth order SRB model well, indicating that correlations in the noise affect the relation between SRB decay and average gate fidelity in a gateset-dependent way. Inset: contributions from σz-noise show that the sequence fidelity can be better than expected for certain gate types and noise channels.
phase estimation [45]. For the underlying gate set, we assume a standard Rabi driving model with IQ control and nearest neighbor exchange. That is, we assume full control of the x-and y-axes of the individual qubits as well as the exchange interaction mediating coupling between two neighboring qubits. This system allows for native access to the minimal gateset G = {X i ( π/2 ), Y i ( π/2 ), CR ij ( π 2 3 )} where CR ij (φ) denotes a controlled rotation by φ about z with control qubit i and target qubit j. Controlled-z rotations by angles π/2 m as required for the QFT can thus be obtained by concatenating 2 3−m minimal gates CR ij ( π 2 3 ).
Despite native access to all necessary gates, we employ QuTiP's implementation [42] of the GRAPE algorithm [112,113] to generate the gates in order to highlight our method's suitability for numerically optimized pulses. For the optimization we choose a time step of ∆t = 1 ns and a total gate duration of τ = 30 ns. For completeness, see Appendix C for details on the optimized gates. We then construct the remaining required gates by sequencing these elementary gates, i.e. the Hadamard gate where B • A denotes the composition of gates A and B such that gate A is executed before gate B. To map the canonical circuit [45] onto our specific qubit layout with only nearest-neighbor coupling, we furthermore introduce SWAP operations to couple distant qubits. These gates can be implemented by three CNOTs, SWAP ij = CNOT ij •CNOT ji •CNOT ij . The CNOTs in turn are obtained by a Hadamard transform of the controlled phase gate, CNOT ij = H j • CR ij (π) • H j . The complete quantum circuit is shown at the top of Fig. 6; for the canonical circuit with all-to-all connectivity refer to Ref. 45. In total, there are 442 elementary pulses, 198 of which are required for the three SWAPs on the first two qubits, so that the entire algorithm would take ∼ 13 µs to run. Note that the circuit could be compressed in time by parallelizing some operations but for simplicity we only execute gates sequentially and do not execute dedicated idling gates.
In order to leverage the extensibility of the filter function approach (see Section III B), we use a Pauli basis for the pulses and proceed as follows: 1. Instantiate the PulseSequence objects for the elementary gates G for the first two qubits and cache the control matrices.
2. Compile all required single-and two-qubit pulses by concatenating the PulseSequences that implement G.
3. Extend the PulseSequences to the full four-qubit Hilbert space.
For our gate parameters and 400 frequency points, this procedure takes around 5 s on an Intel ® Core™ i9-9900K eight-core processor, whereas computing the filter functions naively using Eq. (35) takes around 4 min. The resulting filter functions are shown in Fig. 6 for the noise operators affecting the first qubit; for an in-depth discussion and validation of the fidelities predicted, see the accompanying letter Ref. 16 and its supplementary information. Evidently, the fidelity of the algorithm is most  26) and (48)) and suggest that high frequencies up to the knee at ω ≈ 10 ns −1 cannot be neglected if the cutoff frequency of the noise is sufficiently high or the spectrum does not drop off quickly enough (note the linear scale as opposed to the logarithmic scale for the filter functions). susceptible to DC noise; below roughly ω 10 −3 ns −1 the filter functions level off at their maximum value. In the GHz range there is a plateau with sharp peaks corresponding to the n-th harmonics of the inverse pulse duration ω n = 2πn/τ , where the leftmost belongs to n = 1. The dashed lines show the error sensitivities I α (ω 1 , ω 2 ) := ω2 ω1 dω F (ω) in the frequency band [0, ω] relative to the total sensitivity I α (0, ∞). For a white spectrum, i.e. S(ω) = const., this quantifies the fraction of the total entanglement infidelity that is accumulated up to frequency ω (c.f. Eqs. (26) and (48)). Thus, to obtain a precise estimate of the algorithm's fidelity, five frequency decades need to be taken into account.
These insights demonstrate that our method represents a useful tool to analyze how and to which degree small algorithms are affected by correlated errors, and how this effect depends on the gate implementation. It could thus also be used to choose or optimize gates in an algorithmspecific way.

VI. FURTHER CONSIDERATIONS
Before we conclude, let us address two possible avenues for future work, one for the formalism itself and one for its application.
To extend our approach to the filter function formalism beyond the scope discussed in this work, the most evident path forward is to allow for quantum mechanical baths instead of purely classical ones. Such an extension would facilitate studying for example nonunital T 1 -like processes. In fact, the filter function formalism was originally introduced considering quantum baths such as spin-boson models [18,19] or more general baths [17,23,117], but it remains an open question whether this can be applied to our presentation of the formalism and the numerical implementation in particular. In a fully quantum-mechanical treatment, (sufficiently weak) noise coupling into the quantum system can be modelled via a set of bath operators {D α (t)} α so that H n (t) = α B α (t) ⊗ D α (t) (the classical case is recovered by replacing D α (t) → b α (t)1) [118]. Accordingly, the ensemble average over the stochastic bath variables {b α (t)} α needs to be replaced by the quantum expectation value tr B ( · ρ B ) with respect to the state ρ B of the bath B. One therefore needs to deal with correlation functions of bath operators instead of stochastic variables. An immediate consequence for numerical applications is hence an increased dimensionality of the system, which could be dealt with by using analytical expressions for the partial trace over the bath.
For future applications of our method, it would be interesting to study the effects of noise correlations in quantum error correction (QEC) schemes [119][120][121]. While extensive research has been performed on QEC, noise is usually assumed to be uncorrelated between error correction cycles. In this respect, our formalism may shed light on effects that need to be taken into account for a realistic description of the protocol. As outlined above, we can compute expectation values of (stabilizer) measurements in a straightforward manner from the error transfer matrix. Unfortunately, this implies performing the ensemble average over different noise realizations, therefore removing all correlations between subsequent measurement outcomes for a given noise realization. Hence, the same feature that allows us to calculate the quantum process for correlated noise, namely that we compute only the final map by averaging over all "paths" leading to it, prevents us from studying correlations between consecutive cycles. To overcome this limitation in the context of quantum memory one could invoke the principle of deferred measurement [45] and move all measurements to the end of the circuit, replacing classically controlled operations dependent on the measurement outcomes by conditional quantum operations. Alternatively, to incorporate the probabilistic nature of measurements, one could devise a branching model that implements the classically controlled recovery operation by following both conditional branches of measurement outcomes with weights corresponding to the measurement probabilities as computed from the ensemble-averaged error transfer matrix. An intriguing connection also exists to the quantum Zeno effect, for which quantum systems subject to periodic projective measurements have been identified with a filter function [17,122,123].

VII. CONCLUSION AND OUTLOOK
As quantum control schemes become more sophisticated and take into account realistic hardware constraints and sequencing effects, their analytic description becomes cumbersome, making numerical tools invaluable for analyzing pulse performance. In the above, we have shown that the filter function formalism lends itself naturally to these tasks since the central objects of our formulation, the interaction picture noise operators, obey a simple composition rule which can be utilized to efficiently calculate them for a sequence of quantum gates. Because the nature of the noise is encoded in a power spectral density in the frequency domain, its effects are isolated from the description of the control until they are evaluated by the overlap integral of noise spectrum and filter function. Hence, the noise operators are highly reusable in calculations and can serve as an economic way of simulating pulse sequences.
Building on the results of a separate publication [16], we have presented a general framework to study decoherence mechanisms and pulse correlations in quantum systems coupled to generic classical noise environments. By combining the quantum operations and filter function formalisms, we have shown how to compute the Liouville representation of the exact error channel of an arbitrary control operation in the presence of Gaussian noise. For non-Gaussian noise our results become perturbative in the noise strength. Furthermore, we have introduced the filter_functions Python software package that implements the aforementioned method. We showed both analytically and numerically that our software implementation can outperform Monte Carlo techniques by orders of magnitude. By employing the formalism and software to study several examples we demonstrated the wide range of possible applications.
The capacity for applications in quantum optimal control has already been established above. In a forthcoming publication, we will present analytical derivatives for the fidelity filter function, Eq. (26), and their implementation in the software package [90]. Together with the infidelity, Eq. (48), they can serve as efficient cost functions for pulse optimization in the presence of realistic, correlated noise [43,87]. Since our method offers insight into correlations between pulses at different positions in a sequence, the pulse correlation filter function F (gg ) (ω) with g = g can additionally serve as a tool for studying under which conditions pulses decouple from noise with long correlation times. Such insight would be valuable to design pulses for algorithms. Another interesting application could be quantum error correction in the regime of long-time correlated noise as outlined above in Section VI, where we also briefly touched upon a possible extension of the framework to quantum mechanical baths.
The tools presented here, both analytical and numerical as implemented in the filter_functions software package [41], provide an accessible way for computing filter functions in generic control settings across the different material platforms employed in quantum technologies and beyond. The inner integration is simple to perform and we get if ω + Ω (g) Shifting the limits of integration and performing integration by parts in the case ω + Ω In the main text, we claimed that the contribution of noise sources (α, β) to the total entanglement infidelity I e (Ũ) = αβ I αβ reduces from the trace of the cumulant function K to To show this, we substitute K by its definition in terms of ∆ and Γ according to Eq. (17). This yields for the trace since ∆ is antisymmetric. In order to further simplify the trace tensors on the right hand side of Eq. (A13), we observe that the orthonormality and completeness of the operator basis C defining the Liouville representation of K (c.f. Eq. (14)) is equivalent to requiring that C † C = 1 with C reshaped into a d 2 × d 2 matrix by a suitable mapping. This condition may also be written as because every element C k is Hermitian. Using this relation in Eq. (A13) then yields tr K αβ = − 1 2 kl Γ αβ,kl (2dδ kl − 2 tr(C k ) tr(C l )) = −d tr Γ αβ . (A15) The last equality only holds true for bases with a single non-traceless element (the identity), such as the bases discussed in Section III C. This is because in this case, tr(C k ) = 0 for k > 0 whereas Γ αβ,kl = 0 for either k = 0 or l = 0 since Γ is a function of the traceless noise Hamiltonian for which tr(C 0 H n ) ∝ tr H n = 0 (i.e. the first column of the control matrix is zero, see Eqs. (16) and (21)). Finally, substituting Eq. (A15) into Eq. (A11) we obtain our result In this appendix we lay out in more detail how the fidelity of the optimized S-T 0 qubit gates from Ref. 94 was calculated using filter functions. In two singlet-triplet qubits, angular momentum conservation suppresses occupancy of states with non-vanishing magnetic spin quantum number m s so that the total accessible state space of dimension d = 6 is spanned by {|↑↓↑↓ , |↑↓↓↑ , |↓↑↑↓ , |↓↑↓↑ , |↑↑↓↓ , |↓↓↑↑ }. A straightforward method to single out the computational subspace (CS) dynamics from those on the whole space would be to simply project the error transfer matrixŨ ≈ 1+K with K the cumulant function onto the CS as proposed by Wood and Gambetta [76], that is calculate the fidelity as F e = tr Π cŨ /d 2 c where Π c is the Liouville representation of the projector onto the CS and d c = 4 the dimension of the CS. However, here we use a more involved procedure in order to gain more insight from the error transfer matrix as well as to obtain a better comparison to the fidelities computed by Cerfontaine et al. [94], who map the final 6 × 6 propagator to the closest unitary on the 4 × 4 CS during their Monte Carlo simulation.
To calculate the fidelity of the target unitary on the 4 × 4 CS, we thus construct an orthonormal operator basis C of the full 6 × 6 space that is partitioned into elements which are nontrivial only on the CS on the one hand and elements which are nontrivial only on the remaining space on the other such that C = C c ∪ C . Using such a basis, we can then trace only over CS elements of the error transfer matrixŨ in Eq. (45) to obtain the fidelity of the gate on the CS. Moreover, we retain the opportunity to characterize the gates on the basis of the Pauli matrices.
Since there is no obvious way to extend the Pauli basis for two qubits to the complete space we proceed as follows: For the CS, we pad the two-qubit Pauli basis with zeros on the leakage levels, i.e.
since for unitary operations on the CS we have K 00 ≈ 1 −Ũ 00 = 1 − tr C c 0Ũ C c 0Ũ † = 0. Hence, excludingŨ 00 from the trace corresponds to partially disregarding non-unitary components of the error channel on the computational subspace. Although not the only element that differs compared to the closest subspace unitary, K 00 contains the most obvious contribution, whereas those of other elements are more difficult to disentangle into unitary and non-unitary components.
Similar to the fidelity, we also obtain the canonical filter function shown in panel (b) of Fig. 3 by summing only over columns one through 15 of the control matrix, F ij (ω) = 15 k=1 B ij k (ω) 2 . In fact, including the first column, corresponding to the padded identity matrix C c 0 , in the filter function removes the DCG character of F 12 (ω) and F 34 (ω), which instead approach a constant level of around 20 (note that the filter function is dimensionless in our units) at zero frequency. This is consistent with the fact that the gates were optimized using quasistatic and fast white noise contributions to the fidelity after mapping to the closest unitary on the computational subspace. We have performed Monte Carlo resimulations that support this reading. In Fig. 7 we show the filter functions once including and once excluding the contributions from C c 0 .  1, 1, 1, 0, 0) for the computational subspace. Evidently, including C c 0 removes the DCG character, namely that F ij (ω) → 0 as ω → 0, of the gates but has little effect on the high-frequency behavior. As the pulse optimization minimizes, among other figures of merit, the infidelity of the final propagator mapped to the closest unitary on the computational subspace due to quasistatic and fast white noise, this indicates that excluding C c 0 from the filter function corresponds to partially neglecting non-unitary components of the propagator on the computational subspace. . Note that the optimization is neither very sophisticated nor realistic as the algorithm only maximizes the systematic (coherent) fidelity tr U Q † targ /d and the randomly distributed initial control amplitudes are not subject to any constraints.

Cauchy-Schwarz inequality we then have
is the maximum value that the noise assumes during the pulse, ϑ (g) (t) = Θ(t − t g−1 ) − Θ(t − t g ) is one during the g-th time interval and zero else, and where we approximated the time evolution as piecewise constant.
Then, in order to guarantee convergence of the ME, in terms of the root mean square value δb α . That is, b (m) α = C m b α (0) 2 1/2 = C m δb α for a sufficiently large value C m . Finally, realizing that δb 2 α = dω 2π S α (ω) and by the triangle inequality, where we have introduced the parameter ξ. Thus, the expansion converges if ξ < π/C m . However, we note that in practice the rms noise amplitude δb α will often be infinite, limiting the usefulness of this bound for certain noise spectra.

Infidelity
Again assuming a time dependence B α (t) = s α (t)B α as well as piecewise constant control, we note that for the infidelity we have (c.f. Eq. (45)) where, going from the second to the third line, we have factored out the total power of noise source α from the cross-correlation function, b α (t 1 )b α (t 2 ) = b 2 α (0) b α (t 1 )b α (t 2 ) . Thus, the first order infidelity Eq. (45) is upper bounded by ξ 2 d , the same parameter also bounding the convergence of the ME, and higher orders can be neglected if ξ 2 1. Note that similar arguments can be made for the higher orders of the ME [29]. In particular, the n-th order ME term containing n-point correlation functions of the noise is of order O(ξ n ) as stated in the main text.