Toward Quantum Computing Phase Diagrams of Gauge Theories with Thermal Pure Quantum States

The phase diagram of strong interactions in nature at finite temperature and chemical potential remains largely unexplored theoretically due to inadequacy of Monte-Carlo-based computational techniques in overcoming a sign problem. Quantum computing offers a sign-problem-free approach but evaluating thermal expectation values is generally resource intensive on quantum computers. To facilitate thermodynamic studies of gauge theories, we propose a generalization of thermal-pure-quantum-state formulation of statistical mechanics applied to constrained gauge-theory dynamics, and numerically demonstrate that the phase diagram of a simple low-dimensional gauge theory is robustly determined using this approach, including mapping a chiral phase transition in the model at finite temperature and chemical potential. Quantum algorithms, resource requirements, and algorithmic and hardware error analysis are further discussed to motivate future implementations. Thermal pure quantum states, therefore, may present a suitable candidate for efficient thermal-state preparation in gauge theories in the era of quantum computing.

Digital quantum computers, as well as analog quantum simulators, offer a sign-problem-free approach to computing ground, excited, thermal, and non-equilibrium states of quantum many-body systems including lattice gauge theories (LGTs) .However, because quantum devices naturally implement pure states, quantum computing thermal, i.e., mixed, states is a challenging task.An obvious possibility, thermalizing a non-equilibrium state through unitary time evolution [69,95], may be costly on near-term quantum devices due to long equilibration times involved.Other methods for quantum computing thermal expectation values exist, which generally fall into two categories.One is sampling easily-preparable pure states according to thermal statistics [96][97][98][99][100][101][102], but poor scaling of required samples with system size can lead to prohibitive resource requirements.Another is preparing full (mixed) thermal states on the quantum computer, often relying on resource-intensive algorithms such as quantum phase estimation and incurring significant ancillaryqubit requirements [103][104][105][106][107][108].Hybrid classical-quantum algorithms for near-term applications [97], such as those relying on variational methods [106,[109][110][111][112][113][114][115] or utilizing Jarzynski's equality [116][117][118] are also explored in recent years.
A promising near-term perspective is to resort to the thermal-pure-quantum-(TPQ-) state formulation of statistical mechanics [120][121][122] to obtain thermal expectation values using only a single properly-prepared pure state in the thermodynamic limit.This approach offers a promising path to simulating quantum systems at finite temperature and chemical potential on quantum computers [123,124].Canonical TPQ states, in particular, are formed by imaginary-time evolving pure states drawn from an ensemble of Haar-random states.Thermal expectation values of mechanical observables with TPQ states converge rapidly to the values within the standard (ensemble) formulation of statistical mechanics, improving exponentially with system size.In fact, as is shown in Ref. [120], for sufficiently large systems often a single TPQ state suffices.However, in gauge theories, local symmetries put constraints on the allowed Hilbert space so that the naïve TPQ-state approach cannot be applied.In this manuscript, a protocol for explicitly constructing TPQ states in the physical subspace of gauge theories will arXiv:2208.13112v1[hep-lat] 28 Aug 2022 x,y < l a t e x i t s h a 1 _ b a s e 6 4 = " F x W M G o 9 L S S 7 5 u B 9 K i A n d s Q B x 4 B E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i z c r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g q e N U M W y w W M S q H V C N g k t s G G 4 E t h O F N A o E t o L R z d R v P a H S P J a P Z p y g H 9 G B 5 C F n 1 F j p 4 a 4 n e + W K W 3 V n I M v E y 0 k F c t R 7 5 a 9 u P 2 Z p h N I w Q b X u e G 5 i / I w q w 5 n A S a m b a k w o G 9 E B d i y V N E L t Z 7 N T J + T E K n 0 S x s q W N G S m / p 7 I a K T 1 O A p s Z 0 T N U C 9 6 U / E / r 5 O a 8 M r P u E x S g 5 L N F 4 W p I C Y m 0 7 9 J n y t k R o w t o U x x e y t h Q 6 o o M z a d k g 3 B W 3 x 5 m T T P q t 5 F 1 b s / r 9 S u 8 z i K c A T H c A o e X E I N b q E O D W A w g G d 4 h T d H O C / O u / M x b y 0 4 + c w h / I H z + Q M r p o 2 6 < / l a t e x i t > O n < l a t e x i t s h a 1 _ b a s e 6 4 = " H T Z C 3 w 8 R S y 3 5 4 A e d t h L t K S R i r t 4 = " > A A A B 7 X i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v o q e y K q M e i F 2 9 W s B / Q L i W b Z t v Y b L I k W a E s / Q 9 e P C j i 1 f / j z X 9 j 2 u 5 B q w 8 G H u / N M D M v T A Q 3 1 v O + U G F p e W V 1 r b h e 2 t j c 2 t 4 p 7 + 4 1 j U o 1 Z Q 2 q h N L t k B g m u G Q N y 6 1 g 7 U Q z E o e C t c L R 9 d R v P T J t u J L 3 d p y w I C Y D y S N O i X V S 8 7 a X y e N J r 1 z x q t 4 M + C / x c 1 K B H P V e + b P b V z S N m b R U E G M 6 v p f Y I C P a c i r Y p N R N D U s I H Z E B 6 z g q S c x M k M 2 u n e A j p / R x p L Q r a f F M / T m R k d i Y c R y 6 z p j Y o V n 0 p u J / X i e 1 0 W W Q c Z m k l k k 6 X x S l A l u F p 6 / j P t e M W j F < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 f s 7 Ramsey Interferometry + %&'( < l a t e x i t s h a 1 _ b a s e 6 4 = " X D j P v r b u q r w x 4 k u

Ramsey Interferometry
< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 f s 7 Ramsey Interferometry Ramsey Interferometry Ramsey Interferometry Ramsey Interferometry

One Layer
FIG. 1. Schematic overview of the PTPQ-state preparation.The colored squares in the R circuit denote randomly-chosen single-qubit rotations around X and Y axes of the Bloch sphere by angle π 2 and the T gates.The entangling gates are controlled-Z rotations.The two consecutive single-qubit gates on each qubit are constrained to be different in this construction.Entangling gates in each layer only act on adjacent qubits and this pattern continues, see Ref. [119] for details.H is a Hadamard gate.
be introduced by penalizing non-physical components of the random pure state as they evolve in imaginary time.
By accurately obtaining, via a numerical simulation, the phase diagram of a simple gauge theory (Z 2 LGT in 1 + 1 D with matter) at finite temperature and chemical potential, we demonstrate the utility of TPQ-state approach in studying thermodynamics of gauge theories for the first time.Aiming at quantum-computing applications, associated quantum algorithms, quantum-resource requirements, and robustness to algorithmic and hardware errors are further studied.The results indicate that the TPQ-state approach may be a suitable candidate for efficient phase-diagram studies of QCD in the future.
Thermal Pure Quantum States for Gauge Theories.Canonical TPQ states are defined as [120] with β being the inverse temperature, N the number of degrees of freedom of the (discrete) system, and ψ R ⟩ a Haar-random state.TPQ states approximate thermal expectation values of 'mechanical' operators, i.e., those that are low-degree polynomials of local operators, via with exponential convergence in the system size (as well as in inverse temperature), see Supplemental Material [125] for details.While a single TPQ state suffices as N → ∞, for faster convergence at finite N , a stochastic average over r TPQ states can be preformed, denoted by ⟪⋅⟫ r in the formula.
In gauge theories, Φ R ⟩ may be unphysical, in which case Eq. (1) will not reproduce physical thermal observables.While this issue can be avoided by eliminating the gauge-field degrees of freedom with certain boundary conditions in 1+1 D, such a strategy is not generally applicable.Therefore, we propose 'physical' thermal pure quantum (PTPQ) states by adding a term to the Hamiltonian, where G n are Gauss's law operators at site n, that is [H, G n ] = 0.The function f is chosen such that unphysical components of the state as it evolves in imaginary time are penalized in energy.Such an approach is customary in the context of enforcing Gauss's law in analog and digital quantum simulation of gauge theories, and can be applied to both Abelian and non-Abelian cases [46,84,[126][127][128][129][130][131].
A circuit to prepare PTPQ states on quantum computers is illustrated in Fig. 1.First, a random circuit R, consisting of layers of single-qubit gates and entangling two-qubit gates, is used to prepare an approximate Haarrandom state.Various designs are suggested for such task with studied performance [132], and we adopt the efficient implementation of Ref. [119].This random circuit is followed by a non-unitary operator e −βH 2 acting upon the resulting random state to produce a standard canonical TPQ state.Gauss's law is enforced through action with Q G ≡ e − β 2 ∑ n f (Gn) , the circuit implementation of which depends on the f chosen, see below for the example of Z 2 LGT in 1+1 D. These elements will be further studied in the following.

Thermal chiral phase diagram of Z 1+1
2 .The model that will be studied in the following to demonstrate the value of the TPQ-state approach in gauge theories is Z 2 LGT in 1+1 D coupled to staggered fermions (Z 1+1 2 ).This model is sufficiently simple to allow numerical verifications on classical computers, while it still exhibits a nontrivial phase diagram which is aimed to be reproduced by quantum simulation.The Hamiltonian of the model is The order parameter to be evaluated is the chiral condensate, with ⟨ ΨΨ⟩ ≠ 0 (= 0) in the chiral-symmetry broken (symmetric) phase.In QCD, chiral symmetry is spontaneously broken at small T [133] but restored at large T [134 -138].Chiral symmetry breaking is also observed in 1+1 D Abelian gauge theories [139][140][141].For Z 1+1 2 with staggered fermions, chiral symmetry restoration can be understood intuitively.At asymptotic µ ≫ m, , the state with all fermion sites occupied dominates the thermal Gibbs ensemble and has ⟨ ΨΨ⟩ = 0 according to Eq. ( 5), while at T → ∞, the state is an equal admixture of all basis states, giving rise to a net vanishing condensate.
Figure 2 depicts the phase diagram using PTPQ states with random-circuit depth d = 20 and averaged over r = 10 PTPQ-state realizations.All results are classically computed using exact diagonalization on a system of N = 6 sites corresponding to 11 qubits [125].Chiral symmetry is observed to be restored at high T and µ and is broken otherwise.Despite the small system size, the transition at T → 0 is visibly discontinuous, indicating that a 'true' phase transition is expected in the infinitevolume limit, while only a smooth crossover is observed at µ → 0. Top and right panels of Fig. 2 show a comparison between exact results and PTPQ states.Reasonable agreement is observed, and the accuracy is expected to improve with increasing system size, see Supplemental Material for details [125].Furthermore, as demonstrated in Supplemental Material [125], the standard TPQ states do no reproduce the correct phase diagram unless promoted to PTPQ states.In contrast, regardless of the use of PTPQ states, detecting the confinement-deconfinement transition, where one computes non-local observables, e.g., string tension, is more challenging on small systems due to finite-size effects, as discussed in Supplemental Material [125].
Thermal non-equal-time correlation functions can also be evaluated using the PTPQ-state protocol.These quantities are challenging to access with classical Monte-Carlo techniques even at vanishing chemical potential.As an example, we focus on the thermal current-current correlator, where j n (t) ≡ e iHt j n e −iHt , is the fermionic current operator, and ⟨⋅⟩ β denotes the thermal expectation value at finite β.A quantum circuit for evaluating C n (t) is presented in Fig. 1(c) and involves a standard Ramsey interferometry scheme [144][145][146][147], requiring

Imaginary Time Evolution Random-State Preparation
Ramsey Interferometry + %&'( < l a t e x i t s h a 1 _ b a s e 6 4 = " X D j P v r b u q r w x 4 k u

y t t S I U P C Q x Z p y E t j M h O D S L l T z + t m G F H O Z d p h k z S + a I E y q d / q G H N K I q x J Y R q b m Z B o Q t E G V L Y h + I s v L P W e c / r P n F X T R F H C Y h B M A h y u o w x o A k U H u E Z X u H N U c L + z F t X n G L m C P A + f w B s E + P
Ramsey Interferometry + %&'( < l a t e x i t s h a 1 _ b a s e 6 4 = " X D j P v r b u q r w x 4 k u    an ancilla in a Hadamard superposition controlling j n (t) and j 0 (0), with the PTPQ-state input.The results of this procedure (classically computed) are shown in Fig. 3 for N = 6, T m = 0.2, and µ m = 0.5.The exact results, obtained with exact diagonalization, are compared with a PTPQ computation using r = 1 and r = 20 samples for lower and higher temperatures, respectively, and good agreement is observed.Since real-time evolution is natural for a quantum computer, PTPQ states are expected to be an important tool to compute e.g., transport quantities like conductivity or viscosity.For example, one can use Eq. ( 6) to compute thermal conductivity, σ β,µ ∼ lim ω→0 a ω ∫ dt e iωt ∑ n C n (t), which nonetheless requires long evolution times and large system sizes to be obtained accurately.

Robustness to Quantum-Circuit Implementation Errors.
To investigate inaccuracies occurring in the quantumcircuit implementation of the PTPQ-state preparation in Z 1+1 2 , we recall that this preparation consists of two parts, generation of a Haar-random state, followed by imaginary-time evolution with First, by employing an exemplary current noise model, the effect of device noise in generating a Haar-random state can be investigated using the circuit described in Fig. 1(b).Shown in Fig. 4(a) is the chiral transition of an ideal device versus IBM's 'Sydney' device [142] (simulator) when noise is present in the Haar-random-state preparation, while the imaginary-time evolution is per-formed exactly (shot noise is not considered).Data corresponds to a circuit depth of d = 20 and r = 10 PTPQ realizations for a system of N = 6 fermionic sites.As is seen, device errors included in current simulators have a negligible effect during this part of the algorithm.Whether features like phase transitions remain robust with respect to the actual device error will need to be tested on the available quantum hardware in future studies.
The second part of the algorithm is imaginary-time evolution.Here, an error analysis will specifically depend on the chosen implementation strategy and device.Among the common near-term approaches include the Quantum Imaginary Time Evolution (QITE) algorithm [98] and dilated operator approaches [148,149].Both approaches employ Trotterization for imaginarytime evolution, hence it is beneficial to investigate the Trotterization error in the PTPQ-state preparation.
For simplicity, we consider Trotterization via a firstorder product formula of the form with ∆β ≡ β 2N T , where H = ∑ Γ γ=1 H γ with H γ being each non-commuting term in the Z 1+1 2 Hamiltonian in Eq. ( 4) plus the chemical-potential contribution, after Jordan-Wigner transforming fermionic into spin degrees of freedom.Note that implementing Q G does not introduce any Trotter error as [H n , Q G ] = 0.The Trotter-step (N T ) dependence of the chiral condensate ⟨ ΨΨ⟩ as a function of chemical potential µ at a relatively long imaginary time corresponding to temperature T m = 0.1 is plotted in Fig. 4(b).As is seen, with a small number of steps N T = 10 (∆β = 1), the chiral phase transition is reproduced qualitatively.To achieve percent-level agreement, N T ≳ 100 is required corresponding to step sizes ∆β ≲ 0.1.This suggest that despite large Trotter errors, qualitative features of the phase diagram, such as phase transitions, can be still accessible.
An analytic bound for the multiplicative Trotter error, M (∆β) = S(∆β) e −∆β H − I, has been obtained in Ref. [143], Here, αcomm ≡ ∑ Γ γ1,γ2=1 [H γ2 , H γ1 ] and ⋅ is the spectral norm.The inset of Fig. 4(b) shows the multiplicative Trotter error as a function of ∆β, along with the numerically determined M (∆β) and the mean trace distance DT r between e −∆βH Ψ R ⟩ and S(∆β) Ψ R ⟩ averaged over r = 10 random-state realizations.As can be seen, the numerical results are markedly different from the predicted exponential analytic bound, indicating a power-law dependence on ∆β instead.In fact, the multiplicative Trotter error is seen to scale as ≈ (∆β) 2 , similar to the scaling of the bound derived in Ref. [143] for Trotterized realtime evolution.Indeed, for e −∆βH with ∆β > 0 and a physical Hamiltonian H with lower-bounded spectrum, the exponential contribution in Eq. ( 8) must be replaced by a bounded function to provide a tighter bound.
The influence of device errors on the imaginary-time evolution, through imperfect gates and decoherence, is hardware and algorithm specific and thus difficult to estimate.Nonetheless, one can consider a deviceindependent approach, in which the errors can be parameterized as imaginary-time evolution under an effective Hamiltonian [150] err , with H α err representing 1-or 2-local errors with randomized weights bounded by α.Because both fermion (via a Jordan-Wigner transformation) and gauge variables are spin operators, they will be denoted from now on with the same symbol σ l , where l denotes sites in a chain of length 2N − 1, and even (odd) sites labeling fermions (gauge fields).Motivated by the common bit flip, phase flip, and crosstalk errors in current hardware, we choose Denoting ΨΨ⟩ P T P Q (α) the chiral condensate calculated with the PTPQ states under such a noise model, a relative accuracy can be defined as The chiral condensate under this noise model is plotted in Fig. 4(c) as a function of chemical potential, at T m = 0.1 for α m = 0.1, 0.25, 1, with exact values displayed for reference.For σ x l errors, PTPQ states are still able to capture phase transitions for α m = 0.1, but this ability decays as α m → 1, as the relative accuracy near the transition quickly drops to ∼ 45%.Furthermore, the inset of Fig. 4(c) shows A rel for µ m = 1.5 as a function of α m.A strong difference is observed between σ z l and σ x l errors versus σ z l σ z l+1 errors.σx n and σ z n are both gauge invariant, so the introduced X-type and Z-type noise only violate Gauss's law on half of the qubits.Conversely, ZZtype noise always violates Gauss's law since each term acts on one link and one site.Since the penalty term Q G effectively protects against Gauss's-law-violating errors, PTPQ states are seen to be more robust to the latter.In Supplemental Material [125], we discuss the relation between the link-operator convention and the unitary-noise mitigation.
Scaling and Resource Requirements.While an exact preparation of Haar-random state has a time complexity that scales exponentially in the system size, pseudo-Haarrandom states that are sufficient for most applications can be prepared in polynomial time, as has been established in Refs.[119,[151][152][153].For example, the (pseudo-) random-state preparation circuit for PTPQ states in this work requires only O(N d) two-qubit entangling gates, where N is the number of qubits and d is number of layers.In practice, d ≃ 20 is seen to be sufficient in the examples of this work.In general, the required d for given precision depends on the system size but the scaling is observed to be sublinear [119,153,154].
To cost the imaginary-time evolution, we consider the QITE algorithm [98] because of its prominence among near-term approaches.The QITE algorithm consists of Trotterizing imaginary-time evolution, e.g., Eq. ( 7), and it proceeds to find, using classical processing or variationally, a unitary operator e −i∆βAγ that approximates the non-unitary operator e −Hγ ∆β , such that for given intermediate state Φ⟩ in the Trotter sequence, e −i∆βAγ Φ⟩ ≈ e −Hγ ∆β Φ⟩ N β , where N β ≡ ⟨Φ e −2Hγ ∆β Φ⟩.Note that for PTPQ-state preparation, Q G does not need to be Trotterized and can be performed in a single QITE step for time duration β.
Following this protocol, one incurs two types of errors: the Trotter error ε T ≥ M (∆β) , discussed in the previous section, and the intrinsic QITE error ε Q from approximating a non-unitary evolution with a unitary evolution after N T applications.Here, A γ acts on N q qubits.Increasing N q reduces ε Q but increases the A γ circuit cost, i.e., the total algorithm runtime of the QITE algorithm is ΓN T e O(Nq) [98,125].To estimate this dependence, assume N T ≥ βη where H γ ≤ η is a bound on the operator norm of H γ .Then N q = O 2kC ln 2 √ 2ΓN T ε −1 Q , with k the locality of the Hamiltonian and C the typical correlation length of the state which is acted upon.The Trotter-step dependence can be exchanged with Trotter error ε T using Eq. ( 8), but it should be noted that a more favorable behavior is observed empirically, compared to the predicted exponential step-size dependence, see Fig. 4(b).We add that alternatively, a hybrid interaction-picture QITE-Qubitization simulation proto-col can be employed [65] to improve the scaling when the Hamiltonian terms have large spectral norms, such as Z 1+1 2 in presence of Gauss's law enforcing penalty term.Note that in essence, the QITE circuit complexity is exponential in the correlation length C. While the PTPQ algorithm starts from an infinite-temperature state, where C = 0, this will become prohibitive at a phase transition (of an infinite system), thus highlighting the importance of developing efficient far-term non-unitary evolution schemes.Further details on on the algorithmic scaling is provided in Supplemental Material [125].Conclusions.Physical thermal pure quantum states, introduced in this work as an extension of (canonical) thermal pure quantum states, present a valuable method for quantum computing thermal phase diagram and thermal non-equal-time correlation functions in strongly interacting gauge theories.Via the example of a 1+1 dimensional Z 2 lattice gauge theory coupled to staggered fermions on a small lattice, PTPQ states are shown numerically to reproduce the phase diagram very accurately over all values of temperature and chemical potential, including a chiral phase transition, with favorable resource requirements.Phase transitions without a local order parameter, such as the (de-)confinement transition present in this model, are strongly finite-size dependent, and will be accessible to future quantum computers.Similarly, thermal no-equal-time correlation functions, important for obtaining transport properties of the Quark-Gluon Plasma in ultra-relativistic heavy ion collisions [155] can be accessed using the algorithms of this work on the upcoming quantum-computing devices.
Consider a mechanical observable that can be formally defined by the following bound on its operator norm: O ≤ KN , where K is a constant independent of O and N , and ≪ N denotes the degree of the polynomial of local operators constituting O. N denotes the size of the (discrete) system and here for generality, we do not differentiate between the number of lattice sites and the number of qubits encoding system's degrees of freedom.A Thermal Pure Quantum (TPQ) state has the property that the expectation values of all mechanical observables calculated using a TPQ state converge uni-formly to those obtained from the ensemble formulation of statistical mechanic in the thermodynamic limit.One can then proceed to find properly-prepared states such that they satisfy this condition.
An ansatz for TPQ states is proposed in Ref. [120], where ψ R ⟩ = ∑ i c i i⟩ with {c i } being complex numbers drawn randomly from a unit sphere ∑ i c i 2 = 1, and { i⟩} is any arbitrary orthogonal basis states spanning the full Hilbert space, e.g., a set of product states in the computational basis.To show that this ansatz indeed has the property of a TPQ state, one notes that: i) The thermal expectation value in the ensemble formulation is ⟨O⟩ β,N = Tr(e −βH O) Tr(e −βH ), where trace is defined over any complete orthogonal set of basis states spanning the Hilbert space of a system of size N .One then has the identity in the limit of a sufficiently large number, r, of random states.Here, ⟪⋅⟫ r denotes statistical average over r TPQ states, and properties of random states are used to relate the statistical average of expectation values to trace operation [156].
ii) The quantity D 2 N , defined by can be shown to be bounded as D 2 N ≤ N 2 e β Θ(N ) , see Ref. [120] for a derivation.Here, Θ(N ) indicates that the form is bounded from above and below by k 1 N and k 2 N for two constants k 1 and k 2 at large N .
iii) Finally, generalized Markov's inequality puts a bound on the following probability for any arbitrary > 0. As is seen, this bound exponentially vanishes in the limit of N → ∞, indicating that only r = 1 random-state realization is sufficient to reproduce, with a high probability, the true ensemble expectation values for large N .This completes the proof that the state defined in Eq. (S1) is a TPQ state.Note that for a fixed system size, the higher the temperature, the larger the bound on D 2 N , hence requiring more random-state realizations to achieve a given accuracy on thermal expectation values.This feature is evident from the different number of random states used to reproduce, with comparable accuracy, the non-equal-time correlation functions in Fig. 3 of the main text for two different temperatures.

B. Comparison between Standard and Physical TPQ States
The need for physical TPQ (PTPQ) states in the context of calculating thermal expectation values of LGTs can be demonstrated by evaluating the chiral condensate of Z 1+1 2 both with the standard and the physical TPQ states.The result is shown in Fig. S1 for T m = 0.1, along with the exact values for reference.As in the main text, all results are presented for m = 1 2, m = 1 and a = 1.It is clear that standard TPQ states are not suitable for computing finite-temperature phase diagrams of LGTs, since the errors incurred by drawing the requisite approximate Haar-random state from the entire Hilbert space can lead to incorrect results for the order parameter.This effect is seen to be especially pertinent near the phase transitions.The correct behavior is reproduced by utilizing PTPQ states instead.

C. Chiral Condensate with PTPQ States vs. Exact Values: Convergence with Lattice Size
As described in Sec.A, the quantity D 2 Nq defined in Eq. (S3) quantifies the convergence of thermal expectation values evaluated with (P)TPQ states to the values obtained within the ensemble formulation of statistical mechanic, and is expected to be bounded by a function that decreases exponentially quickly in the number of qubits N q .This quantity for chiral condensate is plotted in Fig. S2   Determining the finite-temperature confinementdeconfinment phase transition and the associated phase diagram is more challenging due to thermal and finitesize effects, however PTPQ states could still be used to explore these phases with sufficient computational resources.At T = 0, confinement can be detected by the decay of the non-local two-point correlation function of two gauge-invariant fermionic operators as a function of distance i among them.Explicitly, an exponential decay of this two-point correlator indicates confinement [157].As depicted in Fig. S3(a), the value of ⟨f † 0 f i ⟩ as a function of i decays exponentially at ≠ 0, where there is expected to be confinement [157].However, a finite temperature can introduce an exponential decay in this correlator even in the absence of confinement, and more robust diagnostics of confinement are required.Additionally, the non-local correlator introduced above constitute a large-order polynomials of local operators, which is prohibitive to the (P)TPQ-state method.
A better suited quantity for exploring confinement at finite temperatures is string tension, which can be expressed as: Here, E(E * ) is the ground-state energy of the system without (with) a pair of static test charges each placed at one end of the lattice, and N is the number of lattice sites.However, finite-size effects appear to be significant at lattice sizes accessible to our classical exact diagonalization, hindering the identification of confined and deconfined phases.For example, Fig. S3(b) plots string tension for a system of N = 6 fermionic sites as a function of temperature at µ m = 0.2.In the infinite-temperature limit, where individual local states are purely disordered and ⟨σ x ⟩ = ⟨σ y ⟩ = ⟨σ z ⟩ = 0 for all n, the system approaches zero energy (dropping constant contributions).Since σ 2 , the addition of each test charge yields a mass energy gain of m 2 , and the test charge at the first lattice site yields an electric field energy gain of .So, in the infinite temperature limit, the string tension approaches m+ N .This vanishes in the bulk limit, indicating a deconfined phase at high temperatures.While the numerical values supports convergence onto this value, clearly identifiable regions of confinement and deconfinement outside of this regime would require simulating larger lattices.Nonetheless, TPQ states can still be used to reproduce the expected behavior.Data points corresponding to the string tension evaluated using the PTPQ states are shown in the same plot, indicating good agreement with the exact values.Therefore, PTPQ-state methodology may be used to construct the confinementdeconfinement phase diagram of Z 1+1 2 and other LGTs on future quantum computers.

E. Link-Operator Convention and Robustness to Unitary Errors
An alternative convention may be used for the link operators such that the electric-field operator is instead represented by σz n , which switches the roles of σx n and σz n in the physical Hamiltonian.While this does not affect exact results, the choice of convention does significantly impact how robust PTPQ states are to different types of introduced unitary error.As was seen in the main text, PTPQ states are robust to error terms that violate Gauss's law, as the penalty term in the PTPQ-state preparation effectively suppresses these errors as well.While the common unitary errors in the quantum device are fixed in the computational basis, the operators that violate Gauss's law in each convention are different, hence one expects the most prominent error types to be different for each convention.
In particular, in the alternate convention to that in the main text, σ x n and σx n both violate Gauss's law, so this type of introduced error is effectively mitigated by Q G .However, neither σ z n nor σz n violate Gauss's law, so PTPQ states are more sensitive to Z and ZZ-type noise.Fig. S4 demonstrates this change of sensitivity to the same types of introduced unitary errors as those examined in Fig. 4(c) of the main text.Therefore, the freedom to choose a link-operator convention can be used to increase PTPQ-state-preparation fidelity depending on the most prominent errors seen in a given hardware.

F. Comments on the QITE-algorithm Runtime Scaling
The runtime analysis of the Quantum Imaginary-Time Evolution (QITE) algorithm has been performed in Ref. [98].The derivation assumes that the norm of each Hamiltonian term is bounded by one-an assumption that can be easily removed upon a slight modification of the results.Let us start by the reverse triangle inequality, e −∆βHγ ≥ I − ∆βH γ ≥ 1 − ∆βη, (S7) where H γ ≤ η for η > 0. Thus, for N T = β 2 ∆β ≥ βη [158], e −∆βHγ ≥ 1 2, matching Eq. (30) of Ref. [98].Then, following the analysis in Ref. [98] with N T ≥ βη, the total algorithm runtime becomes Here, Γ is the number of terms in the Hamiltonian with a given decomposition H = ∑ Γ γ=1 H γ , where each Hamiltonian term acts on at most k neighboring qubits on the underlying interaction-to-qubit graph, C is the upper bound on the correlation length of Φ⟩ (in terms of qubits) [159], where Φ⟩ is a given intermediate state in the Trotter sequence as explained in the main text.Furthermore, d is the spatial dimensionality of the Hamiltonian, and ε Q is the error due to approximating non-unitary evolution with a unitary evolution after N T steps.The total runtime is thus seen to scale quasipolynomially with η.
Finally, an improvement to the algorithmic scaling of QITE in Ref. [98] is to incorporate the improved productformula bounds, for example that derived in Ref. [143] combined with the empirical exponential improvement observed in the main text for the physical Hamiltonian of this work.For Hamiltonians of the form H = ∑ Γ γ=1 H γ and with first-order product formulas, it follows from Eq. ( 43) of Ref. [143] that the Trotter error of timestep ∆β is bounded as in Eq. ( 8) of the main text.Nonetheless, the exact evaluation of the (multiplicative) error for small system sizes reveals that the error in each Trotter step goes as O((∆β) 2 ), hence the number of Trotter steps to reach an overall accuracy ε T after evolving for β 2 = N T ∆β in (P)TPQ-state preparation goes as N T = O(β 2 ε T ).Merging this with the QITE resource requirements discussed in the main text and in Ref. [98], provided that N T ≥ βη, one concludes that each unitary solved for within the QITE algorithm must act on at most N q number of qubits as given in the main text.
n H f n Y 9 w 6 5 0 x m d t A f O J 8 / e i e Z I w = = < / l a t e x i t > | , Ni phys < l a t e x i t s h a 1 _ b a s e 6 4 = " I T m a s p l S Q 4

ΨΨFIG. 2 .
FIG. 2. Chiral phase diagram for the Z2 LGT with fermions in 1 + 1 D with N = 6 and λ m = 4, calculated with PTPQ states with random-circuit depth d = 20 and averaged over r = 10 PTPQ-state realizations.Error bars denote statistical uncertainty from a finite PTPQ sample.
< l a t e x i t s h a 1 _ b a s e 6 4 = " S p 9 8 1 Q e S / 2 9 c f M + O y n 1 E h w i H G e M = " >A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R g x 5 b t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 Cd 4 8 a C I V 3 + R N / + N 2 z Y H r T 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 c g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g p e N U M W y y W M S q E 1 C N g k t s G m 4 E d h K F N A o E t o P x z c x v P 6 L S P J Y P Z p K g H 9 G h 5 C F n 1 F j p v t G / 7 Z c r b t W d g / w l X k 4 q k K P e L 3 / 2 B j F L I 5 S G C a p 1 1 3 e l P x P 6 + T m v D a z 7 h M U o O S z R e F q S A m J t O v S Z 8 r Z E a M L a F M c X s r Y U O q K D M 2 m 6 I N w V t 8 e Z k 0 z y v e Z c W r X 5 S r N 3 k c B T i G E z g D D 6 6 g C n d Q g w Y w Q H i G V 3 h z H p 0 X 5 9 3 5 m L e u O P n M E f y B 8 / k D r p m M 3 A = = < / l a t e x i t > R < l a t e x i t s h a 1 _ b a s e 6 4 = " e t x f F I A o b H k 9 b A b w n Y G z O M h D e n o = " > A A A B 7 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X Y B E 8 l U R E P R a 9 e K x g P 6 A N Z b P Z t G s 3 u 2 F 3 I p T Q / + D F g y J e / T / e / D d u 2 x y 0 9 c H A 4 7 0 Z Z u a F q e A G P e / b W V l d W 9 / Y L G 2 V t 3 d 2 9 / Y r B 4 c t o z J NW Z M q o X Q n J I Y J L l k T O Q r W S T U j S S h Y O x z d T v 3 2 E 9 O G K / m A 4 5 Q F C R l I H n N K 0 E q t H o 0 U m n 6 l 6 t W 8 G d x l 4 h e k C g U a / c p X L 1 I 0 S 5 h E K o g x X d 9 L M c i J R k 4 F m 5 R 7 m W E p o S M y Y F 1 L J U m Y C f L Z t R P 3 1 C q R Gy t t S 6 I 7 U 3 9 P 5 C Q x Z p y E t j M h O D S L 3 l T 8 z + t m G F 8 H O Z d p h k z S + a I 4 E y 4 q d / q 6 G 3 H N K I q x J Y R q b m 9 1 6 Z B o Q t E G V L Y h + I s v L 5 P W e c 2 / r P n 3 F 9 X 6 T R F H C Y 7 h B M 7 A h y u o w x 0 0 o A k U H u E Z X u H N U c 6 L 8 + 5 8 z F t X n G L m C P 7 A + f w B s E + P N Q = = < / l a t e x i t >

<
l a t e x i t s h a 1 _ b a s e 6 4 = " S p 9 8 1 Q e S / 2 9 c f M + O y n 1 E h w i H G e M = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R g x 5 b t B / Q h r L Z T t q l m 0 3 Y 3 Q g l 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H r T 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 c g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g p e N U M W y y W M S q E 1 C N g k t s G m 4 E d h K F N A o E t o P x z c x v P 6 L S P J Y P Z p K g H 9 G h 5 C F n 1 F j p v t G / 7 Z c r b t W d g / w l X k 4 q k K P e L 3 / 2 B j F L I 5 S G C a p 1 1 3 H p 0 X 5 9 3 5 m L e u O P n M E f y B 8 / k D r p m M 3 A = = < / l a t e x i t > R < l a t e x i t s h a 1 _ b a s e 6 4 = " e t x f F I A o b H k 9 b A b w n Y G z O M h D e n o = " > A A A B 7 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X Y B E 8 l U R E P R a 9 e K x g P 6 A N Z b P Z t G s 3 u 2 F 3 I p T Q / + D F g y J e / T / e / D d u 2 x y 0 9 c H A 4 7 0 Z Z u a F q e A G P e / b W V l d W 9 / Y L G 2 V t 3 d 2 9 / Y r B 4 c t o z L 5 P W e c 2 / r P n 3 F 9 X 6 T R F H C Y 7 h B M 7 A h y u o w x 0 0 o A k U H u E Z X u H N U c 6 L 8 + 5 8 z F t X n G L m C P 7 A + f w B s E + P N Q = = < / l a t e x i t > • • • < l a t e x i t s h a 1 _ b a s e 6 4 = " e t x f F I A o b H k 9 b A b w n Y G z O M h D e n o = " > A A A B 7 X i c b V B N S 8 N A E J 3 4 W e t X 1 a O X Y B E 8 l U R E P R a 9 e K x g P 6 A N Z b P Z t G s 3 u 2 F 3 I p T Q / + D F g y J e / T / e / D d u 2 x y 0 9 c H A 4 7 0 Z Z u a F q e A G P e / b W V l d W 9 / Y L G 2 V t 3 d 2 9 / Y r B 4 c t o z

FIG. 4 .
FIG. 4. (a) Chiral condensate ⟨ ΨΨ⟩ as a function of µ m, including simulated device noise [142] on the Haar-random statepreparation circuit shown in Fig. 1(b).(b) Trotterization error in evaluating ⟨ ΨΨ⟩ as a function of µ m for various numbers of Trotter steps N T .Inset: First-order product-formula Trotter error bound from Ref. [143] as a function of ∆β for N = 6 at µ m = 4. Explicitly calculated multiplicative Trotter error M (∆β) is plotted for comparison, as well as the mean trace distance DT r between e −∆βH Ψ R ⟩ and S(∆β) Ψ R ⟩. (c) Unitary errors in ⟨ ΨΨ⟩ as a function of µ m calculated with PTPQ states with introduced unitary Z errors at various relative Hamiltonian error strengths α m.Inset: Relative accuracy A rel as a function of relative Hamiltonian error-term strength α m at µ m = 3 when various types of unitary error are added to the imaginary-time evolution.All results are obtained for d = 20, r = 10, λ m = 8, T m = 0.1, and averaged over results from 50 randomly generated error terms at each value of α m.
FIG. S1.Chiral condensate ⟨ ΨΨ⟩ at T m = 0.1 as a function of (scaled) chemical potential µ m, calculated with standard TPQ states and PTPQ states.Exact results are also plotted for reference.All results are obtained for N = 6 fermionic sites, λ m = 4, r = 10 random-state realizations, and a random-circuit depth d = 20.
as a function of N q .Data from simulated systems up to N = 8, at T m = 0.1, µ m = 0.2, and λ m = 100, averaged over r = 50 PTPQ states with a random-circuit depth d = 50, show an exponential decay in D 2Nq consistent with the expectation.
FIG. S2.Chiral-condensate error convergence with respect to the number of qubits Nq.Results correspond to T m = 0.1, µ m = 0.2, and λ m = 100, averaged over r = 50 PTPQ states with a random-circuit depth d = 50.
FIG. S3.(a) Exponential decay of the ground-state two-point fermionic correlator ⟨f † 0 fi⟩ as a function of i for a system of N = 12 fermionic sites (the Gauss's law is enforced a priori hence allowing to numerically simulate larger lattices), indicating a confined phase at T = 0 and µ m = 0.2.(b) Finitetemperature string tension for a lattice of N = 6 fermionic sites and at µ m = 0.2, calculated using physical TPQ states with λ m = 30, d = 20, and r = 50.Exact results are plotted for reference.
FIG. S4.Relative accuracy A rel defined in Eq. (9) of the main text as a function of relative Hamiltonian error-term strength α m when various types of unitary error are added to the imaginary-time evolution.All results are obtained for d = 20, r = 10, λ m = 8, T m = 0.1, µ m = 3, and averaged over results from 50 randomly generated error terms at each value of α m. Results shown use the alternative convention for link operators, where electric-field operators are represented by σz l