{'fraction_non_alphanumeric': 0.04163016919459684, 'fraction_numerical': 0.011064496795906136, 'mean_word_length': 4.663707571801567, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient classical simulation of open bosonic quantum systems

We propose a computationally efficient method to solve the dynamics of operators of bosonic quantum systems coupled to their environments. The method maps the operator under interest to a set of complex-valued functions, and its adjoint master equation to a set of partial differential equations for these functions, which are subsequently solved numerically. In the limit of weak coupling to the environment, the mapping of the operator enables storing the operator efficiently during the simulation, leading to approximately quadratic improvement in the memory consumption compared with the direct approach of solving the adjoint master equation in the number basis, while retaining the computation time comparable. Moreover, the method enables efficient parallelization which allows to optimize for the actual computational time to reach an approximately quadratic speed up, while retaining the memory consumption comparable to the direct approach. We foresee the method to prove useful, e.g., for the verification of the operation of superconducting quantum processors.

We propose a computationally efficient method to solve the dynamics of operators of bosonic quantum systems coupled to their environments. The method maps the operator under interest to a set of complex-valued functions, and its adjoint master equation to a set of partial differential equations for these functions, which are subsequently solved numerically. In the limit of weak coupling to the environment, the mapping of the operator enables storing the operator efficiently during the simulation, leading to approximately quadratic improvement in the memory consumption compared with the direct approach of solving the adjoint master equation in the number basis, while retaining the computation time comparable. Moreover, the method enables efficient parallelization which allows to optimize for the actual computational time to reach an approximately quadratic speed up, while retaining the memory consumption comparable to the direct approach. We foresee the method to prove useful, e.g., for the verification of the operation of superconducting quantum processors.
Introduction.-In recent years, the potential to revolutionize a spectrum of industries [1], such as finance [2], pharmacy [3], and cybersecurity [4], has given a major impetus towards realizations of practical-use quantum computers [5][6][7][8][9][10][11][12]. The role of classical simulations of such devices in the progress is of great importance. First, accurate classical simulations contribute to the optimization of the control of a quantum processor to maximize the fidelity of the target quantum algorithm. Second, while there has already been claims of achieving the quantum supremacy [13,14], these experiments need to be verified against fully controllable classical simulations. Third, a fair proof of quantum supremacy requires a comparable effort put on the development of classical simulations of the given experimental setup, which may raise the bar for quantum computers.
Currently, some of the most prominent architectures for programmable quantum computers are based on superconducting technology [13,15]. Such devices, among a range of other quantum technological applications [16][17][18], such as those employing quantum optics [14,19] and circuit quantum electrodynamics [20,21], can predominantly be considered as bosonic systems. Recent progress in the development of quantum technology, especially in superconducting circuits [22][23][24], has enabled reaching the limit of very weak interaction between the system under interest and its environment [15,25]. Therefore, the weak interaction limit has become increasingly rele-vant also for classical simulations.
Typically, the experiments with quantum processors terminate with the measurement of the qubits in their computational basis. However, a variety of problems of experimental interest use instead measurements of certain other observables, including certain non-commuting observables [26], correlation functions [27], and dynamical linear response functions [28]. The Heisenberg picture is the natural picture for the implementation of classical simulations of such schemes. The direct approach to solve the dynamics of a given operator in the Heisenberg picture is to transform its adjoint master equation into a set of coupled scalar differential equations corresponding to the number basis, or to some other basis [29][30][31], and to integrate this set numerically. In this Letter, we refer to the former of these methods to as the conventional approach. Alternative techniques include the socalled coherent state representation [32] and the phase space formulation [33,34]. However, in general these approaches are memory-consuming which poses stringent requirements for the utilized computational hardware.
In this Letter, we push the boundaries of classical simulations of weakly dissipative bosonic quantum systems in the Heisenberg picture. The proposed approach is based on expanding the operator under interest in a problem-specific operator basis, mapping the adjoint master equation of the operator into a complex-valued partial differential equation for the generating function arXiv:2203.05526v2 [quant-ph] 11 Mar 2022 Schematic process chart of the proposed method to solve the dynamics of operators of open bosonic quantum systems. The symbols are defined in the main text in the vicinity of the related equations. For clarity, we have used the short-hand notation O for the temporal evolution operator of the generating function G and similar notations for those of Gx, Gy and R, and suppressed the temporal dependencies of the Hamiltonian and the temporal evolution operators. of this expansion, performing a simple and well-defined transformation to this dynamic equation, and solving the transformed dynamic equation numerically. In the process, the operator is mapped to a set of complex-valued functions in a fashion that enables storing it memoryefficiently during the simulation. The approach poses no assumptions on the structure or on the initial conditions of the system.
To demonstrate the power of the proposed method, we consider a network of coupled driven dissipative anharmonic oscillators, routinely used as a model of superconducting quantum processors based, for example, on transmon qubits [35,36]. To this end, we first use a well-known property of Markovian master equations to show that the most commonly used master equation for the anharmonic oscillator breaks down already for arbitrarily weak anharmonicity, and consequently derive the proper master equation starting from the first principles. Subsequently, we use the conventional approach and our approach to numerically solve the proper adjoint master equation of the system and compare the memory consumptions and computation times of the methods.
The method.-In this section, we provide a general description of the proposed method, see Fig. 1. For the sake of clarity, we limit the consideration to single-mode systems. However, the method straightforwardly generalizes to an arbitrary number of discrete bosonic modes as we show in [37]. Moreover, here we do not construct the dynamics for the general dissipators. Such a treatment is given subsequently when we consider the example system of anharmonic oscillators.
We begin by introducing two elementary operators of the system under interest,ô 1 andô 2 , such that (i) they obey the commutation relation [ô 2 ,ô 1 ] = rÎ, where r ∈ C \ {0} andÎ is the identity operator, and (ii) their monomials, {ô k 1ô l 2 } ∞ k,l=0 , form a complete operator basis. Consequently, we expand the Hamiltonian of the system in the Heisenberg picture in this basis aŝ and an arbitrary system operator aŝ Examples of elementary operators that satisfy the above conditions (i) and (ii) are the creation and annihilation operators of a bosonic mode,ô 1 =â † andô 2 =â, for which r = 1, and the corresponding quadratures,ô 1 =p andô 2 =q, for which r = i. The utility in the introduction of these generic elementary operators is clarified below. Let us consider a generic adjoint master equation of the system [38], where the dot denotes the temporal derivative,L n is the Lindblad jump operator of the noise channel n, and {Ô 1 ,Ô 2 } =Ô 1Ô2 +Ô 2Ô1 is the anticommutator. Operator differential equations are in general difficult to solve directly. To trasform Eq. (3) into a c-number differential equation, we insert the expansions of Eqs. (1) and (2) into Eq. (3) and define a complex-valued function F (x, y; t) = ∞ k,l=0 f k,l (t)x k y l that generates the expansion coefficients of operatorÂ(t). A lengthy calculation yields the dynamic equation for the generating function [37], where : * : denotes ordering of the variables to the left with respect to the differentiation operators, H(x, y; t) = ∞ u,v=0 h u,v (t)x u y v , and O n is a multivariate function of the coordinates and the differential operators associated with the noise channel n. For brevity, we omit the general form of O n . However, in Eq. (S21) in [37] we implicitly provide it for the anharmonic oscillator, and the technique used therein to derive it is also applicable for other systems of interest. Note that the choice of the function basis in the definition of F is arbitrary and different choices lead to different forms of the dynamic equation.
To the best of our knowledge, the adjoint master equation has not previously been mapped to the form of Eq. (4). Note that Eq. (4) can be readily solved numerically. However, we take here a more involved approach by defining a new generating function as G(x, y; t) = e r −1 xy F (x, y; t). The idea is to transform the dynamic equation into an almost separable form to enable storing the generating function in terms of univariate functions and a residual function. When discretizing these functions for simulation, the univariate functions and the residual function can be stored more compactly than the untransformed bivariate generating function. The transformed dynamic equation readṡ where O n = e r −1 xy O n e −r −1 xy . The new generating function G(x, y; t) = ∞ k,l=0 g k,l (t)x k y l can be expressed in terms of the original coefficients as The new generating function corresponds to the operator expansionÂ where the operator basis elements are given bŷ which can be verified by inserting Eqs. (6) and (8) into Eq. (7). Next, we introduce three new complex-valued functions according to G(x, y; t) = G x (x; t)G y (y; t) + R(x, y; t), which are governed bẏ where C is the separation constant. Suppose that the new generating function is initially separable, that is, G(x, y; 0) = G x (x; 0)G y (y; 0) and R(x, y; 0) = 0. From Eq. (9c) it is evident that if κ n t is small for all n, then the function R can be approximated as a zero function. Moreover, in the following we find that in the case of a network of damped anharmonic oscillators and larger κ n t, R can be approximated as a single-variable function. This implies that the transformation enables storing the operator memory-efficiently as a set of three single-variable functions intead of one two-variable function. Note that any G(x, y; 0) can be expressed as a linear combination of separable functions, e.g. as Thus, by solving the linear dynamics of each of such initially separable function G (p,q) (x, y; 0), we can reconstruct the dynamics of G(x, y; t). The initial expansion coefficients read where g (p,q) k,l (t) are the expansion coefficients of G (p,q) (x, y; t) in the chosen function basis {x k y l } ∞ k,l=0 . Note that this approach requires solving the dynamics of N 2 thr generating functions, where N thr denotes the truncation of the Taylor expansion in Eq. (10). This number is reduced to N thr if the initial operator is sub-, super-, or diagonal in the chosen operator basis, that is, where d ∈ Z defines whether the initial operator is sub-(d < 0), super-(d > 0), or diagonal (d = 0). In this case, the initial expansion coefficients read Examples of operators that possess this favorable property include the number state projectors, the powers and cumulants ofô 1ô2 , the operator defining the leakage error from a subspace spanned by number states, and the squeezing operator. Note that the dynamics of each G (p,q) (x, y; t) is independent. Thus, the numerical integration of these transformed generating functions can be parallelized in an ideal fashion.
In summary, the method presented in this section can be used to solve the dynamics of an operator of an open bosonic system as follows: (i) Determine the initial transformed expansion coefficients g p,q (0) using Eqs. (2) and (6). (ii) Numerically solve Eqs. (9a)-(9c) with the initial condition G (p,q) (x, y; 0) = x p y q for each nonnegligible g p,q (0). (iii) If desired, reconstruct the original expansion coefficients of the operator using Eq. (11).
Master equation for the anharmonic oscillator.-In this section, we present a master equation for the anharmonic oscillator, which we subsequently solve to demonstrate the power of our method. To justify solving it instead of the conventionally used master equation, we first argue that the latter is not accurate even for low values of anharmonicity.
It has been shown that the von Neumann entropy production rate Π =Ṡ − T −1 EĖ of any Markovian quantum system is non-negative [38], where S = −k B Tr(ρ lnρ) is the von Neumann entropy of the system, T E is the temperature of its environment, E = Ĥ is the mean energy of the system, k B is the Boltzmann constant, andρ is the density operator of the system. Consider the commonly used master equation of the anharmonic oscillator [39][40][41][42][43] where ω is the angular frequency difference between the two least energetic states, η is the anharmonicity, and n = [e ω/(kBTE) − 1] −1 is the mean number of excitations in the environment in the mode with the angular frequency ω. In [37] we show that given a diagonal initial state of the formρ(0) = N ∞ n=0 e −cn |n n|, where N is a normalization constant, for all non-zero values of η there exists c > 0 such that the conventionally used master equation leads to a negative entropy production rate. Consequently, the conventional master equation is not accurate even in the limit of weak anharmonicity.
Due to this observation, in [37] we derive the proper Markovian master equation starting from the first principles, following Refs. [44] and [45]. The master equation readṡ A similar but less explicit form of the master equation of the anharmonic oscillator is given in [46]. In [37], we show analytically that this master equation cures the inconsistency presented above in the limit of small T E .
Example: network of coupled driven dissipative anharmonic oscillators.-Let us consider a network of N bilinearly coupled, classically driven anharmonic oscillators. Adopting the rotating-wave approximation, the Hamil-tonian of the system in the Schrödinger picture readŝ where ω d,j is the drive frequency and Ω j is the complex drive amplitude of the anharmonic oscillator j, H.c. denotes the Hermitian conjugate of the previous term, and J j,j = J j ,j is the coupling constant between the oscillators j and j .
In the numerical simulations below, we approximate the Hamiltonian to be constant during each small time step. Thus, we can simply replace the Hamiltonian of the system transformed into the Heisenberg picture in the adjoint master equation (3) by the above Schrödinger picture Hamiltonian. Moreover, we model the dissipation in each oscillator with the dissipator presented in Eq. (15), and for simplicity, assume T E = 0 and choose the elementary operators asô 1 =â † and o 2 =â. Finally, we define G(x 1 , y 1 , . . . , x N , y N ; t) = F (x 1 , y 1 , . . . , x N , y N ; t) N j=1 e xj yj to obtain the following dynamic equation for the transformed generating functionĠ where where we have used the short-hand notation G ≡ G(x 1 , y 1 , . . . , x N , y N ; t). Similarly to the previous sections, we introduce three new complex-valued functions, which yieldṡ To solve the dynamics of the operator under interest, Eqs. (19a)-(19c) are solved with several initial conditions, G (p1,q1,...,p N ,q N ) (x 1 , y 1 , . . . , x N , y N ; 0) = N j=1 x pj j y qj j . The separation constant C is chosen such that it provides sufficient numerical stability. Here, we use the constant that minimizes the energy of the normalized derivatives, |Ġ x /G x | 2 + |Ġ y /G y | 2 , of the uncoupled system at t = 0. Such a constant depends on the initial conditions and reads We implement numerical solvers based on our approach and on the conventional approach. The physical simulation parameters are chosen to correspond to a typical superconducting quantum processor, and are given in Table I. The operator for which we solve the dynamics isÂ(0) = N j=1â j . For both methods, the numerical integration is implemented using the fourth-order Runge-Kutta method. By tuning the threshold for truncating the representation of the operator, we force the relative errors r to be approximately constant, 100 ppm, for all the simulations. The relative error is defined as where * 2 is the elementwise L 2 norm, M is the representation of the operator in the number basis truncated to three least energetic states, and M ref is that given by the reference method. The reference method is the conventional method executed by storing a larger representation during the simulation, and by using half the time step.
The results in Fig. 2 show that both the number of the stored elements and the number of operations increase with increasing time of the evolution of the system for both of the methods. Furthermore, we find that increasing κ j t to values for which the residual function R cannot be approximated as a zero function, it is nevertheless well-approximated by a function of N variables only, namely, {x j y j } N j=1 . The number of stored elements of the conventional approach scales as the square of that required by our approach as a function of the number of subsystems. This is expected since our approach only requires storing N -variable functions during the simulation whereas the conventional approach requires storing a 2N -variable function. Finally, the number of operations scales similarly for both of the methods as a function of N . This is expected since increasing N increases the number of generating functions, the dynamics of which are to be solved for our method.
Conclusions.-We have introduced a numerically efficient approach to solve the dynamics of certain operators of weakly open bosonic systems. The approach is based on expanding the operators of the system in a problemspecific operator order, mapping the adjoint master equation into a dynamic equation for a generating function of this expansion, applying a transformation to the dynamic equation, and solving it by numerical integration. By considering a network of classically driven damped anharmonic oscillators, we demonstrate that our approach reduces the memory consumption compared to the conventional method quadratically. Moreover, we observe that the computation times of the methods are comparable. However, parallelization of the computationally heavy steps of our method enables reaching a speed-up compared with the conventional approach. In the future, the generality of the proposed approach enables applying it to a range of weakly open bosonic systems under topical interest. Moreover, the technique may be further generalized by considering expansions in more general operator orders and function bases. This may broaden the applicability of the proposed technique even further.
In this section, we derive the dynamic equation for the generating function of an operator of a general bosonic system consisting of N discrete modes. For simplicity, we assume here that the system is closed. However, we demonstrate a derivation of the part of the dynamic equation related to the Lindblad part of the adjoint master equation for the anharmonic oscillator in Eqs. (S19)-(S21). Similarly to the main text, we define for each mode j two elementary operatorsô j,1 andô j,2 that obey the commutation relation [ô j,2 ,ô j,1 ] = r jÎ , and further note that the operators between different modes commute. In terms of these operators, the Hamiltonian of the system in the Heisenberg picture can be expanded asĤ where u = {u 1 , . . . , u N } and v = {v 1 , . . . , v N }, and an arbitrary operator of the system aŝ In the following, the summations range from 0 to ∞ for each index unless otherwise stated. Assuming that the operator does not depend explicitly on time, it evolves according toȦ(t) = i Ĥ (t),Â(t) , which leads to where in the second step we have used the commutation relations of operator monomials [1]. By defining the binomial coefficients as n m = Γ(n+1) Γ(m+1)Γ(n−m+1) we can extend the summations over s j 's to infinity in order to simplify the reorganization of the summations. Consequently, Using the completeness of the operator monomials we directly obtaiṅ Multiplying the both sides by Defining a generating function related to the operatorÂ(t) in the chosen operator order as F (x, y; t) = k,l f k,l (t) N j=1 x kj j y lj j finally gives uṡ where r = (r 1 , . . . , r N ) T , T denotes the transpose, : * : is an operator that orders the variables to the left with respect to the differentiation operators, • denotes the Hadamard product, and H(x, This concludes the derivation. Note that for single-mode systems this equation reduces to the first two lines of Eq. (4) of the main text.