Non-unitary Trotter circuits for imaginary time evolution

We propose an imaginary time equivalent of the well-established Pauli gadget primitive for Trotter-decomposed real time evolution, using mid-circuit measurements on a single ancilla qubit. Imaginary time evolution (ITE) is widely used for obtaining the ground state of a system on classical hardware, computing thermal averages, and as a component of quantum algorithms that perform non-unitary evolution. Near-term implementations on quantum hardware rely on heuristics, compromising their accuracy. As a result, there is growing interest in the development of more natively quantum algorithms. Since it is not possible to implement a non-unitary gate deterministically, we resort to the implementation of probabilistic imaginary time evolution (PITE) algorithms, which rely on a unitary quantum circuit to simulate a block encoding of the ITE operator - that is, they rely on successful ancillary measurements to evolve the system non-unitarily. Compared with previous PITE proposals, the suggested block encoding in this paper results in shorter circuits and is simpler to implement, requiring only a slight modification of the Pauli gadget primitive. This scheme was tested on the transverse Ising model and the fermionic Hubbard model and is demonstrated to converge to the ground state of the system.


I. INTRODUCTION
Imaginary time evolution (ITE) is a routine that is widely used on classical hardware to obtain the ground state (GS) of a system.In this procedure, an initial trial state is evolved through imaginary time, τ = it, by applying the ITE operator, e − Ĥτ .This operator has the property of being a ground state projector in the large τ limit; over time, the excited state contributions from the initial state |Ψ(τ = 0)⟩ decay until only the GS, |ϕ 0 ⟩, remains lim τ →∞ e −( Ĥ−E0)τ |Ψ(τ = 0)⟩ ∝ |ϕ 0 ⟩ . ( Although it is not expected that quantum computers will be able to make the task of ground state determination for generic Hamiltonians efficient [1], they may provide useful speedups in the near future for certain applications using quantum dynamics [2].In particular, owing to the exponential scaling of the Hilbert space with system size, classical techniques for quantum simulation are generally designed to avoid two features in a computation: (1) explicitly storing the many-body wave function, and (2) propagating the wave function by matrix exponentiation.Quantum computers allow for efficient implementation of both of these features.Regarding the former, there are quantum circuits which cannot be efficiently simulated using classical resources [3].The quantum computer thus has access to a richer state space for the ground state search of certain systems.Regarding the latter feature, there are quantum algorithms for efficiently simulating real time evolution.Each algorithm corresponds to a different decomposition of the real time evolution operator; for instance, the operator can be expressed as a truncated Taylor series, which can be implemented via a linear combination of unitaries (LCU) [4][5][6][7], or it can be expressed as a series of alternating signal rotations and signal processing rotations via quantum signal processing (QSP) [8].Due to its simplicity and ease of implementation with the Pauli gadget primitive (refer to Appendix 1) [9], the most widely adopted method is the Trotter decomposition (alternatively referred to as the 'product formula' approach) [10][11][12][13], which evolves each Hamiltonian term for a small time slice.
Imaginary time evolution is much more challenging to implement, owing to the non-unitarity of the ITE operator: since quantum gates are unitary operations, the ITE operator cannot be implemented deterministically on quantum hardware.Many of the previous approaches to this problem have focussed on implementing the ITE operator indirectly using a hybrid quantum-classical approach.These algorithms are thus designed for hardware in the current 'noisy arXiv:2304.07917v3[quant-ph] 29 Oct 2023 intermediate-scale quantum' (NISQ) era, characterised by a limited coherence time and gate fidelity.Whilst the partitioning between classical and quantum hardware saves quantum resources, it typically necessitates the use of heuristics, which can hinder the attainable accuracy of the final ground state.Indeed, the variational quantum eigensolver [14,15], which is among the most promising and widely-used examples of variational algorithms, requires gate errors significantly lower than the error thresholds of fault-tolerant hardware in order to achieve chemical accuracy [16].Within the class of hybrid ITE techniques, two main approaches have emerged: variational ITE (VITE) [17][18][19][20][21] and quantum ITE (QITE) [22][23][24][25].Refer to Section II A for more information regarding these methods.
Unlike their near-term competitors, fault-tolerant algorithms do not rely on heuristics and offer rigorous performance guarantees.To this end, we consider a third approach: the probabilistic ITE (PITE).PITE algorithms are a type of block encoding -that is, the non-unitary operator is embedded inside a larger unitary matrix, which is accessed via post-selection of 'successful' mid-circuit measurements.
Unlike VITE, PITE does not require an Ansatz, and unlike QITE, it does not require the selection of Pauli strings according to the locality of the Hamiltonian.PITE evolves the system exclusively using quantum operations, without the need for classical updates to the state.The main drawback of PITE algorithms arises from the need to implement them as a series of iterative circuits [26], each of which must be applied successfully -that is, each imaginary time step is implemented using at least one successful mid-circuit measurement.
Consequently, as with most block encoding methods, PITE algorithms are characterised by an exponentially decaying probability of success; though they are guaranteed to converge to the GS, convergence is limited by the number of measurement shots required.The practicality of implementing most block encoding methods, including PITE, will be tied to the development of amplitude amplification (AA) schemes to boost the overall success probability of the circuits [27][28][29].AA has been utilised to boost PITE circuits with some success [30].
The ITE operator can be decomposed using similar strategies as the RTE operator; namely, QSVT [31,32], LCU and Trotterisation.PITE algorithms use midcircuit measurements, which naturally renormalise the quantum state throughout the evolution process; this represents an advantage of PITE methods over other approaches which must compute this normalisation, including hybrid methods such as QITE, as well as QSVT.In recent years, a number of PITE algorithms have been proposed [26,[33][34][35] and implemented on hardware [36,37].Amongst these, the LCU proposal from Kosugi et al [35] is the most comparable to the scheme that will be proposed in this paper: both methods are general schemes for constructing the necessary circuits directly given any input Hamiltonian, without the need for system-specific block encodings and without incurring an additional classical cost in computing an initial singular value decomposition of the ITE operator [33].In the LCU proposal, a firstorder approximation of the ITE operator, (1 − Ĥ τ r ) r , where r is the number of time steps, is achieved by the repeated application of controlled forward and backward RTE operators, requiring two controlled Pauli gadgets for each Hamiltonian term per time step.
In this paper, we propose an alternative block encoding of the ITE operator, assuming a Trotterdecomposed form.The algorithm presented can be considered to be the imaginary time counterpart of the well-established Trotterised RTE algorithm, which makes use of the Pauli gadget primitive.There are several advantages to using this approach: firstly, the structure of the Pauli gadget naturally results in diagonal matrices, avoiding the need to incur an initial classical cost in computing the singular value decomposition (SVD) of the ITE operator [33].Secondly, the similarity of the proposed Trotterised ITE algorithm to the Trotterised RTE algorithm also presents advantages.Trotterisation is an active area of research, owing to its simplicity and unique scaling with the commutative structure of the Hamiltonian.In RTE simulations, this can be leveraged to achieve much better empirical performance than more sophisticated algorithms, including LCU-Taylor and QSP [38].By building on the Pauli gadget framework, we ensure that any of the optimisation passes previously designed for RTE Pauli gadgets can be directly transferred to our circuits [39].Further to this, the Trotter decomposition is highly flexible, with many flavours extensively reviewed by the community; as such, there are many options readily available for improving the accuracy of the evolution, which is of particular importance for the task of thermal state preparation -on the other hand, the LCU PITE proposal is only a first-order approximation of the ITE operator.Finally, it is easy to account for non-Hermitian Hamiltonians by concatenating the proposed (probabilistic) imaginary time Pauli gadgets with the ubiquitous (deterministic) real time Pauli gadgets.
The structure of the paper is as follows.Background information is given in Section II, describing the implementation of non-unitary operators on quantum circuits.The proposed algorithm will then be presented in Section III, and the results obtained will be detailed in Section IV.Finally, a discussion of the results will be presented in Section V, and conclusions drawn in Section VI.

A. NISQ ITE Methods
Heuristic ITE approaches approximate the action of the non-unitary ITE propagator with a unitary operator.
VITE represents the unitary operator with a parameterised quantum circuit, whose parameters are optimised classically according to a variational principle [18] which carries out a finite projection of the wavefunction on the tangent space of the exact evolution.The standard VITE optimisation procedure requires the inversion of a matrix whose elements are determined by measurements made on the quantum circuit; due to the sensitivity of matrix inversion to external environmental pertubations, this procedure assumes very low gate errors [40].Matrix inversion can be avoided, and thus the hardware requirements of the standard VITE procedure can be significantly reduced, by assuming a Trotter product form for the Ansatz and variationally optimising each Trotter term sequentially [21].Further, optimising the circuits using classical tensor network methods can result in shallower and more accurate circuits [41].Variational quantum algorithms (VQAs) such as these have substantial drawbacks [42], including the barren plateau (BP) problem, which could prevent convergence for larger problem sizes [43].There are multiple driving factors; in particular, the more expressive the Ansatz is, the more susceptible it is to BPs [43,44].A common approach to mitigate the BP problem is to reduce the expressibility of the Ansatz.In the context of quantum chemistry problems, this corresponds to restricting the span of the Ansatz to a section of the Hilbert space, for instance by exploiting the symmetry of the problem [44].Consequently, VQAs are highly dependent on the choice of Ansatz.
On the other hand, QITE represents the unitary operator in a Pauli basis and solves a linear system of equations for the coefficients of the Pauli tensors; measurements taken from the quantum circuit are used to estimate the matrix elements, and the equations are solved on classical hardware.The derivation is mathematically equivalent to McLachlan's variational principle for VITE [35].However, the domain of these Pauli tensors grows exponentially with the spreading of entanglement.As a result, the matrix dimension for the linear equations, and correspondingly the number of measurements, tends to grow rapidly as the problem size increases.Its scaling can be regularised if restrictive approximations regarding the locality of Hamiltonian terms are adopted [25].

B. Block encoding scheme for non-unitary operation
Block encoding represents a general framework for the implementation of non-unitary operators: a nonunitary matrix A that operates on an n-qubit quantum state |ψ⟩ is embedded as a block inside a larger unitary matrix U A that operates on (n + l) qubits.In theory, any non-unitary matrix A can be block encoded with a single ancilla using a singular value decomposition scheme [45].However, these circuits may be exponentially deep.Application of A is then contingent on successfully post-selecting the l ancillary qubits corresponding to this block, |0⟩ ⊗l ⟨0| ⊗l ⊗ |ψ⟩ ⟨ψ| of the unitary matrix U A .For the purpose of this work, we will iteratively use a single ancilla qubit, i.e. l = 1, block encoding.When l = 1, the block encoded matrix is given by where 1 α ∈ R is a necessary scaling required to ensure that U A is unitary; more precisely, the spectral norm of The spectral norm ∥ • ∥ of a matrix A is the maximum scale by which A can stretch a vector, defined as where σ max (A) is the maximum singular value of A.
Entries given by an asterisk * indicate that there is some freedom in choosing their values -again, the only requirement is that U A is unitary.The action of U A is: The first and second terms represent the components of the state |Ψ⟩ that are parallel and perpendicular to |0⟩ in the auxilliary space, respectively.That is, |⊥⟩ must satisfy the following orthonormality constraints: (⟨0| ⊗ I ⊗n ) |⊥⟩ = 0 and ⟨⊥|⊥⟩ = 1 − 1 α 2 ⟨0|0⟩ ⟨ψ| A † A |ψ⟩, but its system space component is determined by the action of the bottom left element in U A (2) on |ψ⟩.To successfully apply A, one must measure the ancilla qubit to be in the state |0⟩.
According to the postulates of quantum mechanics, following a partial measurement on the ancillary register, the remaining state vector is renormalised, preserving its purity.Thus, the value of 1 α does not affect the successful post-measurement state.To see this, consider the projective measurement operator M0 = (|0⟩ ⟨0| ⊗ I ⊗n ), which only acts on the ancilla qubit, applied to the state |Ψ⟩.The successful ancillary measurement outcome is defined to be The value of 1 α does, however, affect the probability of successfully applying A, The dependence of the success probability on the initial state |ψ⟩ is a consequence of the non-unitarity of 1 α A -roughly speaking, the more non-unitary 1 α A is, the more information that is lost from |ψ⟩ on application of 1  α A, and the smaller the success probability is.Of course, the non-unitarity of 1 α A is affected by the value of the scaling, 1 α : the greater the scaling, the greater the success probability.Different block encodings U A of A may produce different scalings; however, there is a maximum value this scaling can take, 1 α max .Above this value, the block encoding U A is no longer unitary and cannot be implemented with a quantum circuit.This optimal block encoding -that is, the block encoding that applies A with the maximal success probability for any |ψ⟩ -satisfies ∥ 1 α max A∥ = 1, giving When the scaling of a block encoding is not yet optimal, a classically-controlled reversing measurement scheme may be applied to improve the success probability [46].In this scheme, an unsuccessful measurementthat is, the application of the projective measurement operator M1 = |1⟩ ⟨1| ⊗ I ⊗n on the state |Ψ⟩ -triggers the application of a reversing non-unitary operator, R0 , defined such that R0 M1 |Ψ⟩ ∝ |Ψ⟩.The more times the reversing non-unitary is applied, the more the success probability is improved; in the limit of an infinite number of applications of R0 , the optimal success probability is reached, The intuition behind the block encoding approach can be seen by drawing an analogy with the physical process of system-environment interactions.In physical systems, non-Hermicity is generated from entanglement with the environment; in Equation ( 5), we can identify the 'system' to be the n-qubit state |ψ⟩, whereas the 'environment' is represented by the ancillary qubit.First, U A entangles the system with the environment.Next, the local measurement on the ancilla state, before checking the measurement result, produces a mixed state and corresponds to the process of taking the partial trace over the environmental degrees of freedom, resulting in a non-unitary evolution of the system.In particular, however, we are interested in the non-unitary evolution of the component of the system state vector that is entangled with the |0⟩ ancilla state, |//⟩.Initially, we begin in the normalised state |//⟩ = |0⟩ ⊗ |ψ⟩.Application of U A is a continuous real time evolution process, throughout which the magnitude of |//⟩ decreases.At the end of this process, we have |//⟩ = |0⟩ ⊗ 1 α A |ψ⟩ and the magnitude of |//⟩ has decreased by a factor of √ p success (Equations ( 5) and ( 7)).Finally, after looking at the measurement result, the system collapses into the corresponding pure state.This is the second process: discontinuous renormalisation of the state vector (Equation ( 6)).

C. Universal gate set for non-unitary operation with a single ancillary qubit
The cnot gate and all single qubit unitary gates constitute a universal set for the unitary quantum circuit.Terashima et al. [46] proved that, considering only block encodings with one ancillary qubit, l = 1, a possible universal set for non-unitary circuits is given by the cnot gate, all single qubit unitary gates, and a specific single qubit non-unitary gate, N 1 (a), which relies on a single qubit projective measurement.N 1 (a) has the matrix representation for some a ∈ [0, 1].Its circuit construction is given in Figure 1; it is composed of a controlled single qubit unitary gate U 1 (a), whose matrix representation is along with a single ancillary measurement.Again, entries given by an asterisk * indicate that there is some freedom in choosing their values; specifically, a single qubit unitary is described by three degrees of freedom.Note that |a| ≤ 1 is a necessary requirement for U 1 to be unitary.Written explicitly using little endian ordering -that is, in the |ϕ anc ⟩ ⊗ |ψ⟩ basis, where |ϕ anc ⟩ denotes the ancillary qubit state -the unitary corresponding to the circuit in Figure 1, U N1 (a), is a block encoding of N 1 : Upon post-selecting the 0 measurement on the ancilla qubit one obtains the N 1 (a) matrix of elements shown in red: Using this decomposition, successful implementation of a non-unitary circuit corresponds to a measurement shot in which all ancillary measurements from the constituent N 1 (a) gates yield |0⟩.
Each constituent N 1 (a) gate is implemented as a block encoding (12) according to Figure 1, and the overall block encoding we achieve for the ITE operator will determine its scaling 1 α (Section II B).
This is a natural approach to take when the ITE operator has been Trotterised.Given a Hamiltonian expressed in a qubit (Pauli) basis, the Trotter decomposition approximately factors the ITE operator into a product of exponentiated Pauli strings.Since each of these terms can be easily diagonalised with the use of single qubit rotations (refer to Section III A), they can be readily represented in terms of the diagonal N 1 (a) matrices.
The proof presented by Terashima et al [46] could be regarded as a general synthesis strategy for non-unitary operators; however, no attempt has been made to optimise the number of gates for the specific non-unitary operation of ITE.Moreover, the proof relies on an initial singular value decomposition (SVD) of the nonunitary matrix in question, incurring a classical cost that scales exponentially with the number of qubits n, O 2 3n .Instead, we consider a decomposition based on the Pauli gadget framework for Trotterised real time evolution (RTE) [9].This allows us to transfer any optimisation passes previously designed for real time propagation with Pauli gadgets directly to our framework [39].

A. Trotter decomposition
Consider an n qubit Hamiltonian that can be expressed as a sum of M local interactions, Ĥ = M m=1 Ĥm , where M ∼ O (poly(n)) and each Hamiltonian term is Hermitian.Applied to an imaginary time evolution τ = it, product formulae approximate the evolution under the full Hamiltonian, e − Ĥ∆τ , as a sequence of operators, e − Ĥm∆τ , which should each be efficiently implementable on a quantum computer.For instance, the first-order Trotter formula states that [10] where L is the system size.For a longer simulation time, the evolution is divided into r ∈ N Trotter steps where τ r = ∆τ and r is also referred to as the 'Trotter number'.From Equation (15), we can determine the gate count N required for a given evolution time τ and target precision ϵ: The main challenge is to choose r to be as small as possible and still ensure a total simulation error of at most ϵ.This choice is complicated by the fact that the upper bound on the Trotter error in Equation ( 16) can be rather loose.In particular, the error is strongly dependent on the commutator structure of the Hamiltonian; for instance, in the limiting case, all M Hamiltonian terms commute and the error is exactly zero.Consideration of the commutator structure can be used to reduce the theoretical error scaling with respect to L [47]; however, simulations can remain orders of magnitude faster than theory even when many Hamiltonian terms commute [10,38,48].
Grouping terms into mutually commuting partitions can reduce the circuit depth significantly: not only is the Trotter error reduced, allowing for fewer Trotter steps to be taken, but the ordering of terms within a mutually commuting partition can be optimisedwith no effect on the Trotter error -to maximise gate cancellations [39].Partitioning the Hamiltonian into the minimum number of mutually commuting groups is equivalent to the minimum clique cover problem, which is NP-hard [49]; thus, it is typically performed heuristically [50,51].
The product formula approach is highly flexible.For instance, higher order product formulae can be defined [52]; the higher the order, the more accurate the Trotterisation.This presents a trade-off between the order of the decomposition and r: increasing the order reduces the number of Trotter steps r required to achieve a fixed ϵ, whilst also increasing the cost per Trotter step.Recent approaches to Trotter error mitigation have considered randomly permuting the Hamiltonian terms [53,54].Further, classical stochastic ITE (quantum Monte Carlo) methods randomly sample Hamiltonian terms, effectively weighted by the population of walkers located between connected basis states.They rely on the outcome that, in so doing, it is still possible to converge to the ground state [55][56][57].Similarly, the qDRIFT algorithm randomly samples Hamiltonian terms weighted by their coefficients: despite largely forsaking knowledge of the internal commutative structure of the Hamiltonian, it eliminates the explicit dependence of the gate count on the system size (Equation ( 16)) [58][59][60][61].In this paper, we only implement the first-order Trotter formula of Equation ( 14) in our experimental simulations and theoretical examination, however, our techniques can be applied to any of the other decompositions.
Given a fermionic Hamiltonian in second quantised form, the Hamiltonian must first be converted to a qubit representation.That is, fermionic creation and annihilation operators are mapped onto qubit operators (namely, tensor products of Pauli operators, or 'Pauli strings').To this end, we will use the Jordan-Wigner (JW) transformation [62], which maps electronic configurations onto computational basis states; for instance, the state |01⟩ indicates a system with two spin orbitals, in which the first (second) spin orbital is unoccupied (occupied).Thus, for the remainder of this paper, we will consider Hamiltonians expressed in a Pauli basis, where σm is a Pauli string of length n,

B. Pauli gadgets for Trotterised PITE
We now approach the task of block encoding the nonunitary operator, 1 α e − Ĥτ for Ĥ = Ĥ † , via a Trotter decomposition.Each Hamiltonian term Ĥm = Ĥ † m in the Trotter decomposition can be diagonalised by a unitary operator Ûm , yielding: In particular, Pauli strings are unitarily mapped onto other Pauli strings via Clifford operations.The Clifford group can be generated by three gates: the Hadamard gate (H), the phase gate (S), and the cnot gate.We could choose to diagonalise each Pauli string individually: where Λm k = {Z, I} ⊗n and Bm k denotes the basis transformation operator from the σm k -basis to the computational basis: BX = H, BY = SH and BI = BZ = I.Alternatively, we could simultaneously diagonalise a mutually commuting subset of Pauli strings; in this case, the diagonalisation can be constructed efficiently from Clifford operations in a classical preprocessing step [63,64].This approach can reduce the circuit depth significantly when large numbers of terms form a mutually commuting subset, since it allows for more cnot gate cancellations.It is often also applied as a standard procedure for reducing the number of shots required for the evaluation of expectation values [50,51].
Pauli grouping methods are applicable to general qubit Hamiltonians.However, methods that produce a compressed representation of the original Hamiltonian typically yield better performance when applicable [51]; until recently, the leading method for reducing the cost of Hamiltonian simulation for quantum chemistry Hamiltonians was tensor hyper contraction, which uses non-unitary rotations and therefore restricts the method to LCU approaches [65].There are now superior double factorisation techniques which use unitary rotations [66] and can be combined with Trotter decompositions and its non-unitary variants, including the algorithm presented in this work.
For simplicity, we will adopt the first approach (Equation ( 20)) in this work.
Once this change of basis has occurred, the exponentiated diagonal Pauli string must be implemented.Using the ubiquitous phase gadget circuit structure (refer to Appendix 1), this task requires us to find a block encoding for the single-qubit diagonal nonunitary operator e +|γ|∆τ Ẑ , for any γ ∈ R. We note that a circuit structure implementation of the operator e +|γ|∆τ ( Ẑ− Î) has recently been proposed [34]; our work can be seen as an extension of this to include any Hamiltonian expressed in a Pauli basis.
As discussed in Section II B, 1  α is a scaling factor that depends on the block encoding used.The greater this scaling is, the greater the success probability of the block encoding; its maximum allowable value, which still ensures that the block encoding is unitary and thus implementable on a quantum circuit, is discussed in the following section (Section III B 1).

Maximum scaling for the block encoded ITE operator,
First, we consider the exact (non-Trotterised) ITE operator.According to Equation ( 8), the maximum possible scaling for the block-encoded exact ITE operator is given by ) where E 0 is the (unknown) ground state energy.Once the ground state is reached, the block encoded ITE operator cannot project out any more states: continued application does not change the state, but if the scaling is less than the maximal value, the success probability will be less than 1, p success = 1 α 2 e −2E0 .On the other hand, if the scaling is maximal, the ITE operator will be applied deterministically, p success = 1.In other words, this maximal scaling ensures that the ITE operator becomes the identity operator once the ground state is reached.
Note that this scaling is also present in our original definition of the ITE procedure for obtaining the ground state of a system (1).When using ITE methods which do not iteratively renormalise the wavefunction, the prefactor in this definition ensures that the amplitude of the wavefunction does not decay to zero as the ground state is approached.Typically this will require re-scaling of the operator; consider, for instance, classical stochastic ITE methods (quantum Monte Carlo), in which the wavefunction is propagated in imaginary time by sampling the action of the Hamiltonian in some discrete basis populated by 'walkers'.Since the total number of walkers at any one point in the evolution is related to the normalisation of the wavefunction, the ground state energy is estimated Ẽ0 ∼ E 0 throughout the run-time of the algorithm and used to shift the Hamiltonian, e −( Ĥ− Ẽ0)τ , as a means of walker population control [55][56][57]67].On the other hand, PITE algorithms naturally renormalise the intermediate state |ψ i ⟩ at each time step i through the use of partial measurements: rather than affecting the normalisation of the quantum state, due to the cancellation shown in Equation 6.The scaling 1 α affects only the success probability of these measurements (Section II B).Thus, the factor of e E0τ in Equation ( 1) is not required by PITE to prevent the decay of the ground state: 1 α < e E0τ is allowed.Instead, this factor manifests as the optimal scaling 1 α max that minimises the decay of the success probability.We may want to take inspiration from quantum Monte Carlo methods and vary the scaling of the block encoding.However, since this is a non-unitary operation, it is more difficult to implement it using quantum processes.
It is likely that PITE algorithms will need to be run in conjunction with an AA procedure to boost the success probability of the circuits.The PITE block encoding proposed by Liu et al [33] was defined in terms of an extra parameter; they were able to vary this parameter throughout the simulation, alongside Grover's algorithm [27,28], to provide more flexible optimisation of the number of Grover iterates required.In a similar manner, Nishi et al [30] were able to make 1 α a tunable parameter within the LCU-based PITE framework [35], with reduced computational cost as compared with the fixed point search [68,69] or oblivious AA routine [6].As will become apparent, the value of 1 α for the block encoding we propose in Section III B 3 is initially fixed by the circuit implementation we use.However, we note that our framework is also amenable to varying the value of the scaling throughout the simulation, although as discussed, this is not a trivial task and does not lie within the scope of this work.
Whilst it is easy to determine the maximum scaling for the exact ITE operator, it is not obvious how to relate this to the Trotterised form, which would depend on the Trotter error [47].Weak upper and lower bounds for the maximum scaling of the Trotter decomposed operator, defined by Equations ( 15) and ( 17), can be obtained (Appendix 2): where λ = we find that the maximum scaling for the exact ITE operator also lies in this range, 2. Block encoding for e −γ∆τ Ẑ with a single ancillary qubit As a preliminary to block encoding the exponentiated Hamiltonian, we first consider the task of block encoding the simple operator 1 α e −γ∆τ Ẑ , where γ ∈ R. According to Equation (8), this operator is implemented with the maximum possible success probability when the rescaling is given by: The matrix representation of the optimally block encoded operator 1 α max e −γ∆τ Ẑ can be easily expressed in terms of the N 1 (a) matrix structure (Section II C).Consider first the case γ < 0: Similarly, for the case γ > 0, we identify: (28) In order to implement Equations ( 27) and ( 28) using the circuit structure defined in Figure 1, we must identify a single qubit unitary gate U 1 (e −2|γ|∆τ ).Setting ϕ = 2 arccos e −2|γ|∆τ , we find that we can easily implement U 1 (e −2|γ|∆τ ) as a rotation, R x (ϕ).Overall, the implementation of 1 α max Â = e −|γ|∆τ e −γ∆τ Ẑ is obtained with the circuits shown in Figures 2 and 3, provided that the ancillary state is measured to be in the |0⟩ state.where ϕ = 2 arccos e −2|γ|∆τ .The post-selected state is renormalised by the square root of the success probability, psuccess, denoted here ps = 1 + e −4|γ|∆τ .

Block encoding for Trotterised e − Ĥτ
Since the |0⟩ post-selected two qubit gate becomes a diagonal non-unitary single qubit gate, we can use the same machinery of the N -qubit parity gate, otherwise known as the phase gadget (Figure 9 in Appendix 1).Combining the phase gadget structure with the non-unitary single qubit gate (Figures 2 and  3) gives the overall block encoding for the operator e −|γ|∆τ e −γ∆τ Ẑ⊗n (Figure 4).An example circuit for the ITE of the Pauli string X Ŷ Ẑ is given in Figure 5; only a small modification needs to be made compared with the widely-used RTE Pauli gadget (Figure 10).The overall circuit is then formed by concatenating circuits corresponding to each of the exponentiated Pauli strings, with their corresponding basis transformation gates (A.2), and for each of the Trotter steps, according to the Trotter decomposition defined in Equations ( 15) and (17) maximised (Section III B 2): From this, the overall scaling for the Trotterised PITE block encoding is given by where the sum in the second line (31) is carried out over all coefficients of the constituent Pauli strings in Ĥ, excluding any Pauli tensors of identity operators Î⊗n , and λ = M k=1 |c k | is the 1-norm of the Hamiltonian.Equality is achieved when the Hamiltonian does not include the term Î⊗n .Although each individual Trotter term has been optimally block encoded, the scaling for the overall Trotter product will depend on the order of the Trotter terms and may not be optimal (Equation ( 23)).
Successful application of the Trotterised ITE operator requires all mid-circuit ancillary measurements to yield |0⟩ -if a measurement is unsuccessful, the state is projected into the wrong subspace.Given a Trotter decomposition defined according to Equations ( 15) and ( 17), the success probability after r Trotter steps is given by As expected, the success probability decays exponentially with the number of time steps taken and depends heavily on |ψ⟩.The lowest success probability occurs for the state |ψ⟩ worst which undergoes pure amplitude damping under the action of every exponentiated operator in the Trotter decomposition -that is, e −σm k c k ∆τ |ψ⟩ worst = e −|c k |∆τ |ψ⟩ worst for all k, giving rise to the final inequality in Equation (33).This would require |ψ⟩ to be a simultaneous eigenstate of all the Pauli strings comprising Ĥ, which is only possible when all the Pauli strings comprising Ĥ commute; however, consideration of this scenario does provide us with a lower bound for the general case.
It is evident that the success probability for this algorithm can be increased by reducing the 1-norm of the Hamiltonian (33); this task is also of central importance to the reduction of gate complexity for several simulation algorithms, including simulation based on an LCU decomposition [5].Loaiza et al [70] have carried out significant work towards reducing the 1-norm of a molecular Hamiltonian; for instance, they describe several approaches for changing the Pauli basis decomposition of the Hamiltonian (17), including grouping together anti-commuting Pauli strings.When applied to the Trotterised PITE procedure, we must consider that whilst using these methods to reduce the 1-norm of the Hamiltonian would result in a higher success probability, it may also increase the Trotter error, resulting in a higher gate complexity.An alternative approach for reducing the 1-norm is to apply a symmetry shift operator to the Hamiltonian, e − Ĥτ = e − Ŝτ e −( Ĥ− Ŝ)τ , such that the 1-norm of the operator Ĥ − Ŝ is less than that of Ĥ. Provided the state vector satisfies symmetry constraints, for real time evolution t = −iτ , the unitary operator e −i Ŝt applies a phase shift to the state vector and can be ignored.However, this cannot be exploited for imaginary time evolution: shifting the Hamiltonian leaves behind a non-unitary operator, e − Ŝτ , which, as previously discussed in Section III B 1, is non-trivial to implement as a quantum circuit.
In practice, currently one must run a number of circuits to the end of the simulation, and the shots for which all ancillary measurements yielded |0⟩ -that is, the successful shots -are post-selected.Since the probability of successful application decreases exponentially with the number of mid-circuit measurements, the number of successful shots from which to estimate the expectation value also decreases exponentially; a sufficiently high number of shots must be selected to ensure that the final energy is determined with high accuracy.In future applications, the implementation of a quit-if-fail functionality would significantly improve the overall run-time of the algorithm.
The gate complexity and the total number of shots required to obtain an accurate estimate of the ground state energy at the end of the ITE procedure depend on the total evolution time τ for convergence.Consider expressing an arbitrary state |ψ(0)⟩ in the eigenbasis of Ĥ, |ψ(0)⟩ = [71].The unnormalised state after a time τ is and its normalisation is given by To specify how long we need to evolve the state, let τ be the time required to achieve an accuracy of ϵ ≪ 1 in the squared overlap with the ground state, In the large τ limit, only the ground state and first excited state contributions are left: The number of times, r ϵ = τ ∆τ , that a time step of ∆τ must be applied to achieve an error, ϵ, is thus given by As expected, the more easily distinguishable the ground state and first excited state energies are, the faster convergence to the ground state is achieved.Further, the better the initial guess, i.e. the greater the initial overlap with the ground state, the faster convergence is achieved.
Combining Equations ( 33) and ( 39), assuming ∆τ is small enough to be able to neglect the Trotter error, would in principle give a lower bound on the total number of shots required to obtain an error of ϵ in the ground state estimation.
The circuit depth required for convergence to the GS is then determined by the number of Pauli gadgets, rM , where M is the number of local terms comprising the Hamiltonian (17).Each of these Pauli gadgets requires 2n cnot gates; thus each Pauli gadget produces a gate depth that scales at worst linearly in n, and at best logarithmically in n [39].The whole procedure can be implemented using one ancillary qubit by resetting it to |0⟩ following every mid-circuit measurement.

IV. RESULTS
Pytket [72] is used for the construction and compilation of the circuits, and all quantum simulations are performed with the Qiskit Aer simulator [73].For each of the Hamiltonians, the constituent Pauli strings are partitioned into mutually commuting sets to reduce the Trotter error.
Energies and their associated errors are estimated in the following manner.Individual Pauli strings are sampled and the mean energy is constructed as The error ∆E in the mean energy arising from using a finite number of shots is determined using a similar method to Kandala et al [74]: where ⟨∆σ 2 m k ⟩ is the variance on Pauli string σm k and N k is the number of successful shots used in the measurement of ⟨∆σ 2 m k ⟩.
A. Models studied

Ising Model
The transverse field Ising model (TIM) [75] is the simplest spin model that reveals interesting properties of quantum magnetism, such as quantum phase transitions and quantum spin glasses, and can been used to simulate quantum annealing [76].In this paper, given the naïve implementation of the algorithm, we restrict ourselves to the one-dimensional case, for which the spectrum can be found analytically by a JW transformation from the spin model onto free fermions [75] (note that this is the reverse of the JW mapping mentioned elsewhere in this work, which refers to encoding fermionic degrees of freedom into qubits).
The Hamiltonian for the TIM with n sites and periodic boundary conditions is defined as: where Ẑn+1 = Ẑ1 .The Ising model with n sites is represented with n qubits.Under the JW transformation, the TIM Hamiltonian is composed of 2n − 1 Pauli strings.
We prepare the initial state to be an equal superposition of all spin basis states using the Hadamard gate: ⊗n .This represents a uniform prior -that is, we have not encoded any prior information about what we might expect the ground state (GS) to be.Importantly, since we consider the ferromagnetic limit, J > h, we do not benefit from a good initial guess for the GS.

Hubbard model
The fermionic Hubbard model (FHM) is the simplest possible model for correlated electrons; it approximates long-range Coulomb interactions with on-site interactions, but still exhibits a wide range of interesting phenomena, including magnetic ordering, metal-insulator transition and superconductivity [77,78].Further, the FHM exhibits correlations that are difficult to capture by classical methods.To this end, the Hubbard model is widely used as a benchmark for quantum algorithms.Again, we restrict our analysis to the one-dimensional chain under periodic boundary conditions, for which an analytic solution is known [79].
The Hamiltonian for the non-relativistic single-band fermionic Hubbard model in real space is given by: where U > 0 corresponds to repulsive on-site electron-electron interactions, and t < 0 corresponds to a lowering of the kinetic energy of the system by allowing for delocalisation over the sites.We consider the case where electrons are strongly interactingthat is, when U t ≥ 1.The sum over i, j is typically restricted to account for the exponentially decaying overlap of wavefunctions between sites.We adopt the standard model, in which hopping is only considered between nearest-neighbours.
The Hubbard model with m sites is represented by 2m qubits under the JW transformation.An occupation number basis is used to enumerate the states in the Hilbert space; for instance, the 2-site model has basis states |n 1↑ n 1↓ n 2↑ n 2↓ ⟩.The Hamiltonian under the JW transformation is comprised of 7m − 3 Pauli strings, including Î⊗2m .
Since the Hubbard Hamiltonian conserves the total number of spin up and spin down electrons, [ Ĥ, n↑ ] = [ Ĥ, n↓ ] = 0, we can consider the action of the Hamiltonian on a particular sector (n ↑ , n ↓ ) of Hilbert space [80].For a system of m sites, we consider the half-filled sector (m/2, m/2).Figure 6 gives the initial state preparation circuit; the X and cnot gates are used to excite the state into the sector (m/2, m/2), whereas the H and R y gates produce a superposition of the two antiferromagnetic states: We choose θ = π, giving an out-of-phase superposition of the two states.

B. Simulation Results
Figures 7 and 8 show the simulation results for the 4-site TIM and the 2-site Hubbard model respectively; both the expected energy and the fraction of successful measurement shots are plotted and compared to the exact results, obtained from an exact diagonalisation of the Hamiltonian, as well as to the expected evolution of the Trotterised ITE operator.There is good agreement between the expected Trotter results and the simulation results; deviations not accounted for by the error calculations can be attributed to additional stochastic errors.
We observe the expected exponential decay in the total probability of successful application -the longer the evolution time required for convergence to the ground state, the fewer the number of successful shots remaining that can be used in the final estimation of the energy.As expected, the stochastic errors in the mean energy increase throughout the evolution as the number of successful shots decreases.As the system size increases, the final success probability at convergence is reduced.After increasing the number of sites in the Hubbard model to 3 and 4 sites, the success probability at convergence decreases to 10 −5 and 10 −9 respectively, requiring a prohibitive number of shots (roughly 10 9 and 10 13 ) to obtain an accurate estimation of the ground state.This being said, we also observe an initial exponential approach towards the GS energy.

V. DISCUSSION
The PITE algorithm presented in this paper, Trotterised PITE, is shown to systematically recover the ground state for the small systems investigated, giving incentive for further investigation.The most comparable PITE proposal to the method presented in this paper uses a block encoding constructed from a LCU [35].Whilst Trotterised PITE block encodes This circuit structure can be easily generalised to generate a linear combination of the two classically antiferromagnetic states for any number of sites.
the Trotter-decomposed ITE operator, LCU PITE block encodes a first-order approximation to the ITE operator, (1 − Ĥ τ r ) r , where r is the number of time steps -importantly, this is still a ground state projector in the limit τ → ∞ [81].Similar to Trotterised PITE, the LCU algorithm uses one ancillary qubit, which is reset after each mid-circuit measurement, and unlike many of the previous PITE proposals, both methods automatically build the required circuits given any input Hamiltonian expressed in a Pauli basis.Throughout this discussion, we will make comparisons with the LCU PITE method.
The cost of the PITE algorithm arises from two separate sources: (1) the exponential decay in the success probability with the evolution time, necessitating the use of large numbers of shots, from which the successful shots are post-selected, and (2) the depth of the circuits.The first of the two is the limiting factor, as is the case with most block encoding methods.Indeed, the simulation results demonstrate that, given an initial naïve implementation of the Trotterised PITE algorithm, it is prohibitively expensive to apply it to larger systems; for the Hubbard Hamiltonian, it was even infeasible to run simulations for more than 2 sites.For the purposes of GS determination, we only require that repeated application of the operator eventually projects the state onto the GS.Many GS projectors are possible [81].Roughly speaking, the 'better' the GS projector -that is, the fewer the number of applications of the projector, r, required to reach the GS -the more 'information' is thrown away after each application, and the more rapidly the success probability decays.The results in Figures 7  and 8 suggest that, though the success probability for the ITE operator decays exponentially with simulation time, there is an initial exponential approach towards the GS, suggesting that PITE methods could be used as part of a scalable routine for GS determination.For instance, PITE could be used for a shorter length of time, and the resulting state could then be passed to the QCELS method, which only requires a squared overlap with the GS of 0.5 [82,83].This routine would make use of the fast initial decay in excited state contributions from the exponential GS projector, whilst avoiding extra reduction in success probability in the latter, slower part of the convergence.
Since all PITE proposals aim to approximate the ITE operator, they will all produce similar success probabilities at the end of the evolution, the only difference arising from the scaling factors 1 α of these block encodings.This is because the success probability at any given time is only determined by the degree of non-unitarity of the operator -that is, the success probability of a block encoding of the ITE operator is largely determined by the Hamiltonian, and to a lesser extent by the prefactor 1 α , which reduces the success probability by a factor of α 2 (7).The proposed Trotterised PITE method is optimal with respect to the scaling factor for each Hamiltonian term in the product formula.For the exact ITE operator, the scaling factor is not optimal but can be increased by reducing the 1-norm of the qubit Hamiltonian, a task that is of interest to several Hamiltonian simulation procedures; to this end, significant work has already been carried out [70].
The most significant improvements to the success probability, however, would likely be achieved by running the PITE algorithm in tandem with some form of AA procedure [30,33].It is important to appreciate that the design of this AA procedure is a non-trivial task and should be tailored to the specific PITE procedure in order to minimise circuit depths; for instance, Grover's algorithm was implemented with good success by Nishi et al for the LCU PITE proposal [30].
The final success probability in a PITE procedure largely depends on the evolution time required for convergence to the ground state τ , which is a property inherent to the system, since it is fixed by the energy difference between the ground state and first excited state of the system (Equation 39).However, τ is also dependent on the overlap between the initial guess wavefunction and the GS.When using a Jordan-Wigner transformation, which maps each spin orbital onto a qubit, subsequent entanglement of the qubits will encode a Full Configuration Interaction (FCI) wavefunction.It is important to note that, since this work presents a proof of principle implementation of a Trotterised PITE procedure, the initial states used in this work were particularly poor.Typically, the Hartree-Fock (HF) state is used as the initial guess: as the system size grows, the contribution of the HF state to the ground state diminishes [84], and τ must correspondingly increase.Indeed, for strongly correlated systems, there are simulations which suggest that the overlap decreases exponentially with system size [85].Reducing τ requires choosing an initial state with a larger overlap with the GS: for instance, by running an initial classical simulation to screen for the most important amplitudes in the FCI expansion.
Filip et al. used this approach to screen for important amplitudes in a Unitary Coupled Cluster Ansatz, used in a variational quantum eigensolver routine [67].On quantum hardware, we have the additional complexity of needing to find methods of efficiently preparing the desired initial state [86].
The mitigation of the exponential decay in the success probability is the most important obstacle for the implementation of PITE methods.However, we should also consider methods for reducing the circuit depths.Aside from being proportional to τ , the minimum number of time steps required is also inversely related to the maximal time step ∆τ that will still guarantee convergence to the ground state.This is, in turn, determined by the error in the approximation used.Trotterisation is a very popular simulation method due to its simplicity, flexibility and efficiency; as such, it is immediately possible to reduce the circuit depths of the Trotterised PITE method by appealing to the wealth of techniques already developed by the Trotter decomposition community.
For instance, applying higher order decompositions [52] and randomisation protocols [54,[58][59][60][61], grouping together Hamiltonian terms to minimise the number of mutually commuting partitions [10,38,48,63,64], and reordering terms to maximise the number of gate cancellations [39].Further, until recently, the leading method for reducing the cost of Hamiltonian simulation for quantum chemistry Hamiltonians was tensor hyper contraction, which uses non-unitary rotations and therefore restricts the method to LCU approaches [65].There are now superior double factorisation techniques which use unitary rotations [66] and can be combined with Trotter decompositions and its non-unitary variants, including Trotterised PITE.
Trotterised PITE produces shorter circuits than LCU PITE, at the cost of requiring more mid-circuit measurements.LCU PITE uses 2M controlled Pauli gadgets per time step, where M is the number of Pauli strings comprising the Hamiltonian, followed by a single mid-circuit measurement.On the other hand, Trotterised PITE uses M (not controlled) Pauli gadgets per time step, along with M mid-circuit measurements.Since these Pauli gadgets are not controlled, they can be executed with higher fidelity.Note that the requirement for more mid-circuit measurements in our algorithm is not a restriction for some quantum computer architectures, like those based on trapped ions, as they are able to execute mid-circuit measurements on a similar timescale as gate operations.
As mentioned, for the purpose of GS determination, the operator only needs to be a GS projector.However, for the purpose of computing thermal averages, it is important to be able to approximate the ITE operator to high accuracy [22,87].LCU PITE is a block encoding of a first-order approximation to the ITE operator.Thus, regardless of the accuracy with which the constituent controlled RTE operations are performed, the dominant error will scale as (∆τ ) 2 per time step.On the other hand, applying the above methods to the Trotterised PITE method can not only reduce circuit depth, but can also be used to improve the overall error in the approximation to the ITE operator.
The number of Hamiltonian terms may be minimised with a sensible choice of basis.Whilst we have used a second quantised basis throughout this paper, giving rise to a qubit representation of the Hamiltonian comprised of Pauli strings, it has been argued that a first quantised basis may produce more efficient scalings [35].Any ITE algorithm using a Pauli representation for the Hamiltonian can easily be extended to the first quantised basis using the protocol proposed by Kosugi et al. [35,88].
Aside from improving the success probability, and thus reducing the number of shots required, we should also consider reducing the time taken for each shot to run.We note that improvements in hardware to provide a quit-if-fail functionality would remove the need to run failed shots to the end of the circuit, greatly reducing the average run-time across shots.

VI. CONCLUSIONS
We have developed a purely quantum routine for performing probabilistic imaginary time evolution (PITE), based on a Trotter decomposition of the Hamiltonian.The block encoding suggested can be thought of as a modification of the Pauli gadget primitive, which is an efficient and widely used circuit implementation for real time evolution.PITE algorithms avoid many of the limitations of near-term algorithms.Namely, they avoid the restriction on the accuracy that can be achieved as a consequence of using a fixed ansatz in variational methods, including VITE, and any restrictions placed on the locality of the Hamiltonian, as is the case in QITE.
In this paper, we have implemented the Trotterised PITE block encoding with a simple initial protocol: successful shots are post-selected, no amplitude amplification procedure is used, and minimal optimisation for the number of mid-circuit measurements required, and the depth of the circuits, is used.We applied this routine to the task of ground state determination in one-dimensional for the Transverse Ising Model with 4 sites and the fermionic Hubbard model with 2 sites.In particular, we found that the number of shots required by this naïve implementation became prohibitively high even for a 4-site fermionic Hubbard model.We argue that this behaviour is expected given the nature of the PITE approach.The limiting factor in the performance of all PITE algorithms is that they exhibit an exponential decay in the probability of successful application with the number of mid-circuit measurements applied -indeed, this problem is shared by many block encoding methods.Importantly, the algorithm successfully recovers the ground state for small systems, giving incentive to adapt its structure to overcome this limitation.We discussed a multitude of strategies that could be applied to reduce the run-time of the algorithm, of which the inclusion of an amplitude amplification procedure is likely to be the most significant contribution [30,33].Further to this, Trotterised PITE could be used to initially increase the overlap with the ground state as part of a more scalable routine for ground state determination: (1) Trotterised PITE would first be run for a shorter length of time, during which many of the excited state contributions decay, (2) the resulting state is then passed to a quantum subspace diagonalisation method [22] or the QCELS algorithm, which only requires a squared overlap with the GS of 0.5 [82,83].This routine would make use of the fast initial convergence of the exponential GS propagator, whilst avoiding unnecessary reduction in success probability in the latter, slower part of the convergence.ITE methods implemented on quantum circuits have found many applications in recent years.Sokolov et al used ITE for the determination of the ground state of transcorrelated Hamiltonians [89].They obtained promising results, gaining up to four orders of magnitude improvement in the absolute energy error in comparison to non-transcorrelated approaches.More generally, algorithms that perform imaginary time evolution can be used as a subroutine; for instance, in the estimation of low-lying excited states using quantum subspace diagonalisation methods [22].The most comparable method to Trotterised PITE is LCU PITE, which block encodes a first order approximation to the ITE operator.Since Trotterised PITE is a block encoding of the Trotter-decomposed ITE operator, whose accuracy can be readily improved, it is more suited than LCU PITE to applications which require an accurate application of the ITE operator; for instance, in the calculation of finite temperature correlation functions [22,87].Moreover, we note that our algorithm can be readily extended to the real (or imaginary) time Trotterised simulation of non-Hermitian systems, Ĥ = Ĥ1 + i Ĥ2 , where Ĥ1 and Ĥ2 are Hermitian, by concatenating Pauli gadgets for real time propagation (deterministic, Figure 10) and for imaginary time propagation (probabilistic, Figure 5).This could be used for the simulation of open systems (see, for instance, Algorithm I in [90]).
The simplest operator of this form is Ẑ⊗n .Figure 9 exemplifies a quantum circuit realisation of this operator; circuits of this form are commonly known as 'phase gadgets' [9].In this circuit, the 'ladder of cnots' computes the parity of each of the computational basis states and encodes it in the state of the final qubit, setting its state to be |0⟩ for an even parity and |1⟩ for an odd parity.A Z-gate is then applied to this qubit, imparting a phase of +1 for an even parity and −1 for an odd parity, as required.Finally, the parity computation is reversed.where Bm k denotes the basis transformation operator from the σm k -basis to the computational basis.For instance, BX = H and BY = SH, where H and S denote the Hadamard gate and the phase gate, respectively.Equation (A.2) allows us to reduce the problem of simulating any Pauli string to the problem of simulating Ẑ⊗n .Once the phase gadget primitive is combined with the necessary basis transformation gates, the circuit is referred to as a 'Pauli gadget'.According to Equation (A.1), the effect of the RTE operator e −it Ẑ⊗n on each computational basis state is to add a phase conditioned on its parity, e −itpi .This can be achieved by modifying Figure 9 to apply a Z rotation gate, R z (2t), on the parity-storing qubit.
An example quantum circuit for the real time evolution of the Pauli string XY Z is shown in Figure 10.The Pauli gadget is repeated for each of the Pauli strings, using their corresponding basis transformation gates (A.2), and for each of the Trotter steps, according to the Trotter decomposition (15).

Figure 1 .
Figure 1.The non-unitary nature of N1(a) projecting out the magnitude of |0⟩ relative to |1⟩ upon successful measurement of |0⟩ on the ancilla.

M k=1 |c k | is the 1 -
norm of the Hamiltonian.Using the following relation between the spectral norm and the 1-norm of the Hamiltonian expressed in a Pauli basis:

Figure 6 .
Figure 6.Initial state preparation circuit for the 2 site 1D Hubbard model.This circuit produces the state 1

Figure 7 .Figure 8 .
Figure 7. Energy estimate (left panel), relative energy compared to the true GS energy, E0 (centre panel) and success probability (right panel) plotted as a function of the number of time steps for the 4-site TIM Hamiltonian, with J = 0.5, h = 0.1(42), and ∆τ = 0.1.The red curves are the exact results obtained from an exact diagonalisation of the Hamiltonian, the green curves are the expected results from the Trotterised ITE operator, and the black points have been evaluated with the Trotterised PITE algorithm using 10 5 measurement shots.