Phase-Sensitive Quantum Measurement without Controlled Operations

Many quantum algorithms rely on the measurement of complex quantum amplitudes. Standard approaches to obtain the phase information, such as the Hadamard test, give rise to large overheads due to the need for global controlled-unitary operations. We introduce a quantum algorithm based on complex analysis that overcomes this problem for amplitudes that are a continuous function of time. Our method only requires the implementation of real-time evolution and a shallow circuit that approximates a short imaginary-time evolution. We show that the method outperforms the Hadamard test in terms of circuit depth and that it is suitable for current noisy quantum computers when combined with a simple error-mitigation strategy.

Introduction.-Thecomplex phases of quantum amplitudes play an essential role in quantum algorithms [1][2][3][4][5][6] and quantum sensing [7].Many algorithms require measuring the relative phase between two quantum states [8][9][10][11][12][13][14][15][16][17].A common subroutine for this purpose is the Hadamard test, which converts phase information into probabilities by means of interference [18].Despite impressive experimental progress, the Hadamard test remains out of reach for most applications owing to the challenge of implementing the required controlledunitary operation.In this Letter, we propose an alternative method to determine the complex overlap between certain states that uses no ancillary qubits or global controlled-unitary operations.Unlike other ancilla-free schemes [12,19], our approach does not require the preparation of superpositions with a reference state, which are highly susceptible to noise [20][21][22][23][24][25].Instead of interference, our method hinges on the principles of complex analysis.
The proposed approach applies to overlaps of the form of the (generalized) Loschmidt amplitude where H is a local Hamiltonian.Our algorithm requires that |ψ⟩ has a short correlation length and that |ψ ′ ⟩ can be prepared using a unitary circuit.These assumptions are needed to be able to efficiently apply a short imaginary-time evolution to |ψ⟩ [26][27][28][29][30][31][32][33][34] and to perform a projective measurement onto the final state |ψ ′ ⟩ [35].
The absolute value |G(t)| can be obtained by repeatedly evolving |ψ⟩ and averaging over projective measurements onto |ψ ′ ⟩.Here we describe how to obtain the phase.Equation (1) includes several cases of interest.When |ψ ′ ⟩ = |ψ⟩, G(t) is the Fourier transform of the local density of states, which has applications in the study of quantum chaos [36,37], in optimal measurements of multiple expectation values [14], and in estimating energy eigenvalues [9,10,13,15,17].The case when |ψ ′ ⟩ = Ae −iHt ′ |ψ⟩, for a local unitary A, is relevant for probing thermal properties of many-body systems [12,[38][39][40].The key idea underlying our method is to view G as a function of a complex variable z.Assuming that G(z) is analytic and nonzero, the Cauchy-Riemann equations imply that the real-time derivative of the phase of G(z) is equal to the derivative of ln |G(z)| along the imaginary-time direction.We use this relation to obtain the desired phase by carrying out the following three steps on a quantum computer (see Fig. 1).First, a quantum circuit applies an evolution under the Hamiltonian H for a short imaginary time h to the initial state |ψ⟩ [26][27][28][29][30][31].Second, we evolve the system under H for the real time t.Third, we perform a projective measurement onto the state |ψ ′ ⟩ by inverting the circuit that prepares |ψ ′ ⟩ from a computational basis state, followed by measurements in the computational basis.Using these steps, we can estimate |G(t ± ih)|.This yields a finite-difference approximation to the imaginary-time derivative of ln |G(z)|, which is equal to the real-time derivative of the phase.We finally compute the phase of the Loschmidt amplitude by repeating these steps for different values of t and numerically integrating the derivative, starting from a time at which the phase is known.
We show below that our method is efficient if |G(t±ih)| is bounded from below by an inverse polynomial in the system size N .For product states |ψ ′ ⟩ = |ψ⟩, however, the Loschmidt amplitude decays as a Gaussian function on a time scale O(1/ √ N ) [41].In this case, our approach will be inefficient even for short constant times, for which the Loschmidt amplitude can be computed by a polynomial-time, classical algorithm [42].By contrast, no efficient classical algorithm is known for the case |ψ ′ ⟩ = Ae −iHt ′ |ψ⟩.Since the real-and imaginary-time evolution operators commute, our method can be used to compute the phase as a function of t − t ′ , with the expectation value ⟨ψ|e iHt ′ Ae −iHt ′ |ψ⟩ serving as the reference for the integration.This is expected to be classically hard even for small t − t ′ since computing the reference value at times t ′ = poly(N ) is BQP-complete [43].
Although our approach is based on the analytic properties of a function of a continuous variable, we show below that it also works well in the discrete setting of Trotter evolution.Hence, the method applies to both circuitbased quantum computers and to analog quantum simulators supplemented by shallow circuits to implement the imaginary-time evolution.To demonstrate the suitability of the method for near-term quantum devices, we combine it with a simple error-mitigation strategy [44][45][46][47][48][49] and show numerically that the phase can be reliably recovered in a system of N = 24 qubits.Beyond providing a viable alternative to the Hadamard test on near-term quantum computers, our method may be useful in the early fault-tolerant regime as the absence of controlled global operations significantly reduces the circuit depth.
Theoretical approach.-Toformally describe the algorithm, we consider the complex variable z = t − iβ, where t represents real time and β stands for imaginary time or inverse temperature.The generalized Loschmidt amplitude, Eq. ( 1), can be decomposed into its absolute value and phase according to where 0 ≤ r(z) ≤ 1 and ϕ(z) is real.In a system of finite size, the Loschmidt amplitude can be written as a sum of exponentials by expanding the states |ψ⟩ and |ψ ′ ⟩ in the energy eigenbasis.The logarithm ln G(z) is therefore holomorphic everywhere except when G(z) = 0.For an analytic branch of ϕ(z), the Cauchy-Riemann equations applied to ln G(z) = ln r(z) + iϕ(z) give Therefore, if G(t) ̸ = 0 in the interval [t 1 , t 2 ], the phase difference ϕ(t 2 ) − ϕ(t 1 ) can be computed as If the phase ϕ(t 1 ) is known, then ϕ(t 2 ) may be computed from the partial derivative of r(z) along the imaginarytime direction.In practice, we numerically approximate the partial derivative by the mid-point formula where h is a small parameter.This procedure is well defined for r(t ± ih) > 0 in the interval [t 1 , t 2 ].To bound the computational errors, we make the stronger assumption that | ln r(z)| ≤ cN at all points in the complex plane within distance a of the interval [t 1 , t 2 ], for constants c and a > h.In the case of Trotter evolution, we make analogous assumptions for closely related functions [50].We highlight, however, that our approach can be extended to treat zeros in G(t) by separately considering the resulting discontinuities in the phase [50].
The above analysis has reduced the problem to measuring the absolute values r(t ± ih) 2 .It involves nonunitary imaginary-time evolution which cannot be directly applied.However, Motta et al. [27] showed that e ±hH can be simulated by a unitary circuit for short times h if the spatial correlations of |ψ⟩ decay exponentially with correlation length ξ and H is a local Hamiltonian.For each local term H m in the Hamiltonian, it is possible to approximate e ±hHm |ψ⟩ ≈ c ± m V ± m |ψ⟩, where V ± m is a local unitary and c ± m = ⟨ψ|e ±2hHm |ψ⟩ accounts for the normalization.Since h is small, the product over all V ± m /c ± m (in arbitrary order) is a good approximation of e −hH .The operators V ± m act on O((ξ log N ) d ) qubits in d spatial dimensions and the complexity of computing V ± m is quasi-polynomial in N for d > 1.This renders the approach challenging for large ξ.Below, we focus on the simplest case of product states, for which the unitaries V ± m act on the same sites as H m and can be efficiently computed.The resulting circuit has the same structure as a single Trotter step [50].Error analysis.-Wenext analyze the error in the estimated phase arising from the different approximations in our algorithm.The approximation error of the imaginary-time evolution is dominated by the first-order Trotter decomposition, which results in the phase error [50] ∆ϕ The factor t = t 2 − t 1 accounts for the accumulation of errors in the integral in Eq. ( 4).While the real-time evolution can be carried out exactly on analog quantum simulators, digital quantum computers incur an additional Trotter error, leading to the phase error [50] ∆ϕ Here, τ is the time of a single Trotter step, p is the order of the Trotter decomposition [51], and we again included the accumulation of errors in Eq. ( 4).Numerical differentiation and integration give rise to additional errors.They can, however, be safely ignored for practical orders of the Trotter expansion (p ≤ 4) as they are asymptotically at most as big as ∆ϕ ITE and ∆ϕ RTE [50].
We verify these analytic estimates using numerical results for the transverse-field Ising chain, whose Hamiltonian is given by Throughout this work, we set J = 1 and g = 0.5, corresponding to the ferromagnetic phase.Both states |ψ⟩ and |ψ ′ ⟩ are chosen as |↑↑↑ • • •⟩.For the Trotter decomposition, we alternate between the ferromagnetic and transverse field terms.

Method D M
Hadamard test O(t This work O(r Figure 2 shows the error in the phase computed using our approach.The numerical results were obtained by matrix product state simulations with bond dimension 200, for which truncation errors are negligible [50].In Fig. 2(a), we set τ = 0.01 such that the error in the imaginary-time evolution dominates.The phase error collapses onto a single curve upon dividing by N h 2 , which confirms the predicted error due to imaginary-time evolution, Eq. ( 6).Similarly, we set h = 0.01 in Fig. 2(b) to isolate the effect of the real-time Trotter error.The collapse of the data agrees with Eq. ( 7).
In addition to numerical errors, any experiment incurs statistical errors.Given M measurements, a probability p estimated by counting successful outcomes will have a multiplicative error (1 − p)/M p, governed by the standard deviation of the binomial distribution.According to Eq. ( 5), for the measured probabilities p to the final phase for M measurements per time step.The integral in Eq. ( 4) is included in the factor In contrast to the previous errors, the statistical error depends on the magnitude of the measured probabilities.
Comparison with existing methods.-Tocompare our approach to existing methods, we consider the error ∆G in the complex Loschmidt amplitude G.This error is related to the phase error, ∆ϕ, by Here, ∆r is the error from an independent measurement of r, which only requires the Trotterized circuit without imaginary-time evolution.To bound ∆r by ϵ, we need a circuit of depth ) and a number of M r = O 1/ϵ 2 measurements [50].For the term r∆ϕ, we bound the individual contribu- tions to the phase error.For instance, r∆ϕ ITE < ϵ implies that h = O( ϵ/rN t).A similar bound on the real-time evolution gives τ = O((ϵ/rN t 2 ) ). Bounding the statistical error yields the number of measurements M = O(I 2 r 3 t 3 N/ϵ 3 ) for each time step.We note that when r is bounded from below by a constant, the cost of estimating ϕ dominates.
We compare this resource cost to the Hadamard test and sequential interferometry [12].The latter method employs a reference state whose Loschmidt amplitude, including the phase, is known.The details of these two methods are described in the supplemental material [50].Table I summarizes the resource cost for each method.For a constant evolution time t, the circuit depth needed for our algorithm is reduced by a factor O(N 1+1/d ) compared to the Hadamard test with swaps, and by O(N 1/p ) compared to sequential interferometry.This improvement is particularly significant for noisy quantum computers, for which circuit depth is the key limiting factor.
Applications.-Forpractical applications of our protocol, it is important to consider the role of noise.We propose a simple rescaling strategy based on previous work to mitigate the effects of noise [49].The key observation is that errors are unlikely to drive the system towards the target state |ψ ′ ⟩.Hence, the measured probabilities are decreased in a consistent fashion, which can be mitigated by rescaling with the probability of having no noise.This is equivalent to zero-noise extrapolation with an exponential fitting function [44,52,53].Below, we simply use the known noise rate for rescaling.In practice, the rescaling factor can be determined by enhancing the noise or by measuring the survival probability after forward plus backward evolution [49].
As a proof of concept, we apply our approach to compute the local density of states (LDOS) d(E) through the Fourier transform If the initial state has a sufficiently large overlap with the ground state, its LDOS enables determining the groundstate energy.In quantum chemistry, this can hold even for product states, rendering our approach particularly suitable [54].We further note that our approach is compatible with recent proposals that classically process the Loschmidt amplitudes at different times in order to solve the general quantum eigenvalue estimation problem [10,50] with Heisenberg-limited scaling [13,15,17].We apply our approach to compute the LDOS of |ψ⟩ = |↑↑↑ • • •⟩ in an Ising chain of system size N = 24.We numerically carry out the Trotter evolution with Trotter steps τ = h = 0.3 using the Cirq library [55].We add single-qubit depolarizing noise of rate γ = 3 × 10 −3 after each layer of the quantum circuit.We average over 5000 trajectories of a Monte Carlo wavefunction simulation to obtain the probabilities p ± .The results are shown in Fig. 3.Here we have included statistical noise by simulating the experimental sampling procedure (see caption).
Figure 3(a) shows that the depolarizing noise is mitigated well by rescaling r 2 (t) and r 2 (t ± ih) by (1 − γ) N D .The error in the reconstructed Loschmidt amplitude remains small within the range of t in Fig. 3(b).We estimate the LDOS of the initial state by a discrete Fourier transform of the data in Figure 3(b) and similar data for the imaginary part of G(t).The energy resolution is π/t max ≈ 0.31, determined by the maximum time t max = 10.We show the result in Fig. 3(c) for both noisy, error-mitigated (green) and noiseless (orange) Trotter simulations.For reference, we also include the exact result (black line), which is broadened by a Gaussian of width 0.08.For both Trotter simulations, the first point with d(E) > 0.1 appears at E ≈ −7.50, while the exact ground state energy is E 0 ≈ −7.55.
Summary and outlook.-Wepropose a quantum algorithm to estimate the phase of Loschmidt amplitudes applicable to states with short-ranged correlations.It can replace and outperform the Hadamard test for amplitudes that arise from continuous time evolution under a local Hamiltonian.While our analysis focused on generalized Loschmidt amplitudes, the approach can be readily extended to multiple time-evolution operators [50], which renders it applicable to many quantities of physical interest including transport coefficients [56][57][58] and out-of-time-ordered correlators (OTOCs) [59,60].The algorithm requires no ancillary qubits or controlled operations.When combined with a simple error-mitigation strategy, our algorithm may enable phase-sensitive measurements on current noisy quantum devices for system sizes that out of reach for other methods.

Supplemental Material
A. Errors

Zero-free functions
To apply our algorithm to the time interval [t 1 , t 2 ], we require that G(t) is nonzero in this region.In practice, we have to place a lower bound on the magnitude of G(t) to guarantee a bounded error of the algorithm.We will refer to functions that satisfy such a lower bound as zero free, following the terminology introduced in reference [69].Owing to its similarity to the partition function, the Loschmidt amplitude generically takes the form G(z) ∼ e −N g(z) in the limit of large N for some function g(z) [63,65,66,69].This behavior motivates the following definition of a zero-free function.
Definition 1 (Zero-free functions).A sequence of functions f N (z) on the complex plane is called zero free at point z 0 if there exist constants c and a such that Assuming that a function is zero free allows us to bound to the derivatives of ln f N (z) according to the following lemma.
Lemma 1.Consider a sequence of holomorphic functions f N (z), which are zero free at t 0 ∈ R. We let z = t − iβ with t, β ∈ R and log f N (z) = log r N (z) + iϕ N (z) for some analytic branch of the logarithm.Given any constant m, the magnitude of the partial derivatives ∂ m ln r N (z)/∂β m and ∂ m ϕ N (z)/∂t m at z = t 0 is bounded from above by O(N ).
Proof.To prove the lemma for ∂ m ln r N (z)/∂β m , we observe that The right-hand side can be expressed using Cauchy's integral formula as It remains to bound the magnitude of ln f N (ζ) on the circle |ζ − t 0 | = a/2.However, the magnitude depends on the choice of the branch of the logarithm.To overcome this issue, we consider the function fN (z) = e −iϕ N (t0) f N (z).We can bound the Cauchy integral (A3) in terms of ln fN (z) instead of ln f N (z) because the derivatives of the two functions are equal.
To complete the bound, we use the Schwarz integral formula to express ln fN (z + t for any |z| < a.Since Im ln fN (t 0 ) = 0, this yields By assumption, we have By substituting these results into (A2) and (A3), we obtain For any constant value of m, the derivative of ln r N is thus bounded by a quantity O(N ).The proof for ∂ m ϕ N (z)/∂t m follows the same steps, starting from |∂ m ϕ N (t)/∂t m | t=t0 ≤ |d m /dz m ln f N (z)| z=t0 in place of Eq. (A2).
For simplicity, we omit the subscript N below.

Error in numerical derivatives and integration
As a first illustration of the implications of the zero-free condition, let us deal with errors arising from numerical derivatives and integration.For the symmetric finite difference approximation to the derivative the error is given by [62] ∆ϕ where the N dependence directly comes from Lemma 1, assuming that G(z) is zero free.This error will be multiplied by t in our algorithm, following the integration in Eq. ( 4) of the main text.We note that the error has the same scaling as the error ∆ϕ ITE due to the approximate imaginary-time evolution.Therefore, a higherorder finite difference approximation will not improve the asymptotic scaling of the total error.
To compute the phase up to time t, we need to carry out the integral in Eq. ( 4) of the main text.In practice, we can only estimate the rate of change of the phase at a discrete set of points such that the integration is necessarily approximate.We use the Newton-Cotes formula of degree n with a fixed spacing τ between the samples, as is natural for a real-time Trotter step τ .The corresponding numerical integration error is [62] ∆ϕ where s equals n + 2 rounded down to the closest even number (n = 1, s = 2 for the trapezoidal rule; n = 2, s = 4 for Simpson's rule).The dependence on N again follows from Lemma 1.The Trotter error in the real-time evolution, ∆ϕ RTE , typically dominates over the error in the numerical integration, as the coefficient s = 4 for Simpson's rule exceeds the order p of practical low-order Trotter expansions.

Trotter errors
Trotter errors lead to an error in the estimation of the quantities The measured probabilities will have an additive Trotter error that scales as O(N h 2 ) for imaginary-time evolution and as O(N tτ p ) for real-time evolution [51].However, to bound the error on ∂ ln r/∂β, it is necessary to control the multiplicative error in r.This is challenging because r may be exponentially small in the system size.Nevertheless, we show in this section that the above error scalings also apply to the numerical derivative in Eq. (A8) assuming that particular functions are zero free.We require that G(z), its Trotterized version G ITE (z) defined in Eq. (A17), and the function F t0 (z, w) in Eq. (A31) are both zero free.We highlight that if these assumptions are violated, it may nevertheless be possible to compute the phase using the correction method described in Appendix C. When expanded as the Taylor series in h, the finite difference in Eq. (A8) has the form The O(N ) dependence in the last line comes again from applying Lemma 1 to G(z).We will us this form below.
For simplicity, we consider Trotter errors in the imaginary-time and real-time evolution separately.It is straightforward to show that the individual errors add when imaginary-time and real-time evolution are Trotterized at the same time.

a. Imaginary-time evolution
In this section, we provide bounds for the errors incurred by replacing the imaginary-time evolution by a product of local unitaries.We obtain different bounds for product states and for states with exponentially decaying correlations.For product states, the imaginary-time evolution can be realized by local unitaries that have the same support as the local terms appearing the Hamiltonian.The only error in this case arises from Trotterization.For correlated states, more complicated unitaries that act on larger subsystems are required.There are additional errors that depend on the size of the subystem on which the unitaries act.
We consider a decomposition of the Hamiltonian H into Γ = O(N ) local terms, To implement the imaginary-time evolution, we wish to find unitary operators U j (h) that satisfy where c j = ⟨ψ|e −2hHj |ψ⟩ ensures normalization.Although these unitary operators are not unique, we may always choose them such that U j (h) is a smooth function of h with U j (0) = I.
To bound the errors due to the imaginary-time evolution, it is necessary to be more explicit about the dependence U j (h) on h.To this end, we let such that V j (h) = e −iBj h equals U j (h) to first order in h.By differentiating Eq. (A16), it follows that the Hermitian operator B j satisfies This equation always has a (non-unique) solution as it specifies a single column of the matrix representation of B j in an orthonormal basis that includes |ψ⟩ as a basis vector.
For a particular choice of B j , we replace the imaginarytime evolution by Γ j=1 c j V j (h).The corresponding Loschmidt amplitude is given by The unitary operations can be directly realized experimentally and it is straightforward to keep track of the constants c j classically.We note that the ordering of the operators can be chosen arbitrarily as it does not affect the result to lowest order in h.Indeed, one can verify by direct calculation that where To bound the error of the finite-difference A ITE (t, h) = [ln r ITE (t − ih) − ln r ITE (t + ih)]/2h, it is necessary to derive a bound on the higher derivatives of ln r ITE (z).However, G ITE (z) is not holomorphic due to the separate dependence on the real and imaginary parts of z.To be able to apply Lemma 1, we define which is holomorphic in z and satisfies If we assume that f t (z) is zero free at z = 0, then Lemma 1 implies from which it follows that When the initial state |ψ⟩ is a product state, the support of U j (h) can be restricted to the sites on which H j acts.Therefore, B j can be computed efficiently on a classical computer.For more general states, it becomes impractical to solve Eq. (A16) exactly.Nevertheless, if the correlations of the state |ψ⟩ decay exponentially, then U j (h) can be approximated by a local unitary operator Ũj (h) such that According to Theorem 1 in the supplemental information of the paper by Motta et al. [27], it is sufficient for Ũj (h) to act on qubits, where d is the spatial dimension and ξ is the maximum correlation length of the sequence of states k j=1 Ũj (h) |ψ⟩ for 0 ≤ k ≤ Γ.Since the value of h is small, we expect that the correlation length of |ψ⟩ is typically a good estimate of ξ.
Replacing the unitary operators U j (h) by their local counterparts Ũj (h) results in a modified Loschmidt echo rITE (z) whose error is bounded by This leads to an error in the finite difference given by which is to be added to the discretization error of Eq. (A22).For both errors to be bounded by ϵ, we let To obtain Ũj , tomography on N q sites is required, which incurs a computational cost that scales as exp(O(N q )).The cost is polynomial in the system size N and the inverse of desired error 1/ϵ for one-dimensional systems and quasi-polynomial in d ≥ 2 dimensions, provided r(t) is bounded from below by an inverse polynomial in N .

b. Real-time evolution
For the real-time evolution, we consider where U RTE (t) = U p (t/D) D is p-th order Trotter decomposition of the real-time evolution with a total number of D Trotter steps.We define the multiplicative error operator M (t) by The error operator is bounded in operator norm by ∥M (t)∥ = O(N tτ p ) [51], where τ = t/D [51].We define A RTE (t, h) as the approximation to the finite difference where with |ψ(z)⟩ = e −iHz |ψ⟩.In contrast to the imaginarytime evolution, the first-order derivative does not cancel.It contributes the leading-order error, which we analyze in what follows.
We would like to again apply Lemma 1, which leads us to define where w is an independent complex variable from z.With this definition, g(z) = |1 + d dw (ln F t )(z, 0)| 2 .Let us assume that F t0 (z, w) is zero-free for both z and w at z = t 0 and w = 0, i.e., ln |F t0 (z, w)| ≤ cN for all z, w such that |z − t 0 | < a z and |w| < a w .Since ∥M (t 0 )∥ = O(N t 0 τ p ), it is natural to choose a w = O(1/(t 0 τ p )), while keeping a z = O(1) and c = O(1).According to Lemma 1, we get Finally, we have to confirm that the partial derivative in Eq. (A29) does not change the system-size dependence.The leading term at t = t 0 is 1 2 Similar to the proof of Lemma 1, Eq. (A32) can be bounded as Re ∂ ∂β Thus Eq. (A32) still scales as cN/a z a w = O(N t 0 τ p ) and we find that, to leading order, Estimating the magnitude r An individual circuit to measure the magnitude r is required for our method and for sequential interferometry.It consists of Trotterized real-time evolution, which incurs a standard p-th order Trotter error [51] and statistical error where M is the number of measurements.To bound this error within ϵ, we require a circuit depth D r = O(t/τ ) = O(t 1+ The Hadamard test is the standard method to compute the real and imaginary part of ⟨ψ|U |ψ⟩ for a given unitary U and initial state |ψ⟩.The circuit of the Hadamard test is shown in Fig. 4. One first applies a Hadamard gate to an ancillary qubit.It is followed by a controlled unitary c-U acting on the prepared state |ψ⟩ conditioned on the first ancillary qubit.Finally, apply the Hadamard gate again to the ancilla and measure this qubit.The probability of measuring 0 is To measure the imaginary part, we modify the circuit by adding a S † = 1 0 0 −i gate after the first Hadamard gate. It is hence possible to infer G(t) directly from the measured probabilities.Shot noise gives rise to the statistical error where M is the number of measurements.Compared to our proposed method, there is no need of integration, so the Trotter error does not have an additional t dependence: The total error is hence given by Note that this only gives the phase of a single time t whereas our algorithm measures the phase on an interval [0, t].
To bound the error ∆G < ϵ, we need a Trotter step τ = O((ϵ/N t) Original algorithm corresponds to our algorithm described in main text.In the "phase corrected" results, we added a discrete phase jumpy π as well as a phase shift δ, which ensures continuity of the first derivative of G(t).
and thus In practice, with finite time resolution (e.g., the Trotter time), the gradient of the phase becomes singular near t 0 (or even some point where G(t) is almost zero).Even if we have included the (−1) n0 phase jump, this singularity can cause an additional phase factor e iδ when going across t 0 .Thus our algorithm will give a numerical result G(z), where It is possible to correct the phase error δ by requiring that the n 0 -th derivative of G be continuous.In particular, we can directly estimate δ from the expression In the discrete time setting, we evaluate the limits at the closest points on either side of the zero.An example of the algorithm in the presence of small values of r(t) is shown in Fig. 5.While there are no exact zeros in the Trotterized simulation, r(t) becomes very small at the dashed lines in panel (a).In the original algorithm, this leads to large phase jumps because the discretization is too coarse.As shown in panels (b)-(c), it is possible to partially correct these jumps using the method described above, where we impose continuity of the first derivative close at the small values of r(t).

D. Extension to multiple time evolution operators
It is possible to extend our algorithm to multiple time evolution operators.For example, one may consider functions of the form where the times in both evolutions are the same and A can be any unitary operation.Here, G A (t) is again an analytic function of t, so that the phase can directly be extracted by computing the derivative along the imaginary direction.The imaginary-time evolution can be simultaneously applied on both sides under the condition both |ψ⟩ and |ψ ′ ⟩ are short-range correlated.
Our approach can also be generalized for n evolution operators with different time variables: where U j = e −iHj tj and unitaries O j are local for j > 1.One can compute these quantities using the following protocol.
• First switch on only U 1 , i.e., set t 2 = • • • = t n = 0. Our algorithm works since the new intial state • Then switch on U 2 as well.
For each fixed t 1 , perform our algorithm with initial state O 2 • • • O n−1 |ψ⟩, evolution operator e −iH2t2 and final state e iH1t1 |ψ ′ ⟩.The phase at t 2 = 0 has already been determined in the previous step.
• Switch on the rest of the evolution unitaries one by one.

E. Applications
Our algorithm is well suited to tackle a broad spectrum of problems within the domains of quantum chemistry and condensed matter physics.Below we focus on three particularly important examples: (A) probing thermal expectation values, (B) quantum eigenvalue estimation, and (C) computing temporal correlation functions.

Probing thermal properties
Recently, a quantum algorithm was developed to probe finite-energy properties of quantum many-body systems through spectrum filtering [12].In fact, this approach has already been implemented in a proof-of-principle experiment exploring a small instance of the Fermi-Hubbard model [39].The main idea is to define the Gaussian energy filter and apply it to an initial (product) state |ψ⟩ whose mean energy is close to E. We note that it is possible to reach any constant energy density above the ground state by considering the tensor product of "patches", where larger patches allow for lower energies.For a given unitary observable A, the microcanonical ensemble average at the corresponding energy can be obtained from for sufficiently small values of δ.Alternatively, one can filter the whole spectrum as and use classical Monte Carlo method to compute the trace.
The filter P δ (E) can be approximated with a time series and thus the algorithm only requires the measurement of quantities of the form ⟨ψ|e iHt1 Ae −iHt2 |ψ⟩ or ⟨ψ|Ae −iHt |ψ⟩ . (E4) These quantities can efficiently be computed using the strategies outlined in section D.

Quantum eigenvalue estimation
The estimation of quantum eigenvalues, including the ground and excited state energies of a given Hamiltonian, is a crucial problem in quantum chemistry and condensed matter physics.Even though this problem can be solved using the quantum phase estimation algorithm, the resource requirements are prohibitive for current quantum computers.As a more efficient alternative, several recent proposals showed that eigenvalues can be estimated, even with the Heisenberg-limited scaling, by measuring Loschmidt amplitudes at different times followed by classical signal processing [10,13,15,17].Our approach can replace the Hadamard test in these algorithms if the initial state is short-range correlated.
The computation of the local density of states in Fig. 3c of the main text follows the same idea, although we use a Fourier transform instead of more refined signal processing.This simple approach already yields valuable information about the energy eigenvalues, including the ground state energy, despite using a product state as the initial state.In the setting of quantum chemistry, the use of product states can be justified by the fact that the lowest energy Hartree-Fock has a substantial overlap with the ground state for a large class of molecules [54].In all cases, the overlap may be increased by a shallow circuit or an adiabatic evolution that maintains a short correlation length.

Temporal correlation functions
Temporal correlation functions play an important role in condensed matter physics, for instance in the study of transport and diffusion.They typically take the form where ⟨• • •⟩ represents the thermal average.Correlation functions of this type turn out to be nontrivial and hard to compute for many models even at infinite temperature [56][57][58].At infinite temperature the thermal average can be obtained by sampling random product states, with our algorithm providing the complex expectation value for each of them (If B is local, it can be absorbed into |ψ⟩ in the first protocol of section D).Finite temperatures may be accessible by combining the approach with the above spectral filtering algorithm.
Other quantities of interest are the out-of-time-order correlators (OTOCs) [70].OTOCs are widely used as a measure of information scrambling, and are typically defined as for given operators W and V .It has recently been shown that also the imaginary part of OTOCs can contain relevant information [60].Combined with the extension in

FIG. 1 .
FIG. 1. (a)The time derivative of the complex phase ϕ(t) of Loschmidt amplitude G(t) can be estimated from r(t ± ih) = | ⟨ψ ′ |e −iHt e ±hH |ψ⟩ |.The right panel shows the result of our approach with h = 0.1 for the transverse-field Ising chain, Eq. (8), of length N = 40.The solid lines on the lower right correspond to the complex Loschmidt amplitude obtained from the algorithm, while the overlapping dashed lines indicate the exact result.(b) Circuit to measure r(t ± ih).For initial product states, the rescaled imaginary-time evolution has the same brickwork layout as a single real-time Trotter step.

FIG. 2 .
FIG.2.Error in the phase of Loschmidt amplitude, ∆ϕ, computed using our approach for the transverse-field Ising chain, Eq. (8).(a) ∆ϕ/N h 2 as a function of time t with fixed real-time Trotter step τ = 0.01 for different values of the imaginary-time step h and different system sizes N .(b) ∆ϕ/N τ 2 for h = 0.01 and different values of τ and N .The color coding is the same as in (a).

FIG. 3 .
FIG. 3. (a) Absolute value of the Loschmidt amplitudes for an Ising chain of length N = 24 with initial state |↑↑↑ • • •⟩ and Trotter step sizes τ = h = 0.3.The dashed lines include single-qubit depolarizing noise with probability γ = 3 × 10 −3 after each gate.The dash-dotted lines are obtained by the error mitigation described in the text.We quantify the statistical error of the error-mitigated curves by simulating 100 experiments, each of which uses M = 10 6 measurements to estimate the survival probability.The dash-dotted line corresponds to the median of the 100 experiments, while the shaded areas indicate the range between the first and third quartile.(b) Real part of the Loschmidt amplitude computed from the data in (a) using our algorithm.The exact value under continuous time evolution is plotted for reference.The inset shows the absolute difference of the reconstructed values from the exact amplitude.(c) The LDOS obtained through discrete Fourier transform from the data in (b).The vertical, dashed line indicates the exact ground state energy E0 ≈ −7.55.

1 p 1 p /ϵ 1 pFIG. 5 .
FIG. 5.An example of phase correction when G(t) is close to 0.Here the Hamiltonian coefficients are (J, g) = (1, 1) and the initial state is still |ψ⟩ = |ψ ′ ⟩ = |↑↑ • • •⟩.The system size N = 10.The Trotter steps are τ = 0.2 and h = 0.01.(a) The magnitude r(t).The dashed lines indicate times at which r(t) almost vanishes.(b) The phase difference between the estimated and the exact value.(c-d)The real and imaginary parts of G(t).Original algorithm corresponds to our algorithm described in main text.In the "phase corrected" results, we added a discrete phase jumpy π as well as a phase shift δ, which ensures continuity of the first derivative of G(t).

TABLE I .
[50]uit depth D and number of measurements M to estimate the complex Loschmidt amplitude G with additive error ϵ.All protocols use a real-time Trotter decomposition of order p.The Hadamard test is implemented using a single ancilla qubit with swap operations in d spatial dimensions.The latter two methods require M measurements at each intermediate state or time step, but the corresponding values of the phase are also returned.For these approaches, we only consider initial product states and ϵ bounds the error r∆ϕ arising from the uncertainty in the phase.The quantities Ĩ and I depend on the intermediate amplitudes in these sequences, see text and supplemental material[50].