⇦Chinese Physics Letters, 2021, Vol. 38, No. 8, Article code 080301 Variational Quantum Algorithms for the Steady States of Open Quantum Systems Huan-Yu Liu (刘环宇)1,2, Tai-Ping Sun (孙太平)1,2, Yu-Chun Wu (吴玉椿)1,2,3*, and Guo-Ping Guo (郭国平)1,2,3,4 Affiliations 1Key Laboratory of Quantum Information, Chinese Academy of Sciences, School of Physics, University of Science and Technology of China, Hefei 230026, China 2CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China 3Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei 230088, China 4Origin Quantum Computing Hefei, Hefei 230026, China Received 9 April 2021; accepted 4 June 2021; published online 2 August 2021 Supported by the National Key Research and Development Program of China (Grant No. 2016YFA0301700), and the Anhui Initiative in Quantum Information Technologies (Grant No. AHY080000).
*Corresponding author. Email: wuyuchun@ustc.edu.cn Citation Text: Liu H Y, Sun T P, Wu Y C, and Guo G P 2021 Chin. Phys. Lett. 38 080301    Abstract The solutions of the problems related to open quantum systems have attracted considerable interest. We propose a variational quantum algorithm to find the steady state of open quantum systems. In this algorithm, we employ parameterized quantum circuits to prepare the purification of the steady state and define the cost function based on the Lindblad master equation, which can be efficiently evaluated with quantum circuits. We then optimize the parameters of the quantum circuit to find the steady state. Numerical simulations are performed on the one-dimensional transverse field Ising model with dissipative channels. The result shows that the fidelity between the optimal mixed state and the true steady state is over 99%. This algorithm is derived from the natural idea of expressing mixed states with purification and it provides a reference for the study of open quantum systems. DOI:10.1088/0256-307X/38/8/080301 © 2021 Chinese Physics Society Article Text The open quantum system includes the interaction between the system and the environment. This interaction causes the state of the system to dissipate besides its evolution under the Hamiltonian of the system, which can be described by the Lindblad master equation.[1] Solving problems related to open quantum systems is of great theoretical and experimental significance,[2–5] the key to which is to solve the steady state density matrix. This can be achieved by continuously iterating the master equation using classical computation. However, the dimension of the Hilbert space grows exponentially with the size of the system, which even makes it difficult for a medium-sized system. The development of quantum computation gives us confidence that quantum computers can surpass their classical counterparts on certain issues. However, due to the limitation of the number of qubits and the fidelity of quantum operations, we are still far away from the age of fault-tolerant quantum computation. In the near future, we would just get access to the Noisy Intermediate-Scale Quantum (NISQ) devices[6–8] with dozens to hundreds of qubits available. Over the years, an NISQ-friendly algorithm,[9,10] the variational quantum algorithm, has received extensive attention. Usually for a task, we employ parameterized quantum circuits (PQC) to prepare ansatz, and we estimate the cost function through quantum measurement with the quantum processor. Subsequently, a classical optimizer is required to update the parameters based on the measured data to minimize the cost function. These steps are performed in a loop between the two parts until the cost function converges, which returns a set of optimal parameters. Successful examples of variational quantum algorithms include variational quantum eigensolver (VQE)[11–15] and variational quantum simulation (VQS),[16–20] whose purposes are to find the ground state of Hamiltonian $H$ and simulate the state evolution of the systems, respectively. In 2019, some works on the time evolution and steady state ansatz construction for open quantum systems based on the Restricted Boltzmann Machine (RBM)[21,22] architecture were published. In these works, RBMs are used to represent the purification of the mixed state or the pure states in the ensemble.[23–25] A different idea is to represent the $n^2$ elements of the $n$-qubit density matrix with a $2n$-qubit pure state:[26] $\rho=\sum_{ij} \rho_{ij} |i\rangle \langle j| \rightarrow |\rho\rangle = \frac 1C \sum_{ij} \rho_{ij} |i\rangle \otimes |j\rangle$. Then, the function of the Lindbladian $\mathcal{L}\rho$ is transformed into a “Hamiltonian”: $\mathcal{L}|\rho\rangle$. A variational quantum algorithm based on the idea was proposed in the next year,[27] where the Hermitian and positive semi-definite properties of the density matrix can be preserved through special design of the ansatz. However, the inconsistency of the re-normalization condition between the two types of representation encounters problems with quantum measurement. In fact, the natural idea for this problem would be to generate the purification of the mixed state with PQC. However, transforming the Lindbladian into local operators, and therefore efficiently measuring the cost function with quantum hardware, becomes difficult. Here, we make a step further to propose the variational quantum algorithms for the steady states in open quantum systems, where we use PQC to generate the purification of the mixed state. The cost function is based on the Lindblad master equation and it can be decomposed into a sum of polynomial number of terms. Each of these terms can be evaluated using a swap test, Therefore, the cost function can be evaluated using the quantum circuits with linear scale of qubits and polynomial number of quantum gates required. We perform numerical simulations on the one-dimensional dissipative transverse field Ising model to show its reliability. Next, we first introduce the Lindblad master equation and the steady state. Then, we give the method and resource estimation. We will introduce the ansatz, give the decomposition of cost function, introduce the swap test, explain the optimization sketch of parameters, and estimate the quantum gate and qubit complexity. Furthermore, we give the result of the numerical simulation of the algorithm to show its correctness. A summary and discussion are given finally. The Lindblad Master Equation and the Steady State. An open quantum system considers the situation where the system interacts with its environments, in which case the state of the system is a mixed state and it can be represented with a density matrix: $$ \rho = \sum_i p_i |\phi_i\rangle \langle \phi_i|,~~ \tag {1} $$ where $\{|\phi_i\rangle\}$ are pure states in the ensemble and $\{p_i\}$ are their corresponding occurring probabilities. The system evolves with dissipation due to the interaction, which could be described by the Lindblad master equation[1] $$\begin{aligned} \frac{d\rho}{dt} ={}& \mathcal{L}\rho \\ :={}& -i[H,\rho] + \sum_i \gamma_i \Big(c_i\rho c_i^+ - \frac 12 c_i^+c_i\rho-\frac 12 \rho c_i^+c_i\Big), \end{aligned}~~ \tag {2} $$ where $\rho$ and $H$ are the system's state and Hamiltonian, respectively; $\gamma_i$ and $c_i$ are the dissipation rate and the jump operator for the $i$-th channel and $[H,\rho]=H\rho-\rho H$. The steady state satisfies $$ \mathcal{L}\rho_{\scriptscriptstyle{\rm SS}} = 0,~~ \tag {3} $$ which indicates that once the state of the system is $\rho_{\scriptscriptstyle{\rm SS}}$, it will be invariable because of $d\rho/dt=0$. We can view the steady state as $$ \rho_{\scriptscriptstyle{\rm SS}} = \lim_{t\to\infty} \rho(t).~~ \tag {4} $$ Therefore, we can get the steady state by repeated iterations of the master equation from an initial state. However, the complexity of Hamiltonian simulation for open quantum systems makes it unpractical.[28] Therefore, an approximate method is necessary. Here we will focus on the case where the steady state is unique. Indeed, this is acceptable for typical systems with finite dimension of Hilbert space.[29–31] Method—Ansatz. The ansatz selection is a core part of the variational quantum algorithm. With a set of variational parameters $\{\theta\}$ in the PQC $U(\theta)$ and an initial state $|0\rangle$, generally the ansatz is $|\psi(\theta)\rangle=U(\theta)|0\rangle$. An information-based ansatz, such as unitary coupled-cluster method in quantum chemistry simulations and quantum alternating operator ansatz in combinatorial optimization problems, can take the symmetry or other properties of the system into consideration. While in the cases where we know little information about the system, we require that the PQC has a strong expressive power, where the solution can be found when varying the parameters, even though a cost of resource for optimization is needed.
cpl-38-8-080301-fig1.png
Fig. 1. An illustration of the hardware-efficient ansatz used in this work with an example of $2n=4$. In the circuit, each layer consists of single-qubit rotation gates and two-qubit entangling gates. The structure is repeated $M$ times to construct $U(\theta)$ to which the parameter would be different for a specific gate in each layer. The ansatz is obtained by tracing out the bottom $n=2$ qubits.
Here, we use the hardware-efficient ansatz,[11,32,33] which combines easily implementable single- and double-qubit operations as one layer. The layers are repeated to increase the expressive power. Here, the structure we choose to construct the ansatz is $$ U(\theta) = \prod_{i=1}^{M} \prod_{j=1}^{N} U_{\rm ent}(\theta_i)R_z(\theta_{ij1})R_x(\theta_{ij2})R_z(\theta_{ij3}),~~ \tag {5} $$ where $N$ and $M$ are the number of qubits and layers in the circuit, respectively. For the entangling part $U_{\rm ent}$ we choose sequential control-$R_y$ operations. The illustration of the PQC we use is shown in Fig. 1 with a system of two qubits as an example. Considering the $n$-qubit mixed state formulized as Eq. (1), we can view the system and the environment as a closed system, where the system-environment joint state is the purification of the system's state: $$ |\varPhi\rangle = \sum_i \sqrt{p_i} |\phi_i\rangle \otimes |e_i\rangle,~~ \tag {6} $$ such that we can obtain the state of the system by tracing out the environment as $$ Tr_E(|\varPhi\rangle \langle\varPhi|) = \rho. $$ Note that for an $n$-qubit mixed state, we can always use another $n$ auxiliary qubits to construct its purification.[34] Thus, we can first construct a $2n$-qubit pure state $$ |\varPsi(\theta) \rangle= U(\theta)|0\rangle^{\otimes 2n},~~ \tag {7} $$ and then trace out the $n$ auxiliary qubits to obtain the ansatz $\rho(\theta)$. Cost Function. After defining the ansatz of target density matrix, here we introduce the cost function $\mathcal{C}(\theta)$. The steady state satisfies Eq. (3), which indicates that in the Linbblad form for a density matrix $\rho=\sum_{ij}\rho_{ij}|i\rangle\langle j|$, we have $$ \frac{d}{dt}\rho =\mathcal{L}\rho = \sum_{ij} \frac{d}{dt} \rho_{ij} |i\rangle\langle j|.~~ \tag {8} $$ This shows that if $\mathcal{L}\rho=0$, then we have $\frac{d}{dt}\rho_{ij}=0$ for all indices $i$ and $j$, such that the state evolution under the Lindblad master equation will be invariable. Therefore, we define the cost function based on the Frobenius norm of $\mathcal{L}\rho$: $$ \mathcal{C}(\theta) = ||\mathcal{L}\rho(\theta)||_F^2 = \sum_{ij} |[\mathcal{L}\rho(\theta)]_{ij}|^2.~~ \tag {9} $$ It is obvious that $\mathcal{C}(\theta)\geq 0$ for all $\theta$ and $\mathcal{C}(\theta)=0\Leftrightarrow \rho=\rho_{\scriptscriptstyle{\rm SS}}$. We observe that the Frobenius norm of a matrix can be transformed to $$ \mathcal{C} (\theta) = || \mathcal{L} \rho ||_F^2 = {\rm tr}[ (\mathcal{L} \rho)^†\mathcal{L} \rho ]. $$ According to Eq. (2), where the function of Lindbladian can be viewed as a sum of terms, each of which has a unitary operator multiplying on the left and/or right of the density matrix. We can rewrite the function of the Lindbladian as $$\begin{align} \mathcal{L} \rho &= \sum_i f_i U_i \rho V_i,\\ (\mathcal{L} \rho)^† &= \sum_i f_i^* V_i^† \rho U_i^†, \end{align} $$ where $f_i$ are complex numbers and $U_i/V_i$ are unitary operators. Then, the cost function can be decomposed as $$ {\rm tr}[ (\mathcal{L} \rho)^†\mathcal{L} \rho ] = tr\Big(\sum_{ij} f_i^*f_j V_i^† \rho U_i^† U_j\rho V_j\Big). $$ Since ${\rm tr}(AB)={\rm tr}(BA)$, finally the cost function is transformed to the following form: $$ \mathcal{C} (\theta) = \sum_{ij} f_i^*f_j {\rm tr}(\rho U_i^† U_j \rho V_j V_i^†).~~ \tag {10} $$ Swap Test. In Eq. (10), every term has the form of ${\rm tr}(\rho U\rho V)$, which can be estimated using the swap test. A similar circuit has been used in Ref. [16]. The swap test can be applied to obtain the value of inner product terms, which is frequently used to evaluate functions like quantum state fidelity. The quantum circuit of the swap test to evaluate the cost function for this work is shown in Fig. 2. The auxiliary qubit is measured at the end of the circuit and we have the following relation: $$ \Re[ {\rm tr}(\rho_1V\rho_2U) ] = 2p(0) - 1,~~ \tag {11} $$ where $\Re[\cdots]$ denotes the real part of a number. When we want to obtain the imaginary part we can add a phase gate $P=\begin{pmatrix} 1 & 0 \\ 0 & i \end{pmatrix}$ after the first $H$ gate, then the relations will be $$ \Im[ {\rm tr}(\rho_1V\rho_2U) ] = 1-2p(0). $$
cpl-38-8-080301-fig2.png
Fig. 2. The quantum circuit of the swap test to measure terms in the form of ${\rm tr}(\rho U\rho V)$. Desired values are associated with the probability of obtaining $|0\rangle$ when measuring the auxiliary qubit at the end of the circuit. When evaluating the real part, the red phase gate $P$ is not applied and we have $\Re[ {\rm tr}(\rho_1V\rho_2U) ] = 2p(0) - 1$. While for the imaginary part, the gate $P$ is applied with the relation changed to $\Im[ {\rm tr}(\rho_1V\rho_2U) ] = 1-2p(0) $.
Parameter Optimization. In our numerical simulation, we use the gradient-free method, Nelder–Mead (NM) algorithm,[35] which has been used in some variational quantum algorithms.[12,36] The NM method starts with a series of randomly initial points, then obtains a new point to replace the worst point. These steps stop when the change tolerance is satisfied or the maximum iteration times is achieved. The method written in Scipy[37] has a default value of the maximum iteration steps: $200\times d$, where $d$ is the number of parameters. We find that instead of changing this default number, repeating this algorithm for several times obtains a better result. Therefore, in our test, we repeat the algorithms for several times to achieve sufficient optimization. It has to be pointed out that optimizing the parameters to minimize the cost function is quite difficult for variational quantum algorithms, which has been proven to be NP-hard.[38] The Barren plateau phenomenon has been studied for randomly selected PQC structure and initialized parameters in large quantum systems.[39] It was also shown that the phenomenon is related to the locality of the cost function[40] and the expressibility of the ansatz.[41] In addition, finite sampling times and quantum gate errors[42] make the value of measured cost function stochastic, which cause problems with the optimization process. In fact, several optimization methods for stochastic optimization[43] and BP[44] have been proposed and one of our further works will be to focus on the optimization problems. Resource Estimation—Qubits. For an $n$-qubit system, we use $2n$ qubits to construct the purification. In the swap test, we need two copies and one auxiliary qubit. The total number of qubits required to evaluate the cost function is $4n+1$. Single-qubit measurement is enough in the swap test. Number of Terms in the Cost Function. Besides the qubits, the number of terms for measuring is also important for this algorithm, which greatly influences the cost function evaluation efficiency and sampling accuracy. Suppose that the Hamiltonian and $m$ dissipative operators can be expressed as a linear combination of $q_h$ and $q_l$ easily implementable unitary operators [$q_h$ and $q_l$ are assumed to belong to $O[{\rm poly}(n)]$ naturally]. The totally number of terms is then about $$ O[ (2q_h +3mq_l)^2 ]. $$ Here we consider that $m\in O[{\rm poly}(n)]$, then the total number, is also $O[{\rm poly}(n)]$, which can be efficiently measured. Gates and Parameters. In the hardware efficient ansatz with an $n$-qubit system and $M$-layers, there are in total $8nM$ parameters. Since each control-$R_y$ gate can be decomposed into two CNOT gates and two $R_y$ gates, in total $5nM$ single-qubit rotation gates and $2nM$ CNOT gates are required. In the swap test, the $O(n)$ control-$U/V$ gates can be decomposed into single- and double-qubit gates in the same order. While the control-SWAP gate can be realized with $n$ C-SWAP gates, where each one is equal to three Toffoli gates. Therefore, the total number of quantum gates needs scales linearly with the size of the system. All in all, we can see that the resource to measure the cost function is acceptable. Result. We perform numerical simulations on the one-dimensional dissipative transverse field Ising model. The code is written is Python. We use Qiskit[45] and Scipy packages to process the data and optimize the parameters. Cost functions are evaluated using “state-vector simulators” in Qiskit, where no sampling error or gate error is considered. The computational package for dynamics of open quantum systems, QuTiP,[46] is used to numerically obtain the steady state as a benchmark. We use the fidelity between our optimal ansatz and the steady state obtained with QuTiP to compare their difference. The fidelity of two density operators $\sigma$ and $\rho$ is given by $$ F(\sigma,\rho) = \Big({\rm tr} \sqrt{ \sqrt{\sigma} \rho \sqrt{\sigma} }\Big)^2.~~ \tag {12} $$ The Hamiltonian of the model is $$ H =\frac V4 \sum_{i=0}^{L-1} \sigma_i^z \sigma_{i+1}^z + \frac g2 \sum_{i=0}^{L-1}\sigma_i^x,~~ \tag {13} $$ where $L$ is the length of the spin chain with $L+1$ spin particles, labeled as $0,1,\ldots,L$, $\sigma^a = \frac 12 S^a$ with $a=\{x,y,z\}$ is the Pauli operator, $V$ and $g$ models the nearest-neighbor interaction strength in $z$-axis and the amplitude of transverse field along the $x$-axis. The periodic boundary condition is applied, which indicates that $\sigma_L=\sigma_0$. We assume that there is only one quantum channel on every qubit. The jump operator and dissipation rate are defined as follows: $$ c_i = \sigma_i^- = \frac 12(\sigma_i^x - i\sigma_i^y),~~ \gamma_i = \gamma.~~ \tag {14} $$
cpl-38-8-080301-fig3.png
Fig. 3. Data in the optimization process for the dissipative one-dimensional transverse field Ising model. (a) The values of loss function with respect to the iteration times. The total iteration times is about $8\times 10^4$ and the optimal loss function is about $1.8\times 10^{-3}$. (b) The values of fidelity between the ansatz and the stationary state density matrix obtained with QuTiP in the optimization process. The fidelity increases along with the cost function decreasing. The optimal fidelity is over 99.8%, which can show the validity of our method.
We performed our test of this model with the model parameters as $L=4$, $V=0.3$, $g=1$ and $\gamma=0.5$, respectively. We use four auxiliary qubits according to former analysis. The number of layers in the hardware-efficient ansatz is 4. The data in the optimization process is shown in Fig. 3. In Fig. 3(a), we plot the change of loss function with respect to the iteration times. The total iteration times is about $8\times 10^4$ and we reach the optimal loss function at about $1\times 10^{-3}$. We plot the change of fidelity of our ansatz density matrix with the one obtained with QuTiP in Fig. 3(b). Obviously we can see that the fidelity increases in the optimization process and the optimal fidelity is over 99.8%. In Fig. 4, we plot the image pictures of the steady state density matrix obtained with our method and the one from QuTiP. There is no visible difference, which can show the correctness of our method.
cpl-38-8-080301-fig4.png
Fig. 4. Image pictures of the stationary state density matrices obtained with QuTiP and out optimal ansatz. The left and the right parts represent the real and imaginary parts of the density matrix, respectively. There is no visible difference, which can also show the correctness of our method.
Conclusion and Discussion. We have proposed a variational quantum algorithm to find the steady state of open quantum systems. We use hardware-efficient ansatz to generate the purification of the ansatz for the mixed state. The cost function is based on the Lindbladian and is decomposed into a sum of polynomial number of terms, where each can be evaluated using a swap test with the quantum hardware. We perform numerical simulations on the one-dimensional dissipative transverse field Ising model. The optimal fidelity between the ansatz and the one obtained with QuTiP is over 99%, showing the reliability of our method. Just like many other variational quantum algorithms, the limited sampling times and the noise of quantum hardware have caused a lot of problems in the optimization process. Finding an efficient optimization method will be a possible future work for our algorithm. Due to the importance of solving steady states in open quantum systems and the preventive use of this type of loss function in various tasks, more applications would be found with this algorithm. We thank Professor Yong-Jian Han for helpful discussion. The numerical calculations in this work were carried out on the supercomputing system in the Supercomputing Center of University of Science and Technology of China. References On the generators of quantum dynamical semigroupsAn open-system quantum simulator with trapped ionsSimulation of single-qubit open quantum systemsObservation of a Dissipative Phase Transition in a One-Dimensional Circuit QED LatticeQuantum dynamical field theory for nonequilibrium phase transitions in driven open systemsQuantum Computing in the NISQ era and beyondNoisy intermediate-scale quantum (NISQ) algorithmsVariational Quantum AlgorithmsHybrid Quantum-Classical Algorithms and Quantum Error MitigationHardware-efficient variational quantum eigensolver for small molecules and quantum magnetsA variational eigenvalue solver on a photonic quantum processorScalable Quantum Simulation of Molecular EnergiesCloud Quantum Computing of an Atomic NucleusWitnessing eigenstates for quantum simulation of Hamiltonian spectraTheory of variational quantum simulationVariational ansatz-based quantum simulation of imaginary time evolutionVariational Quantum Simulation of General ProcessesHybrid quantum variational algorithm for simulating open quantum systems with near-term devicesEfficient Variational Quantum Simulator Incorporating Active Error MinimizationConstructing exact representations of quantum many-body systems with deep neural networksLatent Space Purification via Neural Density OperatorsVariational Quantum Monte Carlo Method with a Neural-Network Ansatz for Open Quantum SystemsNeural-Network Approach to Dissipative Quantum Many-Body DynamicsVariational Neural-Network Ansatz for Steady States in Open Quantum SystemsConstructing neural stationary states for open quantum many-body systemsVariational quantum algorithm for nonequilibrium steady statesDissipative Quantum Church-Turing TheoremOn the uniqueness of the steady-state solution of the Lindblad–Gorini–Kossakowski–Sudarshan equationStabilizing open quantum systems by Markovian reservoir engineeringNonequilibrium steady state in open quantum systems: Influence action, stochastic equation and power balanceExpressibility and Entangling Capability of Parameterized Quantum Circuits for Hybrid Quantum‐Classical AlgorithmsExpressive power of parametrized quantum circuitsQuantum computation and quantum informationNelder-Mead algorithmQuantum Chemistry Calculations on a Trapped-Ion Quantum SimulatorSciPy 1.0: fundamental algorithms for scientific computing in PythonTraining variational quantum algorithms is NP-hard – even for logarithmically many qubits and free fermionic systemsStructural absorption by barbule microstructures of super black bird of paradise feathersRobust and efficient hydrogenation of carbonyl compounds catalysed by mixed donor Mn(I) pincer complexesConnecting ansatz expressibility to gradient magnitudes and barren plateausNoise-Induced Barren Plateaus in Variational Quantum AlgorithmsAdam: A Method for Stochastic OptimizationAn initialization strategy for addressing barren plateaus in parametrized quantum circuitsQuTiP 2: A Python framework for the dynamics of open quantum systems
[1] Lindblad G 1976 Commun. Math. Phys. 48 119
[2] Barreiro J T, Müller M, Schindler P, Nigg D, Monz T, Chwalla M, Hennrich M, Roos C F, Zoller P, and Blatt R 2011 Nature 470 486
[3] Sweke R, Sinayskiy I, and Petruccione F 2014 Phys. Rev. A 90 022331
[4] Fitzpatrick M, Sundaresan N M, Li A C Y, Koch J, and Houck A A 2017 Phys. Rev. X 7 011016
[5] Marino J and Diehl S 2016 Phys. Rev. B 94 085150
[6] Preskill J 2018 Quantum 2 79
[7]Gambetta J 2019 APS March Meeting 2019 4–8 March 2019, Boston, Massachusetts, Vol. 64, No. 2
[8] Bharti K, Cervera-Lierta A, Kyaw T H, Haug T, Alperin-Lea S, Anand A, Degroote M, Heimonen H, Kottmann J S, Menke T et al. 2021 arXiv:2101.08448 [quant-ph]
[9] Cerezo M, Arrasmith A, Babbush R, Benjamin S C, Endo S, Fujii K, McClean J R, Mitarai K, Yuan X, Cincio L et al. 2020 arXiv:2012.09265 [quant-ph]
[10] Endo S, Cai Z, Benjamin S C, and Yuan X 2021 J. Phys. Soc. Jpn. 90 032001
[11] Kandala A, Mezzacapo A, Temme K, Takita M, Brink M, Chow J M, and Gambetta J M 2017 Nature 549 242
[12] Peruzzo A, McClean J, Shadbolt P, Yung M H, Zhou X Q, Love P J, Aspuru-Guzik A, and O'brien J L 2014 Nat. Commun. 5 4213
[13] O'Malley P J J, Babbush R, Kivlichan I D, Romero J, McClean J R, Barends R, Kelly J, Roushan P, Tranter A, Ding N, Campbell B, Chen Y, Chen Z, Chiaro B, Dunsworth A, Fowler A G, Jeffrey E, Lucero E, Megrant A, Mutus J Y, Neeley M, Neill C, Quintana C, Sank D, Vainsencher A, Wenner J, White T C, Coveney P V, Love P J, Neven H, Aspuru-Guzik A, and Martinis J M 2016 Phys. Rev. X 6 031007
[14] Dumitrescu E F, McCaskey A J, Hagen G, Jansen G R, Morris T D, Papenbrock T, Pooser R C, Dean D J, and Lougovski P 2018 Phys. Rev. Lett. 120 210501
[15] Santagati R, Wang J, Gentile A A, Paesani S, Wiebe N, Mcclean J R, Morley-Short S, Shadbolt P J, Bonneau D, and Silverstone J W 2018 Sci. Adv. 4 eaap9646
[16] Yuan X, Endo S, Zhao Q, Li Y, and Benjamin S C 2019 Quantum 3 191
[17] McArdle S, Jones T, Endo S, Li Y, Benjamin S, and Yuan X 2019 npj Quantum Inf. 5 75
[18] Endo S, Sun J, Li Y, Benjamin S C, and Yuan X 2020 Phys. Rev. Lett. 125 010501
[19] Mahdian M and Yeganeh H D 2020 J. Phys. A 53 415301
[20] Li Y and Benjamin S C 2017 Phys. Rev. X 7 021050
[21] Carleo G, Nomura Y, and Imada M 2018 Nat. Commun. 9 5322
[22] Torlai G and Melko R G 2018 Phys. Rev. Lett. 120 240503
[23] Nagy A and Savona V 2019 Phys. Rev. Lett. 122 250501
[24] Hartmann M J and Carleo G 2019 Phys. Rev. Lett. 122 250502
[25] Vicentini F, Biella A, Regnault N, and Ciuti C 2019 Phys. Rev. Lett. 122 250503
[26] Yoshioka N and Hamazaki R 2019 Phys. Rev. B 99 214306
[27] Yoshioka N, Nakagawa Y O, Mitarai K, and Fujii K 2020 Phys. Rev. Res. 2 043289
[28] Kliesch M, Barthel T, Gogolin C, Kastoryano M, and Eisert J 2011 Phys. Rev. Lett. 107 120501
[29] Nigro D 2019 J. Stat. Mech.: Theory Exp. 2019 043202
[30] Schirmer S G and Wang X 2010 Phys. Rev. A 81 062306
[31] Hsiang J T and Hu B 2015 Ann. Phys. 362 139
[32] Sim S, Johnson P D, and Aspuru-Guzik A 2019 Adv. Quantum Technol. 2 1900070
[33] Du Y, Hsieh M H, Liu T, and Tao D 2020 Phys. Rev. Res. 2 033125
[34] Nielsen M A and Chuang I 2002 Quantum Computation and Quantum Information (Cambridge: Cambridge University Press)
[35] Singer S and Nelder J 2009 Scholarpedia 4 2928
[36] Hempel C, Maier C, Romero J, McClean J, Monz T, Shen H, Jurcevic P, Lanyon B P, Love P, Babbush R, Aspuru-Guzik A, Blatt R, and Roos C F 2018 Phys. Rev. X 8 031022
[37] Virtanen P, Gommers R, Oliphant T E et al. 2020 Nat. Methods 17 261
[38] Bittel L and Kliesch M 2021 arXiv:2101.07267 [quant-ph]
[39] McClean J R, Boixo S, Smelyanskiy V N, Babbush R, and Neven H 2018 Nat. Commun. 9 1
[40] Cerezo M, Sone A, Volkoff T, Cincio L, and Coles P J 2021 Nat. Commun. 12 12
[41] Holmes Z, Sharma K, Cerezo M, and Coles P J 2021 arXiv:2101.02138 [quant-ph]
[42] Wang S, Fontana E, Cerezo M, Sharma K, Sone A, Cincio L, and Coles P J 2020 arXiv:2007.14384 [quant-ph]
[43] Kingma D P and Ba J 2014 arXiv:1412.6980 [cs.LG]
[44] Grant E, Wossnig L, Ostaszewski M, and Benedetti M 2019 Quantum 3 214
[45] Treinish M, Gambetta J, Nation P et al. 2021 Qiskit/qiskit: Qiskit 0.28.0
[46] Johansson J R, Nation P D, and Nori F 2013 Comput. Phys. Commun. 184 1234