Theory of variational quantum simulation

The variational method is a versatile tool for classical simulation of a variety of quantum systems. Great efforts have recently been devoted to its extension to quantum computing for efficiently solving static many-body problems and simulating real and imaginary time dynamics. In this work, we first review the conventional variational principles, including the Rayleigh-Ritz method for solving static problems, and the Dirac and Frenkel variational principle, the McLachlan's variational principle, and the time-dependent variational principle, for simulating real time dynamics. We focus on the simulation of dynamics and discuss the connections of the three variational principles. Previous works mainly focus on the unitary evolution of pure states. In this work, we introduce variational quantum simulation of mixed states under general stochastic evolution. We show how the results can be reduced to the pure state case with a correction term that takes accounts of global phase alignment. For variational simulation of imaginary time evolution, we also extend it to the mixed state scenario and discuss variational Gibbs state preparation. We further elaborate on the design of ansatz that is compatible with post-selection measurement and the implementation of the generalised variational algorithms with quantum circuits. Our work completes the theory of variational quantum simulation of general real and imaginary time evolution and it is applicable to near-term quantum hardware.

The variational method is a versatile tool for classical simulation of a variety of quantum systems. Great efforts have recently been devoted to its extension to quantum computing for efficiently solving static many-body problems and simulating real and imaginary time dynamics. In this work, we first review the conventional variational principles, including the Rayleigh-Ritz method for solving static problems, and the Dirac and Frenkel variational principle, the McLachlan's variational principle, and the time-dependent variational principle, for simulating real time dynamics. We focus on the simulation of dynamics and discuss the connections of the three variational principles. Previous works mainly focus on the unitary evolution of pure states. In this work, we introduce variational quantum simulation of mixed states under general stochastic evolution. We show how the results can be reduced to the pure state case with a correction term that takes accounts of global phase alignment. For variational simulation of imaginary time evolution, we also extend it to the mixed state scenario and discuss variational Gibbs state preparation. We further elaborate on the design of ansatz that is compatible with post-selection measurement and the implementation of the generalised variational algorithms with quantum circuits. Our work completes the theory of variational quantum simulation of general real and imaginary time evolution and it is applicable to near-term quantum hardware.

I. INTRODUCTION
Variational simulation is a widely used technique in many-body physics [1][2][3][4][5] and chemistry [6][7][8]. As the full Hilbert space of a quantum system grows exponentially with respect to the system size, direct classical simulation of a large many-body state is generally impossible. The variational method avoids the exponential space problem by considering trial states from a physically motivated small subset of the exponentially large Hilbert space [9,10]. Variational classical simulation (VCS) works for both static problems in finding the ground state and ground state energy, and dynamical problems in simulating the real and imaginary time evolution of quantum states [3,4,[10][11][12][13][14]. In modern science, the variational method has become a versatile tool for simulating various problems when the target system state can be well modelled classically.
Nevertheless, there also exist many problems that may not be solved classically even with the classical variational method [15][16][17]. This is because there exist highly entangled many-body states that may not be efficiently represented via any classical method. This problem motivates the development of quantum simulation with universal quantum computers [18], because quantum systems can be efficiently encoded or represented by another quantum system and real time evolution of the Schrödinger equation can be efficiently realised via a unitary quantum circuit [19,20]. However, implementing a universal quantum computer requires the coherent and accurate control of millions of qubits [21][22][23]; While the state-of-the-art quantum hardwares can only accurately control tens of qubits [17,24]. Although this number may be extended by one or two orders in the near future, realising a universal fault tolerant quantum computer remains a challenging task.
With noisy intermediate-scale quantum hardware [25], quantum advantages may still be achieved in many tasks with recently proposed hybrid or, more specifically, variational quantum simulation (VQS) methods [26][27][28][29][30][31][32][33][34][35][36][37][38][39]. By dividing the problem into two levels, VQS methods only use the quantum processor to solve the core and classically intractable problem and leave the relatively easy tasks to a classical computer. For solving static problems, the Rayleigh-Ritz variational method can be naturally generalised to the quantum regime, such as the variational quantum eigensolver (VQE) [27][28][29][30][31][32][33][34][35][36][37][38] and the quantum approximate optimisation algorithm (QAOA) [26]. Considering parameterised trial states prepared by quantum circuits [40][41][42], the target static problem is mapped to a cost function of the measurement outcomes of the quantum state, which is further optimised via a classical algorithm. As the state preparation and measurement procedures may be hard to simulate by a classical computer [15][16][17], hybrid or variational quantum simulation can have an intrinsic quantum advantage over VCS and it has been proposed as the standard tool for studying quantum computational chemistry [43,44].
A broader and harder challenge is the simulation of dynamics. Classically, there are three variational principles-the Dirac and Frenkel variational principle [45,46], the McLachlan variational principle [47], and the time-dependent variational method principle [48,49]. Recently, the variational methods have been generalised to simulate real [39] and imaginary [50,51] time evolution of pure quantum systems. By encoding quantum states with a parameterised circuit, the evolution of the state can be mapped to the evolution of the parameters controlling the circuit. However, different variational principles can lead to different evolution equations of parameters. It has not been clear which variational principle or which evolution equation should be used in practice. Furthermore, previous works only focus on the evolution of pure states of closed systems, rather than the more general case of the stochastic evolution of mixed states of open systems. While practical quantum systems always suffer from unavoidable interaction with the environment, it makes almost every quantum system open. Therefore, simulating the evolution of open quantum systems can be useful for studying realistic quantum systems.
In this work, we complete the VQS theory by studying the equivalence and difference of the variational principles and extending the principles to the general stochastic evolution of mixed states. We first review the theory of variational simulation in Sec. II. Then we study the equivalence and difference of the variational principles and the derived evolution equations in Sec. III. We consider general real time stochastic evolution of mixed states in Sec. IV. As an example, we re-derive the evolution equation of parameters for unitary evolution of pure states and find a correction term which applies to previous results and handles the global phase alignment. After applying all the three variational principles to the general mixed state case, we find that McLachlan's variational principle is the most appropriate variational principle that leads to a consistent theory of VQS. In Sec. V, we study variational simulation of imaginary time evolution of mixed states and discuss variational Gibbs state preparation. We further elaborate on the implementation with quantum circuits in Sec. VI. We conclude our work and discuss future directions in Sec. VII.

A. Static problem: Rayleigh-Ritz method
The Rayleigh-Ritz method is designed to solve static problems of finding the ground state energy E 0 of a given Hamiltonian H = j h j σ j . Here, we assume that H is a linear combination of tensor products of local operators σ j with coefficients h j . The ground state energy is the solution of the problem, where the minimisation is over states from the whole Hilbert space. However, the size of the Hilbert space grows exponentially to the system size, it is computationally hard to search the whole Hilbert space. Instead, the Rayleigh-Ritz method only searches states from a subset of the whole Hilbert space to find an approximate solution. For example, classically, one can consider trial states that are linear combinations of a small number of basis states {|φ i }, with a = (a 1 , a 2 , . . . ). Suppose the ground state can be approximated by the trial state |φ( a) with a proper vector of coefficients a, then the ground state energy can be approximated or upper bounded by The minimisation problem can be solved by calculating the partial derivative of E est 0 over a * i , and is equivalent to the minimal solution to where H and S are the matrix of φ i | H |φ j and φ i |φ j , respectively. We refer to the trial state also as the ansatz. In practice, different ansatze have been invented in condensed matter physics and computational chemistry [4,9,[41][42][43]. When the minimisation cannot be solved analytically, one can also numerically solve it with classical optimisation algorithms. With quantum computers, the trial state |φ( θ) can be prepared by applying a sequence of parameterised gates R k (θ k ) to an initial state |0 as |φ( θ) = R N (θ N ) . . . R k (θ k ) . . . R 1 (θ 1 ) |0 . Here, R k (θ k ) is the k th gate controlled by the real parameter θ k and θ = (θ 1 , θ 2 , . . . , θ N ). The ansatz is automatically normalised and the average energy φ( θ)|H|φ( θ) can be obtained by measuring each term φ( θ)|h i |φ( θ) and linearly combining the measurement outcomes. In Sec. VI, we also introduce other possible ways of realising the trial states. To approximate the ground state, we can minimise φ( θ)|H|φ( θ) over the parameter space with a classical optimisation algorithm. For instance, by calculating the gradient of φ( θ)|H|φ( θ) , a local minimum of φ( θ)|H|φ( θ) can be found via the gradient descent algorithm. Such a hybrid algorithm for solving the eigenvalue of the Hamiltonian is called variational quantum eigensolver (VQE) [27][28][29][30][31][32][33][34][35][36][37][38].

B. Real dynamics: three variational principles
Now we introduce three variational principles for simulating real time dynamics of pure quantum states [3,39]. We only review the results here and leave the derivation in Appendix A. To simulate Schrödinger's equation we consider a parametrised trial state |φ( θ(t)) with time dependent parameters. parameters can be complex for classical simulation but it suffices for them to be real for quantum simulation. This is because without loss of generality the trial state can be prepared by a quantum circuit by applying unitary gates R k (θ k ) that is in the form of e iθ k σ with Hermitian operator σ and real parameter θ k . Suppose the state at time t can be represented by the trial state |ψ(t) = |φ( θ(t)) , then we want to approximate the state |ψ(t + δt) at time t + δt by |φ( θ(t + δt)) . According to Schrödinger's equation, the state at time t + δt should be |ψ(t + δt) = |φ( θ(t)) − iδtH |φ( θ(t)) , which, however, may be impossible to be represented by the trial state with any parameters. A variational method is used to project the state |ψ(t + δt) back to the trial state manifold. Then the evolution of the state |ψ under the Schrödinger equation can be mapped onto the trial state manifold and thus we determine the evolution of parameters. a. The Dirac and Frenkel variational principle. Replace |ψ(t) with |φ( θ(t)) , then Schrödinger's equation becomes The L.H.S is in the tangent subspace { ∂|φ( θ(t)) ∂θj }, but the R.H.S. −iH |φ( θ(t)) may not be in this subspace. So the Dirac and Frenkel variational principle [48,49] directly projects the equation onto the subspace by δφ( θ(t)) ∂ ∂t + iH φ( θ(t)) = 0, where δφ( θ(t))| = i ∂ φ( θ(t))| ∂θi δθ i . The evolution of parameters can be solved to be where A and C are In general, A and C are complex, which leads to a complex solutionθ j . A complex solution is not problematic in classical simulation as parameters are complex. However, parameters in quantum simulation are generally real, a complex solution thus leads to a state outside of the ansatz space. Naively, one can take the real or imaginary part of Eq. (8) to make the solution real. Although this seems artificial, we show that the corresponding equations can be obtained from the McLachlan's variational principle and the time-dependent variational principle, respectively.
b. McLachlan's variational principle. The evolution of the trial state is ∂ |φ( θ(t)) /∂t and the evolution according to the Schrödinger equation is −iH |φ( θ(t)) . The McLachlan's variational principle [47] is to minimise the distance between them, δ (∂/∂t + iH) |φ( θ(t)) = 0. (10) Here |ψ = ψ|ψ is the norm of |ψ . When the parameters θ are complex, the solution gives the same evolution of Eq. (8) under the Dirac and Frenkel variational principle. However, assuming θ to be real, the equation becomes which can be obtained by taking the real part of Eq. (8).
Here A R i,j is the real part of A i,j and C I i is the imaginary part of C i . As shown in Sec. III, this equation is equivalent to Eq. (8) when the trial state |φ( θ(t)) is powerful enough to represent the target state |ψ(t) . The advantage of this equation is that it always gives a real solution ofθ, which makes it consistent with VQS with real parameters. Furthermore, we will show shortly that McLachlan's variational principle is the most consistent variational principle that works similarly to general stochastic evolution.
However, it is worth noting that the evolution according to Eq. (11) is not the ultimate one owing to the global phase mismatch between the target state |ψ(t) and the trial state |φ( θ(t)) . Suppose the state at time t can be always represented by the trial state, |ψ(t) = e iθ N +1 (t) |φ( θ(t)) , up to a global phase θ N +1 (t). Although global phase is physically irrelevant, it plays an important role in the variational method. This is because even if |ψ(t) and |φ( θ(t)) are physically equivalent, the derivatives ∂ |ψ(t) /∂t and ∂ |φ( θ(t)) /∂t can still differ significantly. The derivative of the state |ψ(t) is ∂t .
A further point of interest concerns our ability to track the accuracy of the simulation. In practice, after solving the derivatives ofθ in Eq. (11), we can also get the distance between the true evolution and the evolution of the trial state, Here H 2 = φ( θ(t))|H 2 |φ( θ(t)) . Therefore, by measuring A R , C I , and H 2 , we can calculate the distance to verify the variational algorithm. c. The time-dependent variational principle. The Schrödinger equation can be obtained from the Lagrangian L = ψ(t)|(∂/∂t + iH)|ψ(t) . Replacing |ψ(t) with |φ( θ(t)) , the Lagrangian becomes For complex θ i , parameters are θ i , θ * i ,θ i , andθ * i . After applying the Euler-Lagrange equation, the evolution of parameters can be derived and shown to coincide with Eq. (8). For real θ i , parameters of the Lagrangian are θ i andθ i , and the evolution of parameters are which can be obtained by taking the imaginary part of Eq. (8). Here A I i,j is the imaginary part of A i,j and C R i is the real part of C i . As shown in Sec. III, this evolution equation is equivalent to Eq. (8) when the trial state |φ( θ(t)) is powerful enough to represent the target state |ψ(t) . However, as we take the imaginary part of the A matrix, the A I matrix is more likely to be singular, making the evolution unstable.

C. Implementation
In classical simulation, the A matrix and the C vector are calculated classically. Here we show how to evaluate A and C with a quantum circuit introduced by Li and Benjamin [39]. Suppose the parameterised trial state is prepared by |φ( θ) = R |0 with R = R N (θ N ) . . . R k (θ k ) . . . R 1 (θ 1 ). In conventional quantum computing, each unitary gate involves a small subset of qubits, typically one or two. Therefore, the derivative of the unitary can be efficiently decomposed via where σ k,i are also unitary operators and f k,i are complex coefficients. For most single-and two-qubit-gates R k , we find that there are only one or two terms in this expression, and σ k,i is also a one-qubit or two-qubit gate. For example, when R k (θ k ) = e −iθ k σ/2 with a one-or two-qubit Hermitian operator σ, we have ∂R k (θ k )/∂θ k = −i/2 · σ · e −iθ k σ/2 with f k = −i/2 and σ k = σ in Eq. (16). Using the expression (16), we rewrite the derivative of the state as where Then, based on the definition of the Hamiltonian, the elements of A and C can be expressed as and The real and imaginary part of expressions (19) and (20) are in the form where the amplitude a and phase θ are determined by the coefficients, and U is a unitary operator either R † k,i R q,j or R † k,i σ j R. Such a term can be evaluated with the quantum circuit shown in Fig. 1. This circuit needs an ancillary qubit initialised in the state (|0 + e iθ |1 )/ √ 2 and a register initialised in the state |Ψ 0 . The ancillary qubit is measured after a sequence of unitary operations on the register and two controlled unitary operations, where the ancillary qubit is the control qubit. The value is given by Recently, Mitarai and Fujii [53] also showed a way to measure A and C without the ancillary qubit, making the realisation more friendly to near-term quantum hardware.

III. EQUIVALENCE OF THE THREE VARIATIONAL PRINCIPLES
With complex parameters of the trial state, the three variational principles lead to the same evolution equation of parameters. We will see that this is also true for general stochastic evolution and imaginary time evolution of mixed states. However, when considering real parameters as in quantum simulation, we can get three different equations (8), (11), and (15) for the evolution of parameters. With classical ansatz, a necessary condition for the equivalence of the three equations is discussed by Broeckhove et al. [49]. That is these three equations are equivalent when for any parameter θ i , there always exists a parameter θ j , such that However, such a condition is hard to be satisfied in quantum simulation with real parameters. Generally, when |φ( θ(t)) can uniquely represent any state |ψ(t) at time t, Eq. (8) and Eq. (11) both give real solutions and ancilla register redundant gates Quantum circuit for the evaluation of coefficients in the variational pure-state simulator (from Ref. [39]). The ancillary qubit is initialised in the state (|0 + e iθ |1 )/ √ 2 and measured in the |± = (|0 ± |1 )/ √ 2 basis. The probability P+ that the qubit is in the state |+ is ( Here, U k is one of σ k,i , and Uq is one of σq,j or σj (By taking q = N + 1, σj is put on the left side of RN in the product), assuming that k < q. We would like to remark that this circuit is a variant of the circuit proposed in Ref. [52]. This circuit contains N gates on the register, two flip gates (X) on the ancillary qubit, and two controlled unitary gates between the ancillary qubit and the register. Note that gates on the register after the second controlled unitary gate can be omitted. Usually, U k and Uq are unitary operators on only one or two register qubits rather than the entire register. This circuit can be further reduced to a direct measurement on the register qubits without the ancilla [53].
are also equivalent. The necessary and sufficient condition for the equivalence of Eq. (8) and Eq. (11) only requires the solution of Eq. (8) to be real.
Broadly speaking we may expect that Eq. (11) and Eq. (15) each suffice to define the evolution of our parameters. Practically, there are differences however. In particular, we can consider a single qubit system with Hamiltonian σ x and initial state |0 . Suppose the trial state is chosen to be e −iθσx |0 , which can exactly represent the evolution of the state. For such a simple case, the evolution equation (15) based on the time-dependent variational principle cannot work as A I and C R are always zero; nevertheless, Eq. (8) and Eq. (11) can both simulate the evolution. In fact, as the diagonal elements of A are in general real constants, the diagonal elements of A R and A I are real constants and zero, respectively. As the diagonal elements of A I are all zero, it makes the inverse of A I unstable when the off-diagonal elements are also close to 0. This suggests that the evolution equation (15) may not be the ideal choice for VQS.
In the following, we will re-derive the evolution equations from the mixed state perspective. Surprisingly, we will show that Eq. (15) from the time-dependent variational principle cannot be obtained in the mixed state case, whereas a variant of Eq. (11) can be consistently derived. Therefore, our work suggests that McLachlan's principle is the most consistent variational principle in variational quantum simulation.

IV. REAL TIME EVOLUTION
In this section, we extend VQS to mixed states and general stochastic evolutions, which can describe open quantum systems and noisy quantum hardware. Note that this approach differs from that in Ref. [54] where we can variationally simulate the stochastic mater equation.

A. General evolution
The general stochastic evolution of an open quantum system can be modelled by where the superoperator L(ρ) can be decomposed as with unitary operators S i and T i , and coefficients g i . Suppose ρ = ρ( θ(t)) is a parameterised state, we show how to transform the evolution of ρ into the evolution of parameters θ. Here, we mainly focus on the evolution of parameters and refer to Sec. VI for the design of the mixed states ansatz with unitary circuits. a. The Dirac and Frenkel variational principle. Let ρ = ρ( θ(t)), the general stochastic evolution becomes Because ∂ρ( θ(t)) ∂θi }, we can project the equation onto this tangent subspace by following the Dirac and Frenkel variational principle, The evolution of parameters is with elements of M and V defined by In classical simulation, θ can be complex and ∂θj . In this case, M , V , and the solutionθ can all be complex. However, considering quantum simulation with real parameters θ, we have ∂θj . Therefore, M and V are real and the solutionθ is also real. This is different from the pure state case where even if parameters are real, the solutioṅ θ can still be complex. Each term of M and V can be calculated with a quantum circuit that is discussed in Sec. VI.
b. McLachlan's variational principle. We can also minimise the distance between the evolution of the trial state and the true evolution, where ρ = Tr[ρ 2 ]. Then we can get the same result with the Dirac and Frenkel variational principle. After solving the derivatives ofθ, we can also get the distance between the true evolution and the evolution of the trial state, (29) Therefore, by measuring M , V , and Tr[L(ρ) 2 ], we can calculate the distance to verify the variational algorithm.
c. The time-dependent variational principle. The Lagrangian for the general evolution is For general complex parameters, θ and θ * are regarded as independent parameters, and the Euler-Lagrange equation applied to θ * gives the same evolution Eq. (26). However, when parameters are real with θ = θ * , the Euler-Lagrange equation cannot lead to any equation of parameters.
From the above result, we find that the three variational principles are equivalent for complex parameters. However, when parameters are real, the Dirac and Frenkel and the McLachlan variational principles are also equivalent, while the time-dependent variational principle cannot lead to a nontrivial evolution of parameters. Therefore, the McLachlan variational principle always gives a consistent real evolution of parameters in VQS.
In this case, M is defined in Eq. (27) and V is explicitly given by Suppose ρ( θ(t)) is a pure state |φ( θ(t)) , it becomes the Schrödinger equation and we have The evolution equation obtained from directly applying the McLachlan's principle to pure states with unitary evolution is j A R i,jθ j = C I i , as given in Eq. (11).
Compared to the A R matrix and C I vector, there are additional terms that estimate the overlap between the trial state and the derivative of the trial state in M and V . The terms ∂ φ( θ(t))| ∂θi |φ( θ(t)) can be measured with the circuit in Fig. 1 and the average energy φ( θ(t))| H |φ( θ(t)) can be directly measured. As we have discussed in Sec. II, the additional terms take account of the global phase difference between the trial state and the target evolution state. Consider a trial state e iθ N +1 |φ( θ(t)) with N parameters in |φ( θ(t)) and an extra parameter θ N +1 for the global phase. Using Eq. (11), the evolution of the first N parameters of e iθ N +1 |φ( θ(t)) is exactly given by Eq. (33). Because the global phase θ N +1 is irrelevant to measurement, we can thus only evolve the first N parameters with Eq. (33). However, the extra parameter θ N +1 is important for the global phase alignment between the trial state and the exact state as we easily confirm below.
Here, we consider a single-qubit example to illustrate the difference between Eq. (33) and Eq. (11). Suppose the trial state is prepared by |φ( θ(t)) = R z (θ z )R x (θ x ) |0 with two parameters θ z and θ x . Here, R x (θ) = e −iθσx/2 and R z (θ) = e −iθσz/2 , with Pauli matrices σ x and σ z . Suppose the system Hamiltonian is Pauli matrix σ y , we compare the simulation of the time evolution e −iσyt with randomly initialised state, as shown in Fig. 2. We find that the evolution under Eq. (11) generally fails to simulate the exact real time dynamics. Nevertheless, as the trial state can represent all pure qubit state up to a global phase, the evolution under Eq. (33) can indeed simulate the true time evolution e −iσyt . Therefore, Eq. (33) is recommended for variational simulation of real time unitary dynamics of pure states.

V. IMAGINARY TIME EVOLUTION
A. Pure state Next, we consider variational simulation of imaginary time evolution. We first focus on the pure state case studied in Ref. [50]. Replacing real time with imaginary time τ = it, the state at time τ is |ψ(τ ) = e −Hτ |ψ(0) ψ(0)|e −2Hτ |ψ(0) , and the Wick-rotated Schrödinger equation is, where E τ = ψ(τ )|H|ψ(τ ) is the expected energy at time τ . Consider a normalised parametrised trial state |φ( θ(τ )) , the imaginary time evolution of the Schrödinger equation on the trial state space is Therefore, we can evolve parameters to simulate the unphysical imaginary time evolution.
a. The Dirac and Frenkel variational principle. Applying the Dirac and Frenkel variational principle, we have and the evolution of parameters is simplified to with Similar to the real time evolution, even if the parameter θ j is real, the solution of its derivativeθ j may not be real as A and C are complex. b. McLachlan's variational principle Applying the McLachlan's variational principle [47] to imaginary time evolution, we have Whenθ i is complex, the evolution of parameters is the same as Eq. (37). Whenθ i can only take real values, the evolution becomes Due to the normalisation requirement for the trial state |φ( θ(τ )) , i.e., φ( θ(τ ))| φ( θ(τ )) = 1, we have ∂ φ( θ(t))| ∂θi E τ |φ( θ(t)) = 0 and therefore C R i = C R i . c. Time-dependent variational principles The Lagrangian for the Schrödinger equation is When |ψ = φ( θ(t)) has complex θ, we get the evolution of Eq. (37). Suppose θ is real, the evolution becomes where, A I i,j and C R i are the imaginary and real parts of A i,j and C i , respectively. The solution is imaginary and hence not physical.
When parameters are complex, the three variational principles are still equivalent. However, considering real parameters, only the McLachlan's variational principle can guarantee a real solution of the evolution equation.

B. Imaginary time evolution of mixed state
Now, we consider the case that the input state is a mixed state. Under the imaginary time evolution of Hamiltonian H, the state ρ(τ ) at imaginary time τ is where ρ(0) is the initial state. Equivalently, ρ(τ ) follows the time derivative equation where {H, ρ(τ )} = Hρ(τ ) + ρ(τ )H. Consider a parametrised state ρ( θ(t)), the evolution becomes a. The Dirac and Frenkel variational principle To get the evolution of parameters, we project the equation onto δρ( θ) according to the Dirac and Frenkel variational principle. Therefore, we have equal zero and the solution is with M defined in Eq. (27) and Y given by Note that when parameters θ are real, M and Y are also real, and we have In this case, the solution θ is also real note that this is opposite to the pure state case. b. MacLachlan's variational principle We can also minimise the distance between the evolution via parameters and the evolution, The evolution to˙ θ is the same as Eq. (47). Similarly, when˙ θ is real, Y i is reduced to Eq. (49) and we always have real solutions.

(51)
Suppose θ is complex, the evolution is the same as Eq. (47). Instead, when θ is real, we cannot obtain the equation for evaluating˙ θ.
As a special case, suppose ρ( θ(t)) is a pure state |φ( θ(t)) , we have the same M i,j as in real time evolution given in Eq. (33) and Y i = −C R i .

C. Variational Gibbs state preparation
Evolving the maximally mixed state I d /d with imaginary time τ can be used to prepare the thermal state ρ(T ) = e −H/T /Tr[e −H/T ] with temperature T = 1/2τ . Here, d is the dimension of the system. With variational simulation, we cannot simply apply a variational unitary circuit to the maximally mixed state and change parameters in the unitary circuit realise the simulation. This is because I d /d is invariant under unitary, U (θ) · I d /d · U (θ) † = I d /d. Instead, we can input a maximally entangled state |Φ d = 1/ √ d i |ii AE of system AE and evolve the whole system with Hamiltonian H ⊗I d under imaginary time τ . Then, the state of system A at time τ will be the thermal state with temperature T = 1/2τ .

VI. IMPLEMENTATION
In this section, we first show how to realise pure and mixed trial states with quantum circuits possibly assisted with post-selection on measurement outcomes. Then we show how to measure each term of M , V , and Y for simulating general evolutions of mixed states.
A. Trial state implementation a. Pure state. We first consider the pure state case. The conventional way of realising the trial state is to apply a sequence of parameterised gates to an initial state as |φ( θ) = R N (θ N ) . . . R k (θ k ) . . . R 1 (θ 1 ) |0 . In practice, not only the gate, but also the measurement performed to the state can be also parametrised, which however can always be effectively realised by applying a parameterised gate before a fixed measurement. Therefore, we also regard parametrised measurements equivalently as parametrised gates. Next, we consider state preparation via post-selection of measurement outcomes. Specifically, after initialising the target system together with ancillae, we can apply a joint parametrised circuit, measure the ancillae, and post-select the trial state conditioned on the measurement outcome. The quantum circuit of such a procedure is as follows.
As we only consider a pure trial state, the measurement should be a rank one projector, such as |0 0| E without loss of generality. Therefore, the trial state is with the joint state before the measurement |φ( θ) = R( θ) |0 A |0 E , and p( θ) = | φ( θ)| AE |0 E | 2 being the post-selection probability of measuring 0 of the ancillae. b. Mixed state and unitary evolution. For mixed states under unitary evolution, we can design the trial states via two different ways. We can either input a mixed state ρ(0) and apply a unitary circuit R( θ) to prepare the ansatz. Equivalently, we can input a state |0 AE of system AE such that Tr E [|0 AE 0 | AE ] = ρ A (0) with partial trace Tr E of system E, and only apply the circuit to system A such as the following one.

R( θ) |0 AE
In this case, we use one trial state ρ( θ(t)) = Tr E [R( θ) |0 AE 0 | AE R † ( θ)] with one parameter setting θ to directly represent the state ρ(t) at time t. Here, we can also consider more complicated ansatz with measurements such as Eq. (52). Alternatively, we can decompose ρ into pure states and simulate the unitary evolution of each pure state separately. Suppose the initial state ρ(0) has a decomposition ρ(0) = i p i |ψ i (0) ψ i (0)|. In practice, one can randomly input state |ψ i (0) with probability p i and evolve |ψ i (t) to time t under the Hamiltonian H. Then the state ρ(t) at time t would be ρ(t) = i p i |ψ i (t) ψ i (t)|. With VQS, |ψ i (t) is represented by the trial state |φ( θ i (t)) with parameters θ i (t). As the ansatz circuit R( θ) and parameters θ can be different for the evolution of each pure state, evolving the pure states separately could lead to more accurate results than evolving the whole mixed state with one parameter settings. However, as the evolution equation for one trial state is obtained by minimising the distance of the whole state, evolving the pure states separately may not always lead more accurate simulation. Furthermore, it can also be practically hard to decompose the initial state into a mixture of pure states. c. Stochastic and imaginary time evolution. For stochastic and imaginary time evolution evolutions, as they are generally not reversible, we have to introduce ancillae and jointly apply the parameterised circuit. We can also measure a subset of the ancillae to postselect the trial state. However, we can only use one trial state and it cannot be naively simulated by decomposing ρ(0) into pure states as ρ(0) = i p i |ψ i (0) ψ i (0)| and evolving each pure state separately. This is because stochastic and imaginary time evolution is not a linear process acting on the initial state so that ρ(T ) = i p i ρ i (T ) for real time T = t and imaginary time T = τ . Here ρ(T ) and ρ i (T ) are the evolved state of ρ(0) and |ψ i (0) at time T , respectively.
In Ref. [54], we do present an alternative method to simulate the stochastic master equation with different pure states evolving with different parameters. While the simulation involves several new techniques including variational simulation of general unphysical processes, we refer the readers to Ref. [54] for more information.

B. Coefficients evaluation
Now, we consider how to measure the coefficients M , V , and Y in the evolution equations. For unitary evolution of mixed states, we can effective regard the input state as a pure state and only evolve the subsystem of ρ 0 . Suppose the trial state is prepared by directly applying a unitary circuit to the initial state, the M , V , and Y can all be measured with the circuit in Fig. 1. Instead, when the trial state is prepared conditioned on postselecting a measurement outcome as defined in Eq. (52), we need to reconsider how to measure the coefficients. The derivative of the state can be written as The probability p( θ) can be directly measured and the derivative ∂p( θ) ∂θ k can be measured by the circuit in Fig. 1 To evaluate M , V , and Y , we need to measure the following terms for all i and j, (54) It is not hard to check that all these terms can be efficiently measured with the circuit in Fig. 1. Now, we consider the case of general evolutions. Suppose the trial state is prepared by ρ = R N R N −1 · · · R 1 (ρ 0 ), where ρ 0 is a joint state of the system and ancillae and R k are quantum gates that applied ancilla register-1 3. Circuit for the evaluation of differential-equation coefficients of the mixed-state simulator. To evaluate e iθ Tr(ρ1ρ2) , where ρ1 = RN RN−1 · · · R k+1 [S k (R k−1 · · · R2R1ρ0)T † k ], ρ2 = RN RN−1 · · · Rq+1[Sq(Rq−1 · · · R2R1ρ0)T † q ], the ancillary qubit is initialised in the state (|0 + e iθ |1 )/ √ 2 and measured in the |± = (|0 ± |1 )/ √ 2 basis. Here, U k = T k S † k and Uq = TqS † q . The last gate is a controlled swap gate on two registers. In the figure, we have assumed that k < q, and q = N + 1 when the circuit is used to evaluateṼ k coefficients. on the whole system. We express the generator as where S i and T i are unitary operators, and g i are coefficients. Similarly, we write where S k,i and T k,i are unitary operators, and r k,i are coefficients. For example, consider gate R k (θ k ) = e −iθ k σ/2 ρe iθ k σ/2 with a one or two qubit Hermitian operator σ, the derivative is ∂R k (ρ)/∂λ k = −iσ/2 · e −iθ k σ/2 ρe iθ k σ/2 + e −iθ k σ/2 ρe iθ k σ/2 · iσ/2. Using the expression (56), we rewrite the derivative of the mixed state as where Then, using the expression (55), coefficients can be expressed as In Eq. (59), each term is in the from a e iθ Tr(ρ 1 ρ 2 ) , where the amplitude a and phase θ are determined by either r k,i r q,j or r k,i g j . We would like to remark that, in general, ρ 1 and ρ 2 are not reduced density matrices as S and T are separately applied to the left and right sides of the state. Nevertheless, such a term can be evaluated with the quantum circuit shown in Fig. 3. This circuit needs an ancillary qubit initialised in the state (|0 + e iθ |1 )/ √ 2 and two registers initialised in the state ρ 0 . The ancillary qubit is measured after a sequence of operations on registers and some controlled unitary operations, where the ancillary qubit is the control qubit. The value is given by [e iθ Tr(ρ 1 ρ 2 )] = 2P + − 1, where P + is the probability that the qubit is in the state |+ = (|0 + |1 )/ √ 2. We can also consider trial state prepared by post-selection and the matrix elements can be similarly evaluated.

VII. DISCUSSION
In this work, we focus on the theory of variational quantum simulation. We first study the equivalence and difference of the three variational principles. We find McLachlan's variational principle is the most consistent principle for variational quantum simulation with quantum gates controlled by real parameters. Then, we extend variational quantum simulation of pure states to the general evolution of mixed states under both real and imaginary time. We discuss possible realisation of the trial states and show how to efficiently implement the simulation with quantum circuits.
In future works, one can study the design of trials states for specific problems and test our theory for simulating open quantum systems and preparing Gibbs states. It is also interesting to experimentally realise our variational simulation algorithms with current and near-term noisy quantum hardware. As the variational method only uses shallow quantum circuits, error mitigation techniques can be applied to suppress errors [39,[55][56][57][58][59][60][61].

ACKNOWLEDGEMENTS.
This work is supported by the EPSRC National Quantum Technology Hub in Networked Quantum Information Technology (EP/M013243/1). SE is supported by Japan Student Services Organization (JASSO) Student Exchange Support Program (Graduate Scholarship for Degree Seeking Students). QZ acknowledges support by the National Natural Science Foundation of China Grant No. 11674193. YL is supported by NSAF (Grant No. U1730449).
Appendix A: Real time evolution: pure state and unitary evolution

The Dirac and Frenkel variational principle
The Schrödinger equation is, Consider a parametrised trial state |φ( θ(t)) , with θ(t) = (θ 1 (t), θ 2 (t), . . . , θ N (t)), the real time evolution of the Schrödinger equation on the trial state space is Such an equation has a state vector on both sides. While, state ∂|φ( θ(t)) ∂θi can be regarded as the tangent vector of state |φ( θ(t)) at θ(t), we can thus apply a projector to project the right hand side onto the tangent space. Therefore, the equation becomes, For each term of ∂|φ( θ(t)) ∂θi , the coefficient should satisfy Define the matrix elements of M and V as the evolution of parameters is simplified to The same equation can be obtained by applying the Dirac and Frenkel variational principle This is can be verified with and Eq. (A2).
The McLachlan's variational principle [47], applied to real time evolution, is given by Supposeθ i can be complex, then we have Then the evolution is On the other hand, supposeθ i is real, then The corresponding evolution equation for parameters is where A R i,j and C I i are the real and imaginary parts of A i,j and C i , respectively.

Time-dependent variational principles
The Lagrangian for the Schrödinger equation is and the Schrödinger equation is obtained by Suppose |ψ = |φ( θ(t)) , then the Lagrangian is and the evolution is Instead, suppose θ is real, then we have Equivalently, the evolution to parameters is where, A I i,j and C R i are the imaginary and real parts of A i,j and C i , respectively.
Appendix B: Real time evolution: mixed states and general evolution

General evolution
Suppose the initial state is a mixed state ρ, the stochastic evolution is defined by Consider a parameterised trial state ρ( θ(t)), we can get the evolution of θ that simulates the stochastic evolution.
a. The time-dependent variational principle. Let ρ = ρ( θ(t)), the stochastic evolution equation is ∂θiθ i is in the tangent subspace of { ∂ρ( θ(t)) ∂θi }, we need to project the stochastic evolution equation onto this tangent subspace. Following the Dirac and Frenkel variational principle, we have The evolution of parameters is with Note that when parameters θ are real, M and V are also real and the solutionθ j is also real. b. McLachlan's variational principle. We can also minimise the distance between the evolution via parameters and the stochastic evolution via McLachlan's variational principle, where ∂ρ/∂t − L(ρ) =Tr (∂ρ/∂t − L(ρ)) † (∂ρ/∂t − L(ρ)) , Supposeθ is complex, then The evolution toθ is Whenθ is real, then The evolution toθ is Then we can get the same equation for parameters c. The time-dependent variational principle. The Lagrangian for the von Neumann equation is the Schrödinger equation is obtained by Suppose ρ = ρ( θ(t)), then the Lagrangian is Suppose θ is complex, then and the evolution is Instead, suppose θ is real, then we have Then does not lead to a meaningful evolution of parameters.
Proof. Suppose there are N parameters in the ansatz. Consider an additional ancilla that is rotated via a phase gate e iθ N +1 , the evolution of the N + 1 parameters via the pure state equation is equivalent to the evolution of the N parameters via the mixed state equation. Under the pure state equation, the evolution of the N + 1 parameters is Substituting the last equation into the first N equations, we get which it is equivalent to ∂ φ( θ(t))| ∂θ j |φ( θ(t)) θ j = C I i + i ∂ φ( θ(t))| ∂θ i |φ( θ(t)) φ( θ(t))| H |φ( θ(t)) .
This is can be verified with and Eq. (C2).
Proof. Suppose there are N parameters in the ansatz. Consider an additional ancilla that is rotated via a phase gate e iθ N +1 , the evolution of the N + 1 parameters via the pure state equation is equivalent to the evolution of the N parameters via the mixed state equation. Under the pure state equation, the evolution of the N + 1 parameters is N j=1 ∂ φ( θ(τ ))| ∂θ i |φ( θ(τ )) θ j −θ N +1 = 0. (D28) Substituting the last equation into the first N equations, we get N j=1 which it is equivalent to