# Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy for quantum error mitigation

Theo Saporiti, Oleg Kaikov, Vasily Sazonov, Mohamed Tamaazousti

► To cite this version:

Theo Saporiti, Oleg Kaikov, Vasily Sazonov, Mohamed Tamaazousti. Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy for quantum error mitigation. Physical Review A, 2025, 112 (3), pp.032409. 10.1103/slbz-3r4b. hal-04989400

HAL Id: hal-04989400

https://hal.science/hal-04989400v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy for quantum error mitigation

Theo Saporiti,*, Oleg Kaikov, Vasily Sazonov, and Mohamed Tamaazousti

Université Paris-Saclay, CEA, List, F-91120 Palaiseau, France

<figure><img src="image_2.png" /></figure>

(Received 13 March 2025; accepted 29 July 2025; published 3 September 2025)

Mitigation of quantum errors is critical for current NISQ devices. In the present work, we address this task by treating the execution of quantum algorithms as the time evolution of an idealized physical system. We use knowledge of its physics to assist the mitigation of the quantum noise produced on the real device. In particular, the time evolution of the idealized system obeys a corresponding Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of equations. This is the basis for the novel error mitigation scheme that we propose. Specifically, we employ a subset of the BBGKY hierarchy as supplementary constraints in the zero-noise extrapolation (ZNE) method for error mitigation. We ensure that the computational cost of the scheme scales polynomially with the system size. We test our method on digital quantum simulations of the lattice Schwinger model under noise levels mimicking realistic quantum hardware. We demonstrate that our scheme systematically improves the error mitigation for the measurements of the particle number and the charge within this system. Relative to ZNE, we obtain an average reduction of the error by (18.2 ± 0.5)% and (52.8 ± 6.3)% for the respective above observables. We propose further applications of the BBGKY hierarchy for quantum error mitigation.

DOI: 10.1103/slbz-3r4b

## I. INTRODUCTION

The coupling of quantum computers to their surrounding environment is unavoidable, hence efficient methods to reduce quantum noise are highly demanded. While the general theory of quantum error correction offers a framework to achieve fully fault-tolerant computations [1-4], its required qubit overhead remains prohibitively high for today's quantum devices [5]. As an alternative to the currently challenging quantum error correction, quantum error mitigation (QEM) approaches were proposed [6-14]. Despite their fundamental limitations [15-17], they remain the main available tools for the current noisy intermediate-scale quantum (NISQ) devices [18] and for the upcoming early phases of fault-tolerant computing [19-21].

In this work, we empirically investigate how additional information provided by physics can improve the performance of QEM. The cornerstone idea of our approach is that the time-evolved state of any noiseless quantum system, at any time during the computation process, obeys a corresponding Schrödinger equation, so physical laws can be used to verify quantum computations. Unfortunately, this idea alone is of small practical use since a full quantum state tomography of exponentially many measurements is required for the perfect knowledge of the state. To employ this concept in practice, one has to dramatically reduce the number of necessary measurements, for instance by symmetry verification [22] or by using N-representability conditions [23].

We utilize the fact that the full dynamics of the idealized system of $N_Q$ qubits can be obtained from a quantum Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierar-

chy [24-27] of $4^{N_Q}$ equations. In classical computations, to avoid implementing an exponentially large system of coupled dynamical equations, one generally truncates the hierarchy by modeling the high-order correlators or by assuming their vanishing [28-30]. However, for strongly correlated systems or any general computational task, there are no (known) naturally small parameters justifying these truncations [30,31]. In quantum computations, truncations are no longer necessary, as one can directly measure any observable from the quantum device. One can then use corresponding equations from the BBGKY hierarchy to test the correctness of the measurements, and hence of the quantum computations.

In this paper, we use the important fact that the amount of terms in all hierarchical equations is bounded by $poly(N_Q)$, and by focusing on $poly(N_Q)$-large subsets of the hierarchy, only a polynomial in $N_Q$ amount of additional classical resources is needed for the above-mentioned tests. We employ this supplementary information from the BBGKY hierarchy to improve the zero-noise extrapolation (ZNE) method [9,32] in digital quantum simulations. Specifically, we formulate a novel QEM scheme and apply it to the Schwinger model [33] brought to quantum lattice simulations as a $\frac{1}{2}$-spin model in the particular implementation of [34]. The Schwinger model has been widely studied and used as a benchmark toy model for quantum computations, for instance in the recent works given in Refs. [34-41].

The paper is structured as follows. In Sec. II, we derive and describe the BBGKY hierarchy. Section III is dedicated to our QEM method: in Sec. III A we briefly review the ZNE scheme, in Sec. III B we describe how we select the BBGKY equations from the hierarchy, and in Sec. III C we present our QEM technique. In Sec. IV, we apply our method to the lattice Schwinger model. Finally, in Sec. V, we conclude and discuss further potential applications of the method.

<footnote>*Contact author: theo.saporiti@cea.fr</footnote>
## II. THE BBGKY HIERARCHY

We begin by deriving the BBGKY hierarchy and by discussing its properties.

Consider a quantum 1/2-spin model composed of $N_Q$ qubits, each of which is labeled by an index $i \in \{1, \dots, N_Q\} =: S$. Let $A \subseteq S$ represent a subsystem of the spin model, and let

$$ \sigma[A, (\mu_i)_{i \in A}] := \prod_{i \in A} \sigma_i^{\mu_i} \quad (1) $$

define a Pauli string, where $\mu_i \in \{1, 2, 3\}$ and $\sigma_i^{\mu_i}$ is the Pauli operator acting on the $i$th qubit in the $\mu_i$th direction. Assume the model has a Hamiltonian of the form

$$ H := \frac{1}{2} \sum_{i \in S} h_i^{\mu_i} \sigma_i^{\mu_i} + \frac{1}{4} \sum_{\substack{i,j \in S \\ i < j}} V_{ij}^{\mu_i \nu_i} \sigma_i^{\mu_i} \sigma_j^{\nu_i}, \quad (2) $$

where $h_i^{\mu_i}$ is the interaction term of the $i$th spin in the $\mu_i$th direction with an external magnetic field, $V_{ij}^{\mu_i \nu_i}$ is the interaction potential term among the $i$th and $j$th spins of the respective $\mu_i$th and $\nu_i$th directions, and where, from now on, Einstein's summation is implied. Moreover, here and throughout the work, we set $\hbar = c = 1$.

If one injects (1) and (2) into Ehrenfest's theorem and computes all the commutators, one obtains [42]

$$ \begin{aligned} \frac{d}{dt} \left\langle \prod_{i \in A} \sigma_i^{\mu_i} \right\rangle = & \sum_{\substack{i,j \in A \\ i \neq j}} \frac{V_{ij}^{\mu_i \nu_i}}{2} \epsilon_{\mu_j \nu_i \lambda} \left\langle \sigma_j^{\lambda} \prod_{k \in A \setminus \{i,j\}} \sigma_k^{\mu_k} \right\rangle \\ & + \sum_{i \in A} \hbar_i^{\lambda} \epsilon_{\mu_i \lambda \nu} \left\langle \sigma_i^{\nu} \prod_{j \in A \setminus \{i\}} \sigma_j^{\mu_j} \right\rangle \\ & + \sum_{\substack{i \in A \\ j \notin A}} \frac{V_{ij}^{\mu_i \nu_i}}{2} \epsilon_{\mu_i \mu \lambda} \left\langle \sigma_i^{\lambda} \sigma_j^{\nu} \prod_{k \in A \setminus \{i\}} \sigma_k^{\mu_k} \right\rangle, \end{aligned} \quad (3) $$

where $\epsilon_{\mu_i \nu_i \lambda}$ is the three-dimensional Levi-Civita symbol. We call (3) the BBGKY equation of the $\sigma[A, (\mu_i)_{i \in A}]$ Pauli string. This is because if one considers all the Pauli strings of all possible directions [namely, all partitions $A \subseteq S$ of all possible $(\mu_i)_{i \in A}$], then by computing all their associated BBGKY equations (3), the complete exponentially large BBGKY-like hierarchy is generated. More precisely, for a specific Pauli string of length $|A| = n \le N_Q$, the time derivative of that $n$-point correlator is determined by a linear combination of $(n-1)$-point, $n$-point, and $(n+1)$-point correlators, respectively, found in the first, second, and third summations of (3), all selected according to the values $h_i^{\mu_i}$ and $V_{ij}^{\mu_i \nu_i}$ in (2). The right-hand side (rhs) of (3) contains up to $9n(n-1)$-many $(n-1)$-point correlators, up to $9n$-many $n$-point correlators, and up to $27n(N_Q-n)$-many $(n+1)$-point correlators. Importantly, this implies that the amount of correlators in the rhs of (3) is polynomial in $n$ and $N_Q$, and bounded by $81N_Q^2/8$.

As a last remark, note that (2) is constructed out of one-site and two-site interaction terms. A more general Hamiltonian containing higher-order Pauli strings would generate different BBGKY equations, with more terms in their rhs. However, this would not spoil the applicability of our method, as the resulting (3) equations would remain linear combinations of

expectation values. In Sec. III C, we explain why this property is fundamental for the formulation of our mitigation scheme.

## III. MITIGATION TECHNIQUE

In this section, we provide a short review of ZNE in the context of time evolution. Then we present our method, a BBGKY-improved ZNE scheme, which we formulate as a postprocessing linear least-squares (LLS) optimization procedure. Our method incorporates a poly($N_Q$)-large subset of the BBGKY hierarchy evaluated across all time points. We consider the list of $l$ Pauli strings $\{Q_q\}_{q \in \{0, \dots, l-1\}}$, where for brevity $Q_q$ is defined as in (1). Our goal is to mitigate the corresponding measurements obtained from a realistic quantum device. We focus solely on Pauli strings because they form an operator basis for the observables of the system.

### A. The ZNE method

By Trotterization, the evolution time $T$ is discretized into $N$ slices of duration $\Delta t := T/N$, and evolution steps are obtained thanks to a Suzuki-Trotter decomposition scheme of the order of $d_{ST}$ [43–46]. Then, different realizations of the quantum circuit implementing the time evolution are generated, each of them containing local unitary foldings [10] with a frequency of $\eta \ge 0$ insertions per step. Under the assumption that unitary foldings are affected by the same kind of noise as regular evolution steps, this implies an error level at the $s$th step $s \in \{0, \dots, N\}$ relative to the original $\eta = 0$ circuit of

$$ 1 \le \epsilon_{s\eta} := \frac{s + 2\lfloor \eta s \rfloor}{s} \xrightarrow{s \to \infty} 2\eta + 1. \quad (4) $$

Performing the above $m$ times with different noise levels $\eta \in (\eta_1, \dots, \eta_m) := \bar{\eta}$, at each time point $t_s := s\Delta t$, we end up with an experimentally measured set of data points, $(\epsilon_{s\eta}, \langle Q_q \rangle_{s\eta})$, where $\langle Q_q \rangle_{s\eta}$ is the measurement of $Q_q$ at time point $t_s$, estimated with $N_S$ shots, under the $\eta$ noise level [47]. For a given quantity at a given time point, these experimentally measured points can be fitted across the error levels with a chosen model function of $\epsilon$, leading to a zero-noise $\epsilon \to 0$ extrapolated $\langle Q_q^0 \rangle_s$ [9,10]. Common choices of fitting strategies include polynomials of degree $d \le m-1$, (poly-)exponential Ansätze, and Richardson's extrapolation [10].

### B. Selection of BBGKY equations

Given the set of expectation values, $\{\langle Q_q \rangle_{q \in \{0, \dots, l-1\}}\}$, the physical knowledge provided by a BBGKY equation can help their mitigation if and only if its associated Pauli string is hierarchically connected to any of the $\{Q_q\}_{q \in \{0, \dots, l-1\}}$.

For a given generic $\sigma[A, (\mu_i)_{i \in A}]$, the rhs of Eq. (3) provides all of the correlators $\sigma[B, (v_i)_{i \in B}]$, with $B \subseteq S$ of directions $(v_i)_{i \in B}$, that are connected to its time evolution via the hierarchy. In that case, we say that the $\sigma[B, (v_i)_{i \in B}]$ are downstream connected to $\sigma[A, (\mu_i)_{i \in A}]$. We now want to determine the inverse, that is, given a Pauli string $\sigma[B, (v_i)_{i \in B}]$, find all correlators $\sigma[A, (\mu_i)_{i \in A}]$ generating $\sigma[B, (v_i)_{i \in B}]$ in their rhs of (3). In that case, we say that $\sigma[A, (\mu_i)_{i \in A}]$ is upstream connected to $\sigma[B, (v_i)_{i \in B}]$. Overall, we say that a correlator is connected to a Pauli string of interest if it is either downstream or upstream connected to that Pauli string.
<figure><img src="image_3.png" /><figcaption>FIG. 1. In this diagram, *b* is a Pauli string of interest, and *a*₁, *c*₁ are (n + 1)-point correlators, *a*₂, *b*, *c*₂ are n-point correlators, and *a*₃, *c*₃ are (n − 1)-point correlators. All *c*₁, *c*₂, *c*₃ are downstream connected to *b*, and all *a*₁, *a*₂, *a*₃ are upstream connected to *b*.</figcaption></figure>

To find all upstream connected correlators $\sigma[A, (\mu_i)_{i \in A}]$ to a given $\sigma[B, (\nu_i)_{i \in B}]$, there are only three possibilities: denoting $n_A = |A|$ and $n_B = |B|$, $\sigma[B, (\nu_i)_{i \in B}]$ can appear in the BBGKY equation of $\sigma[A, (\mu_i)_{i \in A}]$ either as an $(n_A + 1)$-point, $n_A$-point, or $(n_A - 1)$-point correlator. The selection rules to find $\sigma[A, (\mu_i)_{i \in A}]$ candidates in all three cases are derived in Appendix A. There we also show that, importantly, the upper bound for the amount of Ansätze one has to check to pick up all the upstream connected $\sigma[A, (\mu_i)_{i \in A}]$ is polynomial in $N_Q$. Figure 1 schematically summarizes the two possible kinds of connections.

For our mitigation purposes, we select the BBGKY equa-
tions associated to the subset of the BBGKY hierarchy whose
Pauli strings are connected to any of the {$Q_q$}$_{q \in \{0,...,l-1\}}$ by,
at most, r connections. The subset is obtained by iteratively
computing over itself all of its downstream and upstream con-
nected Pauli strings, as explained in Appendix A, for a total of
r iterations. This generates a list of g BBGKY equations con-
taining $\Lambda \ge l$ correlators, the actual amount of expectation
values {$\langle Q_q \rangle$}$_{q \in \{0,...,l-1,l,...,\Lambda-1\}}$ that have to be measured.

C. Our BBGKY-improved ZNE method

Our method aims at improving the ZNE-obtained measure-
ments ⟨Q<sup>0</sup><sub>q</sub>⟩<sub>s</sub> with better BBGKY-constrained ⟨Q<sup>θ</sup><sub>q</sub>⟩<sub>s</sub> extrapola-
tions. As a first application of the above idea, in this work we
use a least-squares polynomial (LSP) fitting strategy for ZNE.
This is motivated by two reasons: generality of the noise and
convex optimization. Regarding the first argument, we do not
make any assumptions about the nature of the quantum noise,
so that our method is independent of the particular employed
quantum device. By the Weierstrass approximation theorem
[48], any underlying noise dependence of the measurements
can be arbitrarily approached by a polynomial, up to an ac-
cepted error. Therefore, the choice of a LSP is sensible in
the absence of any assumptions about the noise regime. It
is then up to the user to select the degree d of the LSP and
accept its associated Weierstrass approximation error. For the
second argument, the use of LSPs ensures that ZNE can be
formulated as a LLS problem, rendering the fitting (and hence
the mitigation) a convex problem, for which a unique solution
exists and is guaranteed to be found.

In this setting, where ZNE is a LLS minimization
procedure, we want to incorporate the selected BBGKY
equations as linear combinations of ⟨Q<sup>θ</sup><sub>q</sub>⟩<sub>s</sub>. To do so, we

approximate the time derivatives in (3) with derivatives of a
Bernstein polynomial, fitting the N + 1 extrapolated expec-
tation values. We use (derivatives of) Bernstein polynomials
because they are linear in the ⟨Q<sup>θ</sup><sub>q</sub>⟩<sub>s</sub> and because they uni-
formly converge to the functions (derivatives) they are fitting,
with an error of the order O(1/N) [49]. Starting from the
Bernstein polynomial [48]

$$
\langle Q_q^\theta \rangle (t) := \sum_{s=0}^{N} \langle Q_q^\theta \rangle_b b_{sN} \left( \frac{t}{T} \right), \quad (5)
$$

constructed out of the $N + 1$ Bernstein polynomial basis el-
ements $b_{sN}(x)$ of degree $N$, with $x \in [0, 1]$ hence $t \in [0, T]$,
then

$$
\frac{d}{dt} \langle Q_q^\theta \rangle (t) = \sum_{s=0}^{N} \langle Q_q^\theta \rangle_b \beta_{sN} \left( \frac{t}{T} \right), \quad (6)
$$

where we define

$$
\beta_{sN}(x) := \frac{1}{\Delta t} \begin{cases} -(1-x)^{N-1} & \text{if } s=0 \\ x^{N-1} & \text{if } s=N \\ b_{s-1,N-1}(x) - b_{s,N-1}(x) & \text{otherwise.} \end{cases} \quad (7)
$$

Using the approximate time derivatives (6), approxima-
tions of the BBGKY equations (3) can be expressed as linear
combinations of ⟨Q<sup>θ</sup><sub>q</sub>⟩<sub>s</sub>. Thereby, BBGKY equations can be
cast inside the original ZNE LLS procedure as additional
constraints of the minimization problem,

$$
\operatorname{argmin}_{\bar{a}} \frac{1}{2} \|\vec{v}(\bar{a})\|^2 \quad \text{with} \quad \vec{v}(\bar{a}) := M\bar{a} - \bar{y}, \qquad (8)
$$

where $\bar{a}$ will be defined later in (11) and where

$$
M := \left( \frac{M_{\bar{\eta}}}{G} \right) \quad \text{with} \quad M_{\bar{\eta}} := \bigoplus_{q=0}^{\Lambda-1} \bigoplus_{s=1}^{N} \underbrace{\begin{pmatrix} \epsilon_{s\eta_1}^d & \cdots & 1 \\ \vdots & \ddots & \vdots \\ \epsilon_{s\eta_m}^d & \cdots & 1 \end{pmatrix}}_{=:M_{s\bar{\eta}}} . \tag{9}
$$

Here, $M_{s\bar{\eta}}$ is a Vandermonde-like matrix, $G$ is filled with $g(N+1)$ lines of appropriate entries encoding the corresponding BBGKY equations at every time point $t_s$, and the target vector $\bar{y} := [(\bar{y}_{01\bar{\eta}}, \dots, \bar{y}_{0N\bar{\eta}}), \dots, (\bar{y}_{\Lambda-1,1\bar{\eta}}, \dots, \bar{y}_{\Lambda-1,N\bar{\eta}}), \bar{g}]$ sequentially groups all $s \ge 1$ experimental measurements $\bar{y}_{qs\bar{\eta}} := (\langle Q_q \rangle_{s\eta_1}, \dots, \langle Q_q \rangle_{s\eta_m})$ together with the $g(N+1)$-long $\bar{g}$ vector encoding the $s=0$ expectation values of the BBGKY equations. We do not mitigate any $\langle Q_q^\theta \rangle_0$ because they can be numerically computed with arbitrary precision at $t=0$, and hence they are known *a priori*. Notice that by setting $g=0$, one recovers the original ZNE procedures, in which case $\langle Q_q^\theta \rangle = \langle Q_q^0 \rangle$ and $M$ decouples into $N$ independent ZNEs, each minimizing the error on

$$
\langle Q_q \rangle_{s\eta} \approx \langle Q_q^0 \rangle_s + \sum_{\delta=1}^{d} a_{qs\delta} \varepsilon_{s\eta}^{\delta}, \quad (10)
$$

where all LSP coefficients $a_{qs\delta}$ are packed into

$$
\bar{a} := [(\bar{a}_{01}, \dots, \bar{a}_{0N}), \dots, (\bar{a}_{\Lambda-1,1}, \dots, \bar{a}_{\Lambda-1,N})], \quad (11)
$$

with $\bar{a}_{qs} := (a_{qsd}, \dots, a_{qs1}, \langle Q_q^\theta \rangle_s)$. In particular, $\langle Q_q^\theta \rangle_s$ can be extracted from $(\bar{a})_p$ at index $p = p(q, s) = (d+1) + (s-$
<figure><img src="image_4.png" /><figcaption>FIG. 2. Schematic depiction of our method: each slice of data (red crosses) is fitted with a LSP (dotted gray lines), producing a series of intermediate ZNEs (dark red points), which in turn are fitted with a Bernstein polynomial (purple line), giving access to time derivatives at each time point (purple squares). These are then used to improve the ZNEs thanks to the BBGKY equations, producing different LSPs (continuous gray lines) leading to BBGKY-improved ZNEs (black circles). The horizontal axis represents time, the vertical axes the expectation values, and the diagonal axes the error levels, where the dashed black lines represent the $\epsilon \ge 1$ bound.</figcaption></figure>

$1)(d+1)+q(d+1)N$. Finally, notice that $M$ is a rectangular matrix of polynomial size $[mN\Lambda + g(N+1)] \times (d+1)N\Lambda$ and that $G$ couples together all extrapolated $\langle Q_q^\theta \rangle_s$ in two ways: across all time points, as in (6), and according to their connections in the hierarchy, as in (3). Figure 2 graphically represents our method, and in Appendix C we give an example of an $M$ matrix.

## IV. RESULTS OF THE MITIGATION

We now briefly overview the lattice Schwinger model and the quantities we want to mitigate with our method. We then show and discuss the obtained numerical results.

### A. The lattice Schwinger model

The Schwinger model [33] describes one-dimensional quantum electrodynamics. This continuous model can be brought to its lattice Hamiltonian formulation via Kogut-Susskind construction [50]. Then, the original degrees of freedom can be recast into quantum $\frac{1}{2}$-spins [51] with open boundary conditions [34]. We test our method on the latter, whose (dimensionless) Hamiltonian is

$$ H := -\frac{m}{g} \sqrt{x} \sum_{i \in S} (-1)^i \sigma_i^3 + \sum_{i=1}^{N_Q-1} \left( \frac{N_Q}{4} - \frac{1}{2} \left[ \frac{i-1}{2} \right] + l_0(N_Q - i) \right) \sigma_i^3 \\ + \frac{x}{2} \sum_{i=1}^{N_Q-1} (\sigma_i^1 \sigma_{i+1}^1 + \sigma_i^2 \sigma_{i+1}^2) + \frac{1}{2} \sum_{\substack{i,j \in S \\ i<j}} (N_Q - j + \lambda) \sigma_i^3 \sigma_j^3, \quad (12) $$

where $m/g$ is the lattice mass over coupling ratio, $x = (N_Q/V)^2$ with $V$ the (dimensionless) lattice volume, $l_0$ the background electric field, $\lambda \gg 1$ a Lagrange multiplier to restrict simulations within the vanishing total charge sector, and a final constant term was disregarded. The Hamiltonian (12) is in the appropriate form of (2), and we are interested in the Pauli strings of the quantities

$$ Q := \frac{1}{2} \sum_{i \in S} \sigma_i^3 \quad \text{and} \quad P := \frac{N_Q}{2} - \frac{1}{2} \sum_{i \in S} (-1)^i \sigma_i^3. \quad (13) $$

These are, respectively, the electric charge operator and the particle number operator. Moreover, $[Q, H] = 0$ and $[P, H] \ne 0$, meaning that they represent two distinct behaviors over which we can test our method: $\langle P \rangle$ will vary in time, while $\langle Q \rangle$ will stay constant. In the following, we will often employ the abuses of notation $Q_q = Q, P$ with $q = Q, P$ to indicate the mitigation of the above linear combinations of Pauli strings.

### B. Numerical framework

We assess the effectiveness of our method against ZNE with, respectively, the following 2-norms:

$$ L_q^\theta := \sqrt{\Delta t \sum_{s=0}^{N} [\langle Q_q^\theta \rangle_s - \langle Q_q \rangle(t_s)]^2}, \quad (14) $$

$$ L_q^0 := \sqrt{\Delta t \sum_{s=0}^{N} [\langle Q_q^0 \rangle_s - \langle Q_q \rangle(t_s)]^2}, \quad (15) $$

which quantify the accumulated error of the extrapolations against the exact-diagonalization (ED) evolution $\langle Q_q \rangle(t)$ over all time points $t_s$. All of our computations are performed in QISKIT 1.3 [52] within a simulated quantum device whose realistic noise model is generated in real time from the back-end physical properties of the IBM Brisbane quantum processor [53]. In the following, we fix the parameters $N_Q = 4$, initial state $|0101\rangle$, $N_S = 10240$, $N = 20$, $T = 4$, $\bar{\eta} = (0, 1, 1.5, 2)$, $d_{ST} = 1$, $d = 2$, $\lambda = 100$, and $V = 30$ [34]. Unless otherwise stated, we fix $r = 0$, while the remaining parameters $m/g$ and $l_0$ will vary throughout these simulations. We systematically check through idealized noiseless estimations ($N_S \to \infty$) of simulations that, for measurements of $\langle Q \rangle$ and $\langle P \rangle$, their total Trotter errors of the order of $O[N(e\Delta t)^{d_{ST}+1}]$, where $e$ is one unit of energy, are, respectively, smaller or equal to $1.9 \times 10^{-16}$ and 0.02 with respect to ED evolutions. Therefore, as the Trotter error is negligible, only the shot noise of the order of $O(1/\sqrt{N_S})$ is considered for the computation of the errors associated to $L_q^\theta$ and $L_q^0$. These are, respectively, denoted $\Delta L_q^\theta$ and $\Delta L_q^0$ or, respectively, presented as $L_q^0 \pm \Delta L_q^0$ and $L_q^0 \pm \Delta L_q^0$, and they estimate the standard deviation of the
<figure><img src="image_5.png" /><figcaption>FIG. 3. Top panel: Improvement of the ZNE mitigation for the particle number $\langle P \rangle$ with respect to the ED evolution. Bottom panel: Advantage of our method compared to the top panel. In both panels, $\Delta L_p^0 \in [0.00891, 0.00920]$ and $\Delta L_p^\theta \in [0.00562, 0.00774]$.</figcaption></figure>

two kinds of mitigations. Throughout the work, all standard deviation estimations are computed with the conventional propagation of errors formula, assuming independent variables and truncating at the leading order.

### C. Numerical results

Figure 3 shows, for the $P$ observable, a parameter scan of $10 \times 10$ blocks over the region $(l_0, m/g) \in [0, 1.5]^2$ where the error of ZNE with respect to the ED dynamics (top panel) and the improvement of our method with respect to ZNE (bottom panel) are computed. The top panel gives us the scale of the ZNE error, and the bottom panel tells us by how much that error was reduced within our BBGKY-improved scheme. A systematic improvement over the entire parameter region is manifest by the presence of only negative $L_p^\theta - L_p^0$ values and, taking average values as in Table I, it is approximately $(18.2 \pm 0.5)\%$. In all simulations $\Delta L_p^\theta < \Delta L_p^0$, meaning that our method reduces the estimated standard deviation of ZNE.

<table><caption>TABLE I. First column: Average errors of ZNE against ED. Second column: Average errors of our method against ED. Third column: Absolute average improvements of our method against ZNE. Fourth column: Relative average improvements of our method against ZNE. The averages for the $q = P, Q$ observables are computed, respectively, from the data of Figs. 3 and 5.</caption><thead><tr><th>q</th><th>L<sup>0</sup><sub>q</sub> (ZNE)</th><th>L<sup>θ</sup><sub>q</sub> (Ours)</th><th>L<sup>0</sup><sub>q</sub> - L<sup>θ</sup><sub>q</sub></th><th>1 - L<sup>θ</sup><sub>q</sub>/L<sup>0</sup><sub>q</sub></th></tr></thead><tbody><tr><td>P</td><td>1.975 ± 0.009</td><td>1.659 ± 0.007</td><td>0.316 ± 0.011</td><td>(18.2 ± 0.5)%</td></tr><tr><td>Q</td><td>0.128 ± 0.008</td><td>0.060 ± 0.007</td><td>0.068 ± 0.011</td><td>(52.8 ± 6.3)%</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIG. 4. Evolution of the particle number $\langle P \rangle$ in the $(l_0, m/g) = (0, 0.15)$ regime. The red crosses are measurements at increasing noise levels (shot noise not shown, of average value 0.01), the ZNEs are represented in dark red (estimated standard deviations shown as a band), the BBGKY-improved ZNEs are represented in black (estimated standard deviations shown as a band), the purple line is the Bernstein polynomial associated to the latter, and the dashed gray line is the ED evolution.</figcaption></figure>

The diagonal line in both panels can be explained by a system-dependent artifact caused by the alignment of the ED dynamics to the saturation, or flattening, of the measurements as in Fig. 4.

Figure 4 shows the time evolution of the $(l_0, m/g) = (0, 0.15)$ block of Fig. 3. We see that after ~15 Trotter steps, the measurements saturate and recovering the original dynamics becomes challenging. Nevertheless, thanks to the additional BBGKY constraints, we see that the Bernstein polynomial correctly tries to match the time derivatives of the ED evolution.

Figure 5 shows the contents of Fig. 3, but for the $Q$ observable. Again, a systematic improvement over the entire parameter region is observed with our method and, from Table I, we see that it is approximately $(52.8 \pm 6.3)\%$ with respect to ZNE. Once more the reduction $\Delta L_Q^\theta < \Delta L_Q^0$ of estimated standard deviation by our method is confirmed on all simulations. Here, the bottom-left artifact region of the bottom panel can be explained by the small values of $(l_0, m/g)$, reducing the importance of the BBGKY equations in the LLS minimization.

Figure 6 shows the time evolution of the $(l_0, m/g) = (0, 0)$ block of Fig. 5. Here, no saturation phenomenon occurs because all measurements remain equally noisy. Nevertheless, thanks to the additional BBGKY constraints, we see again that the Bernstein polynomial correctly tries to match the null time derivative of the conserved quantity $Q$.

Table I summarizes as averages the errors, the absolute and the relative improvements of the previous two parameter scans, and their associated standard deviation estimations, shown in Figs. 3 and 5. Again, with our method, we see a systematic improvement of the two mitigations with respect to ZNE, for both the nonconserved $P$ and the conserved $Q$, although we observe a larger relative improvement in the mitigation of $\langle Q \rangle$. This is because noise concentrates around the ED evolution in Fig. 6 so, in the minimization of the
<figure><img src="image_7.png" /><figcaption>FIG. 5. Top panel: Improvement of the ZNE mitigation for the electric charge $\langle Q \rangle$ with respect to the ED evolution. Bottom panel: Advantage of our method compared to the top panel. In both panels, $\Delta L_Q^0 \in [0.00831, 0.00836]$ and $\Delta L_Q^\theta \in [0.00706, 0.00717]$.</figcaption></figure>

LLS problem, the Bernstein polynomial is less penalized in deviating from the ZNE to match the null time derivative of $\langle Q \rangle$.

Figure 7 shows the error committed by ZNE and our BBGKY-informed method over the ED evolution of $\langle Q \rangle$ in the ($l_0, m/g$) = (0, 0) regime, using measurements estimated with an increasing number $N_S$ of shots. The top panel displays a comparable convergence of both mitigations to corresponding limiting errors as $N_S$ increases. Moreover, the bottom panel highlights the similarities in the magnitudes of the corresponding estimated standard deviations. Therefore, as our BBGKY-informed scheme behaves similarly to ZNE against

<figure><img src="image_8.png" /><figcaption>FIG. 6. Evolution of the electric charge $\langle Q \rangle$ in the ($l_0, m/g$) = (0, 0) regime. The same symbols and colors of Fig. 4 are employed, and the same average shot noise of the raw measurements is observed.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7. Top panel: Errors for the mitigation of the charge $\langle Q \rangle$, in the ($l_0, m/g$) = (0, 0) regime at different $N_S$, with ZNE (data represented in dark red) and our method (data represented in black). Bottom panel: Magnitudes of the above panel's error bars.</figcaption></figure>

statistical noise, the inherited stability and robustness of the LLS ZNE scheme are recovered in our method as well. Comparable results are obtained for the mitigation of $\langle P \rangle$ in the ($l_0, m/g$) = (0, 0) block.

We now study how the size of the selected subset of BBGKY equations affects the results of the mitigation. This is displayed in Fig. 8, where the mitigation of $\langle P \rangle$ in the ($l_0, m/g$) = (0, 0) block is repeated for different maximal connections radii $r$. The top panel displays the errors of the ZNE and BBGKY-improved mitigations, while the bottom panel quantifies how many BBGKY equations $g$ cover the dynamics of the measured $\Lambda$ quantities.

A ratio of $g/\Lambda = 1$, which is reached at $r = 6$ and remains so at $r = 7$, means that the subset of equations is fully

<figure><img src="image_10.png" /><figcaption>FIG. 8. Top panel: Errors for the mitigation of the particle number $\langle P \rangle$, in the ($l_0, m/g$) = (0, 0) regime at different $r$, with ZNE (data represented in dark red) and our method (data represented in black). Bottom panel: Determination of the hierarchical subset.</figcaption></figure>
determined with respect to its unknowns. In the present case, $\Lambda = 126 < 256 = 4^{N_0}$ because the BBGKY hierarchy splits into four independent hierarchies of sizes 1, 1, 126, and 128. In particular, the two first hierarchies are composed, respectively, of solely the identity and $\sigma_1^3\sigma_2^3\sigma_3^3\sigma_4^3$, which are conserved quantities of (12). The third hierarchy involves $\{\sigma_i^3\}_{i \in S}$, which are required to build the Q and P observables.

We see that the inclusion of additional BBGKY constraints reduces the error with respect to ZNE immediately starting from $r=0$, in accordance with Fig. 3. In our experiments with the Schwinger model, no clear pattern regarding the reduction or the increase of absolute errors can be stated for $r > 0$, either in this ($l_0, m/g) = (0, 0)$ simulation or in other different points of the parameter scan. We conjecture that a similar behavior is valid for other systems and therefore the optimal regime of our method lies close to a small $r=0$ hierarchical subset. This is because, given the ratio of the number of lines in the $M$ matrix distributed among ZNE and the BBGKY equations, the LLS problem prevents the Bernstein polynomial from deviating too much from the original ZNEs even as $g$ increases. Moreover, the number of time points, $N$, may be too small for the Bernstein polynomial to benefit from its uniform convergence, leading to a systematic non-negligible additional approximation error from (6). Comparable results and conclusions are observed for the mitigation of $\langle Q \rangle$ in the ($l_0, m/g) = (0, 0)$ block.

## V. CONCLUSIONS

In this paper, we considered executions of quantum algorithms as time evolutions of idealized systems. Their noiseless dynamics are governed by physical laws, which can be employed to mitigate the quantum errors arising from their physical realizations. For this purpose, we derived a corresponding BBGKY-like hierarchy and selected a poly($N_Q$)-large subset of its equations. We proposed a novel QEM scheme encapsulating these supplementary physics-informed constraints into the ZNE procedure. A LLS problem is thereby obtained, whose required classical computational resources scale polynomially in $N_Q$. We numerically investigated the effectiveness of our method on digital quantum simulations, mimicking realistic quantum hardware noise, of the lattice Schwinger model.

By applying our method to the lattice Schwinger model, we assessed our BBGKY-informed QEM scheme against ZNE by comparing the mitigation of quantum observables to known ED evolutions. It was found that in the considered regions of the parameter scan and under the selected input parameters, our method systematically improves the QEM of ZNE measurements. Moreover, with our method, it was found that the range of relative improvement of the error with respect to ZNE spans from $(18.2 \pm 0.5)\%$ to $(52.8 \pm 6.3)\%$, depending on whether or not the mitigated quantity is conserved. It was also found that our method exhibits the same stability and robustness properties of conventional ZNE. Finally, it was found that the maximal connections radius of the selected hierarchical subset should be close to $r=0$, as such a subset of BBGKY equations already guarantees a systematic improvement over ZNE.

Further expansion of this work includes considering Hamiltonians encoding more general multispin interactions, performing imaginary-time evolution, and evolution with time-dependent Hamiltonians. The latter would pave the way to the mitigation of adiabatic time evolution. This would also allow the mitigation of quantum circuits, if implemented as Trotterized time evolutions of a time-dependent system. A more detailed analysis of the stability and robustness of our method could be investigated and different nonlinear fitting strategies other than LSPs could be implemented and tested. Finally, the supplementary physical knowledge of the BBGKY hierarchy could be used to help not only ZNE, but other mitigation schemes as well, either on digital or analog quantum machines. This could be done either as an entirely postprocessing procedure or directly affecting the parameters of quantum computation, as in variational methods.

## ACKNOWLEDGMENTS

This work has received support from the French State managed by the National Research Agency under the France 2030 program with Reference No. ANR-22-PNCQ-0002. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.

## DATA AVAILABILITY

The data that support the findings of this article are openly available [54].

## APPENDIX A: UPSTREAM CONNECTED CORRELATORS

Here we derive an algorithm composed of three subroutines to obtain all $\sigma[A, (\mu_k)_{k \in A}]$ correlators upstream connected to the target $\sigma[B, (v_k)_{k \in B}]$. First, define the function $\bar{\epsilon}(\mu, v) := \lambda|\epsilon_{\mu v\lambda}|$, where Einstein's summation is still implied, to be interpreted as that index such that $\epsilon_{\mu v\bar{\epsilon}(\mu, v)} \neq 0$. Let, also, $n_A = |A|$ and $n_B = |B|$.

For the $n_A = n_B$ case, we observe in the second summation of (3) that for $\sigma[B, (v_k)_{k \in B}]$ to appear in the expectation value, it must be $A = B$. Regarding the directions $(\mu_k)_{k \in A}$, selecting an $i \in A$ ($n_B$ choices) and letting $\mu_k = v_k$ for all $k \in A \setminus \{i\}$, we see that if we want $v$ to pick up the desired $v = v_i$ value, then, because of the Levi-Civita symbol, it must be $\mu_i \neq v_i$ (two choices). Then, $\lambda = \bar{\epsilon}(\mu_i, v_i)$, and hence $\sigma[B, (v_i)_{i \in B}]$ appears in the second summation only if

$$h_i^{\bar{\epsilon}(\mu_i, v_i)} \neq 0. \qquad (A1)$$

In total, $2n_B$ possible $\sigma[A, (\mu_k)_{k \in A}]$ Ansätze have to be checked against the magnetic field.

For the $n_A = n_B - 1$ case, we observe in the third summation of (3) that for $\sigma[B, (v_k)_{k \in B}]$ to appear in the expectation value, it must be $A = B \setminus \{j\}$ for a selected $j \in B$ ($n_B$ choices). Regarding the directions $(\mu_k)_{k \in A}$, selecting a specific $i \in A$ ($n_B - 1$ choices) and letting $\mu_k = v_k$ for all $k \in A \setminus \{i\}$, we see that if we want $\lambda$ to pick up the desired $\lambda = v_i$ value, then, again because of the Levi-Civita symbol, it must be $\mu_i \neq v_i$ (two choices). Then, $\lambda = \bar{\epsilon}(\mu_i, v_i)$, and hence $\sigma[B, (v_k)_{k \in B}]$
appears in the third summation only if

$$V_{ij}^{\bar{\epsilon}(\mu_i, \nu_i)v_j} \neq 0. \qquad (A2)$$

Notice that the index $v$ will necessarily pick $v = v_j$ because it is contracted independently. In total, $2n_B(n_B - 1)$ possible $\sigma[A, (\mu_k)_{k \in A}]$ Ansätze have to be checked against the interaction potential.

For the $n_A = n_B + 1$ case, we observe in the first summation of (3) that for $\sigma[B, (v_k)_{k \in B}]$ to appear in the expectation value, it must be $A = B \cup \{i\}$ for a selected $i \in S \setminus B$ ($N_Q - n_B$ choices). Regarding the directions $(\mu_k)_{k \in A}$, selecting a specific $j \in B$ ($n_B$ choices) and letting $\mu_k = v_k$ for all $k \in A \setminus \{i, j\}$ and $\mu_i \in \{1, 2, 3\}$ (three choices), we see that if we want $\lambda$ to pick up the desired $\lambda = v_j$ value, then, once again because of the Levi-Civita symbol, it must be $\mu_j \neq v_j$ (two choices). Then, $v = \bar{\epsilon}(\mu_j, v_j)$, and hence $\sigma[B, (v_k)_{k \in B}]$ appears in the first summation only if

$$V_{ij}^{\mu_i \bar{\epsilon}(\mu_j, v_j)} \neq 0. \qquad (A3)$$

In total, $6n_B(N_Q - n_B)$ possible $\sigma[A, (\mu_k)_{k \in A}]$ Ansätze have to be checked against the interaction potential.

To summarize, there is a polynomial in $n_B$ and $N_Q$ amount of Ansätze to check against the coefficients of (2), and the number of checks is overall bounded by $9N_Q^2/4$.

## APPENDIX B: RANDOM SHIFT OF ERROR LEVELS

Every time (4) is employed to measure a $(\epsilon_{s\eta}, \langle Q_q \rangle_{s\eta})$ data point, a small random shift $\epsilon_{s\eta} \rightarrow \epsilon_{s\eta} + \chi$ is performed, where $\chi$ is a realization of the normally distributed random variable $X \sim N(0, 1/N_S)$. This is to avoid cases where two different noise levels $\eta_1 \neq \eta_2$ produce the same error level, for example, $\eta_1 = 1$ and $\eta_2 = 1.5$ leading to $\epsilon_{1\eta_1} = \epsilon_{1\eta_2} = 3$, thereby introducing artifacts in the LSP fitting. The variance of $X$ is chosen to be $1/N_S$ to mimic the standard deviation estimator when $N_S \rightarrow \infty$ and also because, in that limit, the $\chi$ shift should not affect $\epsilon_{s\eta}$, as quantum errors had an infinite $N_S \rightarrow \infty$ amount of possibilities to arise.

## APPENDIX C: EXAMPLE OF AN M MATRIX

To illustrate the construction of $M$ with a simple yet nontrivial example, consider the mitigation of $l = \Lambda = 2$ quantities after $N = 2$ Trotter steps with the $g = 2$ fictitious BBGKY equations,

$$\frac{d}{dt}\langle Q_0 \rangle = 0 \quad \text{and} \quad \frac{d}{dt}\langle Q_1 \rangle = V\langle Q_0 \rangle, \qquad (C1)$$

where $V \neq 0$ is constant and $Q_0, Q_1$ are, respectively, n-point and $(n+1)$-point correlators. Then, in this example, the $\vec{v}$ vector is given by

<figure><img src="image_11.png" /></figure>

(C2)
[1] P. W. Shor, Scheme for reducing decoherence in quantum computer memory, Phys. Rev. A **52**, R2493(R) (1995).

[2] A. M. Steane, Error correcting codes in quantum theory, Phys. Rev. Lett. **77**, 793 (1996).

[3] A. R. Calderbank and P. W. Shor, Good quantum error-correcting codes exist, Phys. Rev. A **54**, 1098 (1996).

[4] M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information: 10th Anniversary Edition* (Cambridge University Press, Cambridge, 2010).

[5] A. Chatterjee, A. Ghosh, and S. Ghosh, Quantum prometheus: Defying overhead with recycled ancillas in quantum error correction, arXiv:2411.12813.

[6] N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, Dynamical decoupling for superconducting qubits: A performance survey, Phys. Rev. Appl. **20**, 064027 (2023).

[7] J. J. Wallman and J. Emerson, Noise tailoring for scalable quantum computation via randomized compiling, Phys. Rev. A **94**, 052325 (2016).

[8] E. van den Berg, Z. K. Minev, and K. Temme, Model-free readout-error mitigation for quantum expectation values, Phys. Rev. A **105**, 032620 (2022).

[9] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. **119**, 180509 (2017).

[10] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, Digital zero noise extrapolation for quantum error mitigation, in *2020 IEEE International Conference on Quantum Computing and Engineering (QCE)* (Denver, CO, USA, 2020), pp. 306–316.

[11] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Phys. Rev. X **8**, 031027 (2018).

[12] A. Strikis, D. Qin, Y. Chen, S. C. Benjamin, and Y. Li, Learning-based quantum error mitigation, PRX Quantum **2**, 040330 (2021).

[13] C. Kim, K. Park, and J. Rhee, Quantum error mitigation with artificial neural network, IEEE Access **8**, 18853 (2020).

[14] Z. Cai, X. Xu, and S. C. Benjamin, Mitigating coherent noise using pauli conjugation, npj Quantum Inf. **6**, 17 (2020).

[15] R. Takagi, S. Endo, S. Minagawa, and M. Gu, Fundamental limits of quantum error mitigation, npj Quantum Inf. **8**, 114 (2022).

[16] R. Takagi, H. Tajima, and M. Gu, Universal sampling lower bounds for quantum error mitigation, Phys. Rev. Lett. **131**, 210602 (2023).

[17] Y. Quek, D. Stilck França, S. Khatri, J. J. Meyer, and J. Eisert, Exponentially tighter bounds on limitations of quantum error mitigation, Nat. Phys. **20**, 1648 (2024).

[18] D. Pomarico, M. Pandey, R. Cioli, F. Dell'Anna, S. Pascazio, F. V. Pepe, P. Facchi, and E. Ercolessi, Quantum error mitigation in optimized circuits for particle-density correlations in real-time dynamics of the Schwinger model, Entropy **27**, 427 (2025).

[19] Y. Suzuki, S. Endo, K. Fujii, and Y. Tokunaga, Quantum error mitigation as a universal error reduction technique: Applications from the NISQ to the fault-tolerant quantum computing eras, PRX Quantum **3**, 010345 (2022).

[20] Z. Zimborás et al., Myths around quantum computation before full fault tolerance: What no-go theorems rule out and what they don't, arXiv:2501.05694.

[21] A. Zhang et al., Demonstrating quantum error mitigation on logical qubits, arXiv:2501.09079.

[22] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. E. O'Brien, Low-cost error mitigation by symmetry verification, Phys. Rev. A **98**, 062339 (2018).

[23] S. E. Smart and D. A. Mazziotti, Quantum-classical hybrid algorithm using an error-mitigating N-representability condition to compute the Mott metal-insulator transition, Phys. Rev. A **100**, 022517 (2019).

[24] N. N. Bogoliubov, Kinetic equations, J. Phys.-USSR **10**, 265 (1946).

[25] M. Born and H. S. Green, A general kinetic theory of liquids. I. The molecular distribution functions, Proc. R. Soc. London A **188**, 10 (1946).

[26] J. G. Kirkwood, The statistical mechanical theory of transport processes I. General theory, J. Chem. Phys. **14**, 180 (1946).

[27] J. Yvon, The statistical theory of fluids and the equation of state, La théorie statistique des fluides et l'équation d'état, Actuittés scientifiques et industrielles: Hydrodynamique, acoustique: Théories mécaniques (Hermann & cie, Paris, 1935).

[28] S. Chari, R. Inguva, and K. Murthy, A new truncation scheme for BBGKY hierarchy: Conservation of energy and time reversibility, arXiv:1608.02338.

[29] L. Pucci, A. Roy, and M. Kastner, Simulation of quantum spin dynamics by phase space sampling of Bogoliubov-Born-Green-Kirkwood-Yvon trajectories, Phys. Rev. B **93**, 174302 (2016).

[30] R. Paškauskas and M. Kastner, Equilibration in long-range quantum spin systems from a BBGKY perspective, J. Stat. Mech.: Theory Expt. (2012) P02005.

[31] D. Lacroix, Y. Tanimura, S. Ayik, and B. Yilmaz, A simplified BBGKY hierarchy for correlated fermions from a stochastic mean-field approach, Eur. Phys. J. A **52**, 94 (2016).

[32] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X **7**, 021050 (2017).

[33] J. Schwinger, Gauge invariance and Mass. II, Phys. Rev. **128**, 2425 (1962).

[34] T. Angelides, P. Naredi, A. Crippa, K. Jansen, S. Kühn, I. Tavernelli, and D. S. Wang, First-order phase transition of the Schwinger model with a quantum computer, npj Quantum Inf. **11**, 6 (2025).

[35] M. Honda, E. Itou, Y. Kikuchi, L. Nagano, and T. Okuda, Classically emulated digital quantum simulation for screening and confinement in the Schwinger model with a topological term, Phys. Rev. D **105**, 014504 (2022).

[36] G. Pederiva, A. Bazavov, B. Henke, L. Hostetler, D. Lee, H.-W. Lin, and A. Shindler, Quantum state preparation for the schwinger model, PoS LATTICE2021, 047 (2022).

[37] A. Yamamoto, Toward dense QCD in quantum computers, PoS LATTICE2021, 122 (2022).

[38] B. Chakraborty, M. Honda, T. Izubuchi, Y. Kikuchi, and A. Tomiya, Classically emulated digital quantum simulation of the Schwinger model with a topological term via adiabatic state preparation, Phys. Rev. D **105**, 094503 (2022).

[39] D. Ghim and M. Honda, Digital quantum simulation for spectroscopy of schwinger model, PoS LATTICE2023, 213 (2024).

[40] O. Kaikov, T. Saporiti, V. Sazonov, and M. Tamaazousti, Phase diagram of the Schwinger model by adiabatic preparation of states on a quantum simulator, arXiv:2407.09224.
[41] M. D'Anna, M. Krstic Marinkovic, and J. C. P. Barros, Adiabatic state preparation for digital quantum simulations of QED in 1 + 1D, arXiv:2411.01079.

[42] T. Cox and P. C. E. Stamp, Partitioned density matrices and entanglement correlators, Phys. Rev. A **98**, 062110 (2018).

[43] H. F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. **10**, 545 (1959).

[44] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient quantum algorithms for simulating sparse Hamiltonians, Commun. Math. Phys. **270**, 359 (2007).

[45] N. Hatano and M. Suzuki, Finding exponential product formulas of higher orders, Lect. Notes Phys. **679**, 37 (2005).

[46] W. Magnus, On the exponential solution of differential equations for a linear operator, Commun. Pure Appl. Math. **7**, 649 (1954).

[47] We systematically add a small random $O(1/\sqrt{N_S})$ shift to $\varepsilon_{S\eta}$ for the reasons explained in Appendix B.

[48] S. Bernstein, Démonstration du théorème de Weierstrass fondée sur le calcul des probabilités, Commun. Kharkov Math. Soc. **XIII**, 2 (1913).

[49] M. S. Floater, On the convergence of derivatives of Bernstein approximation, J. Approx. Theory **134**, 130 (2005).

[50] J. B. Kogut and L. Susskind, Hamiltonian formulation of Wilson's lattice gauge theories, Phys. Rev. D **11**, 395 (1975).

[51] P. Jordan and E. P. Wigner, About the Pauli exclusion principle, Z. Phys. **47**, 631 (1928).

[52] A. Javadi-Abhari, M. Treinish, K. Krsulich, C. J. Wood, J. Lishman, J. Gacon, S. Martiel, P. D. Nation, L. S. Bishop, A. W. Cross, B. R. Johnson, and J. M. Gambetta, Quantum computing with Qiskit arXiv:2405.08810.

[53] The simulations for *P* were conducted on February 4, 2025, while those for *Q* on February 5, 2025.

[54] T. Saporiti, Dataset for “BBGKY hierarchy for quantum error mitigation”, Zenodo (2025), doi:10.5281/zenodo.16737954.