# Which Electronic Structure Method to Choose in Trajectory Surface Hopping
Dynamics Simulations? Azomethane as a Case Study

Thomas V. Papineau  [  Denis Jacquemin  [  Morgane Vacher  [
<EMAIL_ADDRESS>
###### Abstract

Non-adiabatic dynamics simulations have become a standard approach to explore
photochemical reactions. Such simulations require underlying potential energy
surfaces and couplings between them, calculated at a chosen level of theory,
yet this aspect is rarely assessed. Here, in combination with the popular
trajectory surface hopping dynamics method, we use a high-accuracy XMS-CASPT2
electronic structure level as a benchmark for assessing the performances of
various post-Hartree-Fock methods (namely CIS, ADC(2), CC2 and CASSCF) and
exchange-correlation functionals (PBE, PBE0, CAM-B3LYP) in a TD-DFT/TDA
context, using the isomerization around a double bond as test case. Different
relaxation pathways are identified, and the ability of the different methods
to reproduce their relative importance and timescale is discussed. The results
show that multi-reference electronic structure methods should be preferred,
when studying non-adiabatic decay between excited and ground states. If not
affordable, TD-DFT with TDA and hybrid functionals, and ADC(2) are efficient
alternative, but overestimate the non-radiative decay yield and thus may miss
deexcitation pathways.

UN, Nantes]Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France
UN, Nantes]Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France
[IUF]Institut Universitaire de France (IUF), F-75005 Paris, France UN,
Nantes]Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France

Light-induced phenomena, such as photochemical reactions, have widespread
applications ranging from biology to energy conversion 1, 2, 3, 4, 5, 6. The
photo-induced events involved in these phenomena generally occur on the
femtosecond timescale, making them ultrafast phenomena, a challenge for
experiments. In order to obtain theoretical insights into these processes,
nonadiabatic dynamics simulations are the methods of choice 7, 8, 9. Among
them, the most popular is trajectory surface hopping (TSH)10, 11, an on-the-
fly mixed quantum-classical method. TSH uses the quantum propagation of an
electronic wavepacket, inducing “hops” of independent semiclassical
trajectories between the electronic states. It is an on-the-fly dynamics
method, meaning that the potential energy surfaces (PESs) are generated, using
a selected electronic structure method, as needed as the nuclear trajectories
move along.

During non-radiative relaxations, the molecule encounters regions with
strongly non-adiabatically coupled electronic states, such as the so-called
conical intersections 12, 13, 14, 15. In general, the most suited electronic
structure methods for describing such intersections are multi-configurational
wavefunction methods 16, 17, 18 – although other methods were shown to be
adequate in specific cases19. However, the multi-configurational methods come
with a high computational cost, scaling exponentially with the size of the
active space. This limits the theoretical studies to small molecular systems
(few tens of atoms) and/or short simulation time (up to a ps). To overcome
this bottleneck, the alternative has been to use computationally cheaper
single-reference electronic structure methods such as the second-order
algebraic diagrammatic construction (ADC(2)) 20, 21, 22, 23, 24, 25, or linear
response time-dependent density functional theory (TD-DFT) 20, 26, 27, 28, 29,
30, sometimes with the Tamm-Dancoff approximation (TDA) 31 or its spin-flip
variant 32. Previous works studied the effect of the electronic structure
methods onto the crossing topologies33.

A common approach in the literature to benchmark electronic structure methods
is to compare energies and geometries of key structures or along specific cuts
of the potential energy surfaces19, 34. But globally, the impact of the choice
of the electronic structure methods on the non-adiabatic dynamics has been
little studied35, 36, 37, 38, 39, since simulating and comparing dynamics with
several electronic structure approaches is often computationally beyond reach.
Very recently, the nonadiabatic dynamics of cyclopropanone was extensively
studied with different electronic structure methods and nonadiabatic dynamics
algorithms, showing a much larger sensitivity to the former than the latter38.
The goal of the present Letter is to benchmark the suitability for non-
adiabatic dynamics simulations, in particular internal conversion, of a series
of different electronic structure methods: (extended) multi-state complete
active space second order perturbation theory (MS-CASTP2 40, 41, 42 and XMS-
CASPT2 43), state-averaged complete active space self-consistent field (SA-
CASSCF) 18, 16, second-order approximate coupled cluster singles and doubles
(CC2) 44, ADC(2) 45, 46, configuration interaction singles (CIS) 47, 48, TD-
DFT 49, 50 with several functionals – PBE,51 PBE052 and CAM-B3LYP,53 as
respective representatives of local, global hybrid and range-separated hybrid
models – with and without TDA 54. The formal scaling of the CC2, ADC(2) and
CIS methods is $N^{5}$, with $N$ the number of atoms, while the TD-DFT
approaches formally scale to $N^{4}$. Our aim is to pinpoint a light
electronic structure scheme allowing the study larger molecules of chemical
interest (several hundreds of atoms) and/or on a longer scale (tens to
hundreds of ps).

To assess the performances of a large panel of electronic structure methods,
this Letter focuses on the cis-to-trans photoisomerization of the azomethane
molecule (Scheme Which Electronic Structure Method to Choose in Trajectory
Surface Hopping Dynamics Simulations? Azomethane as a Case Study). This
photochemical reaction is known to occur rapidly in approximately 200 fs upon
light excitation to its first electronic state ($S_{1}$, $n\pi^{*}$) 55, 56.
It proceeds through a rapid rotation about the central CNNC moiety, reaching a
$S_{1}$/$S_{0}$ conical intersection located around CNNC = 90∘, roughly
halfway between the cis and trans isomers, at approximately $t$ = 40 fs.
Previous simulations have made the distinction between “standard”
trajectories, which decay through this first conical intersection, and
“rotator” trajectories which decay through a second symmetry equivalent
intersection reached after a 270∘ total rotation.55, 56 The effect of the
initial conditions on the longer-timescale dynamics was recently investigated
in that system57. Cis-trans isomerization being one of the main types of
photochemical reactions, this photoisomerization reaction provides an ideal
playground for studying the behavior of different electronic structure methods
in describing excited state population decay and geometrical changes. The
small size of the chosen system allows us to use a high-accuracy reference
method accounting for both static and dynamic electronic correlations, namely
XMS-CASPT2, and to propagate the electronic wavefunction with explicitly
calculated non-adiabatic coupling vectors (NACV) (see Methods section for
details).

cis-to-trans photoisomerization of azomethane upon excitation to its first
electronic excited state.

Figures 1a and 1b show the electronic and nuclear dynamics, respectively,
obtained with the reference electronic structure method. Nuclear geometries
are assigned to trans or cis isomers depending on their CNNC dihedral angle:
$>$ or $<$ 90∘, respectively. One can distinguish four time regimes across the
simulation: I) an initial delay of about 30 fs during which all molecules
remain cis and in the first excited state; II) the crossing of a first conical
intersection region up to approximately 60 fs, inducing a non-adiabatic decay
of about 20-25% from $S_{1}$ to $S_{0}$, associated with a complete
isomerization from cis to trans; III) a second plateau in the electronic
populations lasting approximately 40 fs, while most molecules evolve in the
trans configuration; and IV) the crossing of a second strong-coupling region
with further non-adiabatic decay until the end of the simulation, during which
a fraction of the molecules come back to the cis configuration.

Figure 1: Time evolution of (a) electronic populations of $S_{0}$ (blue) and
$S_{1}$ (yellow) states, and (b) cis (teal) and trans (purple) isomer
populations during the XMS-CASPT2 (NACV) reference dynamics. Crashed
trajectories are represented with the red area. The four different time
regimes, namely I, II, III and IV, are delimited by dashed lines.

Figure 2a shows, with a polar graph, the time evolution of the CNNC dihedral
angle – the representative coordinate of the isomerization – with the time
indicated by the increasing distance from the center, for all 100
trajectories. The active state on which each nuclear trajectory is propagated
is indicated by the color of the curve. Analyzing both electronic and
geometrical information in such a way allows identifying the four main
deexcitation pathways presented in Figure 2b: (i) trajectories staying on
$S_{1}$ until the end of the simulation time; (ii) trajectories deexciting
toward $S_{0}$ at the first conical intersection (dihedral of 90∘), ending in
trans; (iii) trajectories staying excited on $S_{1}$ at the first conical
intersection, then decaying toward $S_{0}$ at the second intersection
occurrence (270∘), thus ending in cis; (iv) trajectories deexciting at the
first conical intersection and having enough kinetic energy to cross the
isomerization barrier in the ground state, thus ending in cis. In the language
introduced above, pathways (ii) and (iv) describe “standard” trajectories,
while pathway (iii) describes “rotator” trajectories. In the reference XMS-
CASPT2 dynamics, pathway (i) dominates with almost half of the trajectories;
pathways (ii) and (iii) represent each a bit less than 20% of the
trajectories; pathway (iv) is trifling.

Figure 2: (a) Polar representation of all 100 trajectories of the XMS-CASPT2
(NACV) reference dynamics. For each trajectory, the evolution of the CNNC
dihedral angle is represented by the rotation around the plot, while the
elapsed simulation time is depicted by the increasing distance from the
center. The active state is indicated by the blue ($S_{0}$) or yellow
($S_{1}$) color. Crashes of trajectories are marked with a red cross. (b)
Separation into the four main deexcitation pathways.

Let us now turn our attention to the dynamics obtained with the other methods.
First, staying with CASPT2, two other electronic propagation algorithms are
compared: the Local Diabatization (LD),58 and a simpler energy-gap threshold-
based approach (thrs), in which hops toward $S_{0}$ are enforced when the
$S_{1}-S_{0}$ energy gap goes below 0.1 eV. XMS-CASPT2 with LD was performed
using the BAGEL code59. XMS-CASPT2 and MS-CASPT2 with LD were also achieved
using OpenMolcas60, 61: there, the electronic wavefunction overlap required in
the LD algorithm is calculated approximately, taking into account the CASSCF
part of the CASPT2 wavefunctions only, i.e., neglecting the perturbation part.
It is emphasized that this approximate LD scheme, noted “lLD” hereafter, uses
the correct CASPT2 state ordering and can describe (if necessary, based on the
CASPT2 results) the electronic states as linear combinations of CASSCF
wavefunctions. Despite the development of NACV for methods like TD-DFT62, 63,
64, 65, the threshold-based procedure remains the most commonly used for TD-
DFT and ADC(2) (to describe $S_{1}\rightarrow S_{0}$ transitions). The other
benchmarked electronic structure methods are CASSCF, CC2 and ADC(2), CIS, TD-
PBE, TD-PBE0, TD-CAM-B3LYP, TDA-PBE, TDA-PBE0 and TDA-CAM-B3LYP. Additional
computational details about each method are provided in the Methods section.
Importantly, some electronic structure methods are more prone to failure,
leading to crashes, than others. This typically appears due to convergence
errors in the vicinity of electronic degeneracies. In the present work, all
crashed trajectories are kept and considered as a special case.

Figure 3 compares the time duration of the four different regimes identified
in the reference XMS-CASPT2 (NACV) dynamics (Figure 1) for all benchmarked
methods. All multi-configurational methods, including CASSCF, accurately
reproduce both the four regimes observed in the reference and their duration.
The agreement between NACV and LD schemes was already demonstrated in an
extensive benchmark studies of the different electronic propagation schemes56.
The results of lLD indicate that considering only the CASSCF part of the
CASPT2 wavefunction for calculating overlaps is a reasonable approximation in
the present case, despite a slight underestimation of the duration of the
first decay (regime II). Importantly, the threshold-based XMS-CASPT2 dynamic
shows little difference when compared to the XMS-CASPT2 (NACV) dynamics,
justifying the use of this procedure for CIS and TD-DFT methods.

Figure 3: Durations of the four regimes identified during the non-adiabatic
dynamics calculated with the reference XMS-CASPT2 (NACV) method (reminded at
the top). White areas indicate time windows where less than 5% of the initial
$S_{1}$ population is remaining (because of deexcitations and/or crashes),
making the regime no longer identifiable. Red areas indicate a regime with
more than 50% of the trajectories having crashed.

Regarding the single-reference methods, CC2 appears to be highly unstable near
the $S_{1}/S_{0}$ intersections in azomethane, resulting in systematic
failures in their vicinity. This result is in agreement with a previous study
devoted to the non-adiabatic dynamics among quasi-degenerate excited states in
adenine.20 In contrast, CC2 was able to predict expected characteristics of
the $S_{1}/S_{0}$ conical intersection in a model of rhodopsin protein
chromophore19. In the present work, all tested functionals without applying
TDA, as well as TDA-PBE, show a similarly unstable behavior. This results in
prematurely ended dynamics after regime II for these methods. The use of TDA
proved more resilient and yielded no significant failure. This observation is
in agreement with a previous study combining TSH and TD-DFT.66 ADC(2)
reproduces well the four time regimes and their duration. CIS predicts a
continuous decay in the $S_{1}$ population until the end of the simulation,
although with a lower slope during regime IV. For all tested functionals
(applying or not TDA), the initial delay preceding the first decay (regime I)
is significantly too short – PBE0 performs slightly better in that regard.
This is compensated by a longer regime II, making the end of this regime
appearing at ca. 60 fs irrespective of the chosen electronic structure method
(when there is not too many crashes). Both TDA-PBE0 and TDA-CAM-B3LYP
reproduces well the four different time regimes; the timescales obtained with
TDA-PBE0 are in better agreement with those of the reference. While CASSCF is
the best performing method in reproducing reference’s timescales, ADC(2) is
notably efficient. Moreover, TDA along with an hybrid functional provides
acceptable results despite the too short initial delay.

Figure 4 compares the yield of the non-adiabatic decay by extracting the
electronic population of both states at the limits of the regimes I-IV
identified with the reference XMS-CASPT2 (NACV) results. The exact times
limiting the different regimes depend on the electronic structure method (see
Figure 3). In particular, CIS and ADC(2) (thrs) have no regime III. The decay
yield is well reproduced by XMS-CASPT2 with LD and lLD electronic propagation
schemes. It is slightly faster in the case of the threshold-approach. This can
be explained: the threshold-based hopping procedure does not allow
trajectories to hop back to the excited state (back-hopping phenomenon). As
seen with the timescale analysis, most of the trajectories crashed with CC2,
PBE (with and without TDA) and with the hybrid functionals without TDA. All
other tested electronic structure methods present larger decays than the
reference one. While the use of the threshold-approach can explain part of the
too important decay, it cannot account for all of it. The too fast decay
predicted here in azomethane by TDA-PBE0 and TDA-CAM-B3LYP contrasts with the
failure of TD-DFT to describe the non-adiabatic decay in adenine studied in
previous works.35, 20 Ignoring MS-CASPT2, CASSCF is here again the best
performing method in comparison to XMS-CASPT2. While showing different amounts
of failure, the single-reference methods ADC(2) (LD), CIS, TDA-PBE0 and TDA-
CAM-B3LYP predict similarly too large non-adiabatic decay yields.

Figure 4: Electronic state populations ($S_{0}$ in blue, $S_{1}$ in yellow) at
the end of each identified time regime, obtained with all electronic structure
methods. Crashed trajectories are represented by a red area.

Figure 5 presents the relative importance of the main relaxation pathways
identified in the reference dynamics (Figure 2) for each studied electronic
structure method. It is noted that for some methods, mainly DFT-based ones, a
significant amount of trajectories do not follow any of the four typical
relaxations pathways and thus remain unsorted. All multi-reference
wavefunctions methods reproduce qualitatively the main relaxation pathways,
with MS-CASPT2 and CASSCF overestimating the importance of the “standard”
trajectories, that is, pathways (ii) and (iv). This is in agreement with the
too efficient electronic population decay through the first encountered
conical intersection (Figure 4). This feature is even more present in the
single-reference post-HF wavefunction and DFT methods: with them, the pathways
that stay on the excited state after the first crossing are incorrectly
negligible and pathway (ii) dominates largely the dynamics. One could be
tempted to assign this result to the threshold approach. However, the XMS-
CASPT2 (thrs) results indicate that it is not the main reason and that this
feature is rather an intrinsic failure of the single-reference approaches.

Figure 5: Relative importance of the four identified relaxation pathways
(i-iv) obtained with the different electronic structure methods.

In this Letter, we have tested a range of different electronic structure
methods for simulating the cis-to-trans photoisomerisation reaction of
azomethane in combination with the surface hopping dynamics method: multi-
reference (XMS-CASPT2, MS-CASPT2 and CASSCF) and single-reference (ADC(2),
CC2, CIS) wavefunction methods as well as TD-DFT based approaches (PBE, PBE0
and CAM-B3LYP functionals with and without TDA). Different electronic
propagation and hopping algorithms were considered: using NACV, LD or the
energy-threshold approach. During the first 150 fs after photoexcitation to
the $S_{1}$ state, the molecule has time to reach the first 90∘ conical
intersection after ca. 40 fs, and the 270∘ conical intersection from 120 fs
and onwards. Analyzing the electronic and geometric time evolutions, four
different time regimes and relaxation pathways could be identified in the
reference XMS-CASPT2 (NACV) dynamics.

Multi-reference methods reproduce closely the timescale and yield of the
induced non-adiabatic dynamics, as well as the different relaxation pathways.
CASSCF slightly overestimates the decay through the first conical
intersection, and thus the amount of “standard” trajectories and corresponding
pathways. CC2, TD-PBE, TDA-PBE, TD-PBE0 and TD-CAM-B3LYP appear to be too
unstable in the vicinity of electronic degeneracies and thus, no trustworthy
statistical results could be obtained. ADC(2) (LD), TDA-PBE0 and TDA-CAM-B3LYP
qualitatively reproduce the timescale of the decay (despite crashes of up to
40% of the trajectories) but predicts an almost total decay through the first
crossing and thus poorly describe the “rotator” trajectories. CIS, in addition
to underestimating the “rotator” trajectories, does not reproduce the plateau
in the electronic population evolution after crossing the first intersection.
The too important non-adiabatic decay obtained with these methods cannot be
attributed to the threshold approach.

In short, the present trajectory surface hopping work on azomethane shows that
multi-reference electronic structure methods are the preferred choices to
simulate non-adiabatic dynamics. If not affordable, ADC(2), TDA-PBE0 and TDA-
CAM-B3LYP are alternatives which reproduce well the timescale of the reference
dynamics but overestimate the non-adiabatic decay yield when going through a
crossing and thus may miss some deexcitation pathways. Interestingly, the
suitability of the different electronic structure methods could hardly be
predicted from a static analysis of the potential energy surfaces (see SI for
energetics at Franck-Condon geometry and scan along the CNNC dihedral angle
coordinate). This suggests that the common approach to benchmark and validate
electronic structure methods before simulating dynamics may not be sufficient,
as concluded in a recent work38. Of course, the present conclusions were
obtained in the specific case of the cis-to-trans photoisomerisation reaction
of azomethane upon $S_{1}$ excitation, and thus cannot be assumed to be
general for all molecules, all electronic states, nor all types of photo-
induced processes.

## Methods

The surface hopping simulations were carried out with the SHARC code67, 68
using the fewest switches algorithm: 100 trajectories were simulated for 150
fs with a time step of 0.5 fs. Initial conditions for the different sets of
100 trajectories were generated with SHARC tools, sampling geometries and
velocities from the Wigner distribution (without temperature broadening). This
was done by using harmonic frequencies calculated at the relevant optimized
ground-state geometry using the respective electronic structure level. The
initial conditions are thus in principle not identical for the various tested
electronic structure methods, although in practice no significant difference
was observed (see SI). This choice was made in order to reflect a realistic
non-adiabatic study. We, nevertheless, simulated dynamics at ADC(2) and TDA-
PBE0 levels of theory with initial conditions generated with XMS-CASPT2 method
(see SI): the results show that the initial conditions have little effect (in
contradiction with a previous study on that system57) and confirm that the
difference come from the electronic structure. To conserve the total energy
after a hop, the full velocity vector was rescaled. An energy based
decoherence correction was applied with a decay factor of 0.1 Hartree69. For
all electronic structure methods benchmarked herein, the selected basis set is
_def2_ -SVP70.

The reference XMS-CASPT2 (NACV) calculations have been achieved with the BAGEL
code59 and a modified version of the SHARC-BAGEL interface, using an active
space consisting of 6 electrons in 4 orbitals: the $\pi$, $\pi^{*}$ and the
two N lone pairs. No shift was used, and IPEA was deactivated. BAGEL was also
used to perform the XMS-CASPT2 (LD) simulations. XMS-CASPT2 (lLD), MS-CASPT2
and CASSCF simulations were performed with OpenMolcas60, 61, using the same
active space and RI approximation with analytical Cholesky decomposition. CC2
and ADC(2) electronic structure calculations have been done with the Turbomole
7.3 code71, using the defaults parameters requested by the SHARC-RICC2
interface (RI approximation with default selection of auxilliary basis set,
and frozen core orbitals). Finally, CIS, TD-DFT and TDA-DFT calculations were
made with Gaussian 16 rev. A03,72 using the tight SCF convergence criteria and
fine integration grid. CIS calculations have been achieved by requesting TDA-
HF calculations. The TD- and TDA-PBE0 dynamics have also been performed using
Turbomole 7.3 code71 (using a modified version of the SHARC-RICC2 interface)
and Orca 5.0.373, 74. In Turbomole, the RI approximation was activated and the
m3 integration grid was used. Within Orca, the RIJCOSX approximation was
enabled, along with the _def2_ /J and _def2_ -SVP/C auxilliary basis sets. The
tight SCF convergence criteria was requested and the defgrid2 integration grid
was selected.

The authors thank the ANR for financial support in the framework of the
BiBiFlu project. The project is also partly funded by the European Union
through ERC grant 101040356 (M.V.). Views and opinions expressed are however
those of the author(s) only and do not necessarily reflect those of the
European Union or the European Research Council Executive Agency. Neither the
European Union nor the granting authority can be held responsible for them.
The authors thank the CCIPL/Glicid mesocenter installed in Nantes and GENCI-
IDRIS (Grant 2021-101353) for generous allocation of computational time.

Active space orbitals. Comparison of initial internal coordinates for all
electronic structure methods. Test of the initial conditions on the dynamics.
Total energy variations along the dynamics for all benchmarked electronic
structure methods. Time evolutions of electronic populations and isomer
populations defining the different time regimes for all benchmarked electronic
structure methods. Energetics at Franck-Condon geometry and scan along the
CNNC dihedral angle coordinate, with all tested electronic structure methods.
The modified SHARC interfaces can be made available upon reasonable request.

## References

  * Wald 1968 Wald, G. Molecular Basis of Visual Excitation. _Science_ 1968, _162_ , 230–239
  * Liu et al. 1990 Liu, Z. F.; Hashimoto, K.; Fujishima, A. Photoelectrochemical Information Storage Using an Azobenzene Derivative. _Nature_ 1990, _347_ , 658–660
  * Polosukhina et al. 2012 Polosukhina, A.; Litt, J.; Tochitsky, I.; Nemargut, J.; Sychev, Y.; Kouchkovsky, I. D.; Huang, T.; Borges, K.; Trauner, D.; Gelder, R. N. V. et al. Photochemical Restoration of Visual Responses in Blind Mice. _Neuron_ 2012, _75_ , 271–282
  * Zhang et al. 2013 Zhang, J.; Zou, Q.; Tian, H. Photochromic Materials: More Than Meets The Eye. _Adv. Mater._ 2013, _25_ , 378–399
  * Richards 2016 Richards, V. Molecular Machines. _Nature Chem._ 2016, _8_ , 1090–1090
  * Baroncini et al. 2020 Baroncini, M.; Silvi, S.; Credi, A. Photo- and Redox-Driven Artificial Molecular Motors. _Chem. Rev._ 2020, _120_ , 200–268
  * Vacher et al. 2016 Vacher, M.; Bearpark, M. J.; Robb, M. A. Direct Methods for Non-Adiabatic Dynamics: Connecting the Single-Set Variational Multi-Configuration Gaussian (vMCG) and Ehrenfest Perspectives. _Theor. Chem. Acc._ 2016, _135_ , 1–11
  * Curchod and Martínez 2018 Curchod, B. F. E.; Martínez, T. J. Ab Initio Nonadiabatic Quantum Molecular Dynamics. _Chem. Rev._ 2018, _118_ , 3305–3336
  * Crespo-Otero and Barbatti 2018 Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics. _Chem. Rev._ 2018, _118_ , 7026–7068
  * Tully and Preston 1971 Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H+ with D2. _J. Chem. Phys._ 1971, _55_ , 562–572
  * Tully 1990 Tully, J. C. Molecular Dynamics with Electronic Transitions. _J. Chem. Phys._ 1990, _93_ , 1061–1071
  * Teller 1937 Teller, E. The Crossing of Potential Surfaces. _J. Phys. Chem._ 1937, _41_ , 109–116
  * Robb et al. 1995 Robb, M. A.; Bernardi, F.; Olivucci, M. Conical Intersections as a Mechanistic Feature of Organic Photochemistry. _Pure Appl. Chem._ 1995, _67_ , 783–789
  * Yarkony 2001 Yarkony, D. R. Conical Intersections: The New Conventional Wisdom. _J. Phys. Chem. A_ 2001, _105_ , 6277–6293
  * Domcke and Yarkony 2012 Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. _Ann. Rev. Phys. Chem._ 2012, _63_ , 325–352
  * Roos et al. 1980 Roos, B. O.; Taylor, P. R.; Sigbahn, P. E. A Complete Active Space SCF Method (CASSCF) using a Density Matrix Formulated Super-CI Approach. _Chem. Phys._ 1980, _48_ , 157–173
  * Yeager and Jor̈gensen 2008 Yeager, D. L.; Jor̈gensen, P. Convergency Studies of Second and Approximate Second Order Multiconfigurational Hartree-Fock Procedures. _J. Chem. Phys._ 2008, _71_ , 755–760
  * Olsen 2011 Olsen, J. The CASSCF method: A Perspective and Commentary. _Int. J. Quantum Chem._ 2011, _111_ , 3267–3272
  * Tuna et al. 2015 Tuna, D.; Lefrancois, D.; Wolański, Ł.; Gozem, S.; Schapiro, I.; Andruniów, T.; Dreuw, A.; Olivucci, M. Assessment of Approximate Coupled-Cluster and Algebraic-Diagrammatic-Construction Methods for Ground- and Excited-State Reaction Paths and the Conical-Intersection Seam of a Retinal-Chromophore Model. _J. Chem. Theory Comput._ 2015, _11_ , 5758–5781
  * Plasser et al. 2014 Plasser, F.; Crespo-Otero, R.; Pederzoli, M.; Pittner, J.; Lischka, H.; Barbatti, M. Surface Hopping Dynamics with Correlated Single-Reference Methods: 9H-Adenine as a Case Study. _J. Chem. Theory Comput._ 2014, _10_ , 1395–1405
  * Xie et al. 2019 Xie, W.; Sapunar, M.; Došlić, N.; Sala, M.; Domcke, W. Assessing the Performance of Trajectory Surface Hopping Methods: Ultrafast Internal Conversion in Pyrazine. _J. Chem. Phys._ 2019, _150_ , 154119
  * Braun et al. 2022 Braun, G.; Borges, J., Itamar; Aquino, A. J. A.; Lischka, H.; Plasser, F.; do Monte, S. A.; Ventura, E.; Mukherjee, S.; Barbatti, M. Non-Kasha Fluorescence of Pyrene Emerges from a Dynamic Equilibrium Between Excited States. _J. Chem. Phys._ 2022, _157_ , 154305
  * Mansour et al. 2022 Mansour, R.; Mukherjee, S.; Pinheiro, M.; Noble, J. A.; Jouvet, C.; Barbatti, M. Pre-Dewar Structure Modulates Protonated Azaindole Photodynamics. _Phys. Chem. Chem. Phys._ 2022, _24_ , 12346–12353
  * Mansour et al. 2022 Mansour, R.; Toldo, J. M.; Barbatti, M. Role of the Hydrogen Bond on the Internal Conversion of Photoexcited Adenosine. _J. Phys. Chem. Lett._ 2022, _13_ , 6194–6199
  * Ferté et al. 2023 Ferté, A.; Houssin, A.; Albouy, N.; Merritt, I. C. D.; Vacher, M. ESIPT in the Pyrrol Pyridine Molecule: Mechanism, Timescale and Yield Revealed using Dynamics Simulations. _Phys. Chem. Chem. Phys._ 2023, _25_ , 9761–9765
  * Ye et al. 2020 Ye, L.; Xu, C.; Gu, F. L.; Zhu, C. Functional and Basis Set Dependence for Time-Dependent Density Functional Theory Trajectory Surface Hopping Molecular Dynamics: Cis-Azobenzene Photoisomerization. _J. Comput. Chem._ 2020, _41_ , 635–645
  * Cao and Chen 2020 Cao, J.; Chen, D.-c. The Excited-State Relaxation Mechanism of Potential UVA-Activated Phototherapeutic Molecules: Trajectory Surface Hopping Simulations of both 4-Thiothymine and 2,4-Dithiothymine. _Phys. Chem. Chem. Phys._ 2020, _22_ , 10924–10933
  * Pápai et al. 2021 Pápai, M.; Li, X.; Nielsen, M. M.; Møller, K. B. Trajectory Surface-Hopping Photoinduced Dynamics From Rydberg States of Trimethylamine. _Phys. Chem. Chem. Phys._ 2021, _23_ , 10964–10977
  * Liu and Zhu 2021 Liu, Y.; Zhu, C. Trajectory Surface Hopping Molecular Dynamics Simulations for Retinal Protonated Schiff-Base Photoisomerization. _Phys. Chem. Chem. Phys._ 2021, _23_ , 23861–23874
  * Zobel et al. 2023 Zobel, J. P.; Kruse, A.; Baig, O.; Lochbrunner, S.; Bokarev, S. I.; Kühn, O.; González, L.; Bokareva, O. S. Can Range-Separated Functionals be Optimally Tuned to Predict Spectra and Excited State Dynamics in Photoactive Iron Complexes? _Chem. Sci._ 2023, _14_ , 1491–1502
  * Talotta et al. 2020 Talotta, F.; Boggio-Pasqua, M.; González, L. Early Relaxation Dynamics in the Photoswitchable Complex trans-[RuCl(NO)(py)4]2+. _Chem. Eur. J._ 2020, _26_ , 11522–11528
  * Minezawa and Nakajima 2019 Minezawa, N.; Nakajima, T. Trajectory Surface Hopping Molecular Dynamics Simulation by Spin-Flip Time-Dependent Density Functional Theory. _J. Chem. Phys._ 2019, _150_ , 204120
  * Gozem et al. 2014 Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; Huix-Rotllant, M.; Ferré, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A. et al. Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection. _J. Chem. Theory Comput._ 2014, _10_ , 3074–3084
  * Marsili et al. 2021 Marsili, E.; Prlj, A.; Curchod, B. F. E. Caveat when using ADC(2) for studying the photochemistry of carbonyl-containing molecules. _Phys. Chem. Chem. Phys._ 2021, _23_ , 12945–12949
  * Barbatti et al. 2012 Barbatti, M.; Lan, Z.; Crespo-Otero, R.; Szymczak, J. J.; Lischka, H.; Thiel, W. Critical appraisal of excited state nonadiabatic dynamics simulations of 9H-adeninea). _J. Chem. Phys._ 2012, _137_ , 22A503
  * Mai et al. 2019 Mai, S.; Atkins, A. J.; Plasser, F.; González, L. The Influence of the Electronic Structure Method on Intersystem Crossing Dynamics. The Case of Thioformaldehyde. _J. Chem. Theory Comput._ 2019, _15_ , 3470–3480
  * Chakraborty et al. 2021 Chakraborty, P.; Liu, Y.; McClung, S.; Weinacht, T.; Matsika, S. Time Resolved Photoelectron Spectroscopy as a Test of Electronic Structure and Nonadiabatic Dynamics. _J. Phys. Chem. Lett._ 2021, _12_ , 5099–5104
  * Janoš and Slavíček 2023 Janoš, J.; Slavíček, P. What Controls the Quality of Photodynamical Simulations? Electronic Structure Versus Nonadiabatic Algorithm. _J. Chem. Theory Comput._ 2023, _19_ , 8273–8284
  * Barneschi et al. 2023 Barneschi, L.; Kaliakin, D.; Huix-Rotllant, M.; Ferré, N.; Filatov(Gulak), M.; Olivucci, M. Assessment of the Electron Correlation Treatment on the Quantum-Classical Dynamics of Retinal Protonated Schiff Base Models: XMS-CASPT2, RMS-CASPT2, and REKS Methods. _J. Chem. Theory Comput._ 2023, _19_ , 8189–8200
  * Andersson et al. 1990 Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-Order Perturbation Theory with a CASSCF Reference Function. _J. Phys. Chem._ 1990, _94_ , 5483–5488
  * Andersson et al. 1992 Andersson, K.; Malmqvist, P.; Roos, B. O. Second-Order Perturbation Theory with e Complete Active Space Self-Consistent Field Reference Function. _J. Chem. Phys._ 1992, _96_ , 1218–1226
  * Finley et al. 1998 Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. The Multi-State CASPT2 Method. _Chem. Phys. Lett._ 1998, _288_ , 299–306
  * Shiozaki et al. 2011 Shiozaki, T.; Győrffy, W.; Celani, P.; Werner, H.-J. Extended Multi-State Complete Active Space Second-Order Perturbation Theory: Energy and Nuclear Gradients. _J. Chem. Phys._ 2011, _135_ , 081106
  * Christiansen et al. 1995 Christiansen, O.; Koch, H.; Jørgensen, P. The Second-Order Approximate Coupled Cluster Singles and Doubles Model CC2. _Chem. Phys. Lett._ 1995, _243_ , 409–418
  * Schirmer 1982 Schirmer, J. Beyond the Random-Phase Approximation: A new Approximation Scheme for the Polarization Propagator. _Phys. Rev. A_ 1982, _26_ , 2395–2416
  * Schirmer and Trofimov 2004 Schirmer, J.; Trofimov, A. B. Intermediate State Representation Approach to Physical Properties of Electronically Excited Molecules. _J. Chem. Phys._ 2004, _120_ , 11449–11464
  * Bene et al. 1971 Bene, J. E. D.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. X. Molecular Orbital Studies of Excited States with Minimal and Extended Basis Sets. _J. Chem. Phys._ 1971, _55_ , 2236–2241
  * Foresman et al. 1992 Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. Toward a Systematic Molecular Orbital Theory for Excited States. _J. Phys. Chem._ 1992, _96_ , 135–149
  * Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. _Phys. Rev. Lett._ 1984, _52_ , 997–1000
  * CASIDA 1995 CASIDA, M. E. _Recent Advances in Density Functional Methods_ ; Recent Advances in Computational Chemistry 0; WORLD SCIENTIFIC, 1995; Vol. Volume 1; pp 155–192
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Adamo and Barone 1999 Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: the PBE0 Model. _J. Chem. Phys._ 1999, _110_ , 6158–6170
  * Yanai et al. 2004 Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). _Chem. Phys. Lett._ 2004, _393_ , 51–56
  * Hirata and Head-Gordon 1999 Hirata, S.; Head-Gordon, M. Time-Dependent Density Functional Theory within the Tamm–Dancoff Approximation. _Chem. Phys. Lett._ 1999, _314_ , 291–299
  * 55 Sellner, B.; Ruckenbauer, M.; Stambolić, I.; Barbatti, M.; Aquino, A. J. A.; Lischka, H. Photodynamics of Azomethane: A Nonadiabatic Surface-Hopping Study. _114_ , 8778–8785
  * Merritt et al. 2023 Merritt, I. C. D.; Jacquemin, D.; Vacher, M. Nonadiabatic Coupling in Trajectory Surface Hopping: How Approximations Impact Excited-State Reaction Dynamics. _J. Chem. Theory Comput._ 2023, _19_ , 1827–1842
  * Pieroni et al. 2023 Pieroni, C.; Becuzzi, F.; Creatini, L.; Granucci, G.; Persico, M. Effect of Initial Conditions Sampling on Surface Hopping Simulations in the Ultrashort and Picosecond Time Range. Azomethane Photodissociation as a Case Study. _J. Chem. Theory Comput._ 2023, _19_ , 2430–2445
  * Granucci et al. 2001 Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. _J. Chem. Phys._ 2001, _114_ , 10608–10615
  * 59 BAGEL, Brilliantly Advanced General Electronic-structure Library. http://www.nubakery.org under the GNU General Public License.
  * Fdez. Galván et al. 2019 Fdez. Galván, I.; Vacher, M.; Alavi, A.; Angeli, C.; Aquilante, F.; Autschbach, J.; Bao, J. J.; Bokarev, S. I.; Bogdanov, N. A.; Carlson, R. K. et al. OpenMolcas: From Source Code to Insight. _J. Chem. Theory Comput._ 2019, _15_ , 5925–5964
  * Li Manni et al. 0 Li Manni, G.; Fdez. Galván, I.; Alavi, A.; Aleotti, F.; Aquilante, F.; Autschbach, J.; Avagliano, D.; Baiardi, A.; Bao, J. J.; Battaglia, S. et al. The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry. _J. Chem. Theory Comput._ 0, _0_ , null
  * Send and Furche 2010 Send, R.; Furche, F. First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance. _J. Chem. Phys._ 2010, _132_ , 044107
  * Tavernelli et al. 2010 Tavernelli, I.; Curchod, B. F. E.; Laktionov, A.; Rothlisberger, U. Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm–Dancoff approximation and beyond. _J. Chem. Phys._ 2010, _133_ , 194104
  * Li and Liu 2014 Li, Z.; Liu, W. First-order nonadiabatic coupling matrix elements between excited states: A Lagrangian formulation at the CIS, RPA, TD-HF, and TD-DFT levels. _J. Chem. Phys._ 2014, _141_ , 014110
  * Ou et al. 2015 Ou, Q.; Bellchambers, G. D.; Furche, F.; Subotnik, J. E. First-order derivative couplings between excited states from adiabatic TDDFT response theory. _J. Chem. Phys._ 2015, _142_ , 064114
  * Tapavicza et al. 2008 Tapavicza, E.; Tavernelli, I.; Rothlisberger, U.; Filippi, C.; Casida, M. E. Mixed time-dependent density-functional theory/classical trajectory surface hopping study of oxirane photochemistry. _J. Chem. Phys._ 2008, _129_ , 124108
  * Mai et al. 2023 Mai, S.; Avagliano, D.; Heindl, M.; Marquetand, P.; Menger, M. F. S. J.; Oppel, M.; Plasser, F.; Polonius, S.; Ruckenbauer, M.; Shu, Y. et al. SHARC3.0: Surface Hopping Including Arbitrary Couplings — Program Package for Non-Adiabatic Dynamics. https://sharc-md.org/, 2023
  * Mai et al. 2018 Mai, S.; Marquetand, P.; González, L. Nonadiabatic Dynamics: The SHARC Approach. _WIREs Comput. Mol. Sci._ 2018, _8_ , e1370
  * Granucci and Persico 2007 Granucci, G.; Persico, M. Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping. _J. Chem. Phys._ 2007, _126_
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * 71 TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com.
  * Frisch et al. 2016 Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H. et al. Gaussian 16, Revision C.01. 2016; Gaussian Inc. Wallingford CT
  * Neese et al. 2020 Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The ORCA quantum chemistry program package. _J. Chem. Phys._ 2020, _152_ , 224108
  * Neese 2022 Neese, F. Software update: The ORCA program system—Version 5.0. _WIREs Comput. Mol. Sci._ 2022, _12_ , e1606


# Scalable mitigation of measurement errors on quantum computers

Paul D. Nation<EMAIL_ADDRESS>Hwajung Kang  Neereja Sundaresan  Jay M.
Gambetta  IBM Quantum, Yorktown Heights, NY 10598 USA

###### Abstract

We present a method for mitigating measurement errors on quantum computing
platforms that does not form the full assignment matrix, or its inverse, and
works in a subspace defined by the noisy input bit-strings. This method
accommodates both uncorrelated and correlated errors, and allows for computing
accurate error bounds. Additionally, we detail a matrix-free preconditioned
iterative solution method that converges in $\mathcal{O}(1)$ steps that is
performant and uses orders of magnitude less memory than direct factorization.
We demonstrate the validity of our method, and mitigate errors in a few
seconds on numbers of qubits that would otherwise be intractable.

##  I Introduction

Recently, rapid developments in the fabrication, control, and deployment of
quantum computing systems has brought qubit counts to $\sim 100$, where it
might be possible to show advantage over classical computation methods in one
or more limited cases [1, 2, 3]. However, such breakthroughs are hampered by
noise and errors that conspire to limit the effectiveness of quantum computers
at tackling problems of appreciable scale. To counteract these effects,
researchers have turned to mitigation methods that approximately cancel
quantum gate [4, 5, 6, 7, 8, 9, 10] and measurement assignment [11, 12, 13,
14, 15, 16, 17, 18, 19] errors. For short-depth quantum circuits that can be
executed on current generation hardware, measurement errors play an out-sized
role, and their correction is critical to many near-term experiments [20, 21,
22, 23, 24, 25, 26, 27].

In the canonical situation where initialization noise is minimal, measurement
errors over $N$-qubits can be treated classically and satisfy

$\vec{p}_{\rm noisy}=A\vec{p}_{\rm ideal},$ (1)

where $\vec{p}_{\rm noisy}$ is a vector of noisy probabilities returned by the
quantum system, $\vec{p}_{\rm ideal}$ is the probabilities in absence of
measurement errors (but still includes e.g. gate errors), and $A$ is the
$2^{N}\times 2^{N}$ complete assignment matrix (A-matrix) where element
$A_{\tt row,\tt col}$ is the probability of bit-string $\tt col$ being
converted to bit-string $\tt row$ by the measurement error process [see App.
(A) for examples]. While computing $A$ requires executing $2^{N}$ circuits, it
is often the case that errors on multiple qubits can be well approximated
using at most $\mathcal{O}(N)$ calibration circuits; the A-matrix can be
approximated efficiently.

Equation (1) has a solution $\vec{p}_{\rm ideal}$ readily found using direct
LU-factorization. However, direct methods necessarily generate quasi-
probability distributions due to finite sampling that contain negative values,
but still sum to one, that are incompatible with the requirement of
$\vec{p}_{\rm ideal}$ being a probability vector. Consequently, a bounded-
minimization approach solving $||A\vec{p}_{\rm ideal}-\vec{p}_{\rm
noisy}||^{2}_{2}$, where $\vec{p}_{\rm ideal}$ is constrained to be positive,
is often used in place of a direct solution [12, 13, 14, 16, 28]. Although
physically appealing, the run times of these methods are orders of magnitude
longer than those of direct techniques. Alternatively, it has been shown that
quasi-probabilities can be used provided that one mitigates expectation values
[29, 4, 11]. As proven in Ref. [11], these quasi-probabilities provide an
unbiased estimate for the expectation value $\xi$ of an operator $O$, with a
spectral radius of one, that is diagonal in the computational basis

$\xi=\sum_{i=0}^{2^{N}-1}\left[OA^{-1}\vec{p}_{\rm noisy}\right]_{i}.$ (2)

Near-term algorithms such as the ubiquitous Variational Quantum Eigensolver
(VQE) [30, 20] and quantum machine learning [22, 27] rely on the computation
of expectation values, making the correction of measurement errors in these
quantities an important step along the road to quantum advantage.

System | Avg. assignment error (%)
---|---
Alibaba 11Q [31] | 7.4
Google 53Q Sycamore [32] | 3.2
IBM 27Q Falcon_R5.11 [33] | 1.1
IONQ 11Q [34] | 0.4
Quantum Inspire 5Q Starmon [35] | 4.0
Rigetti 32Q Aspen-9 [34] | 6.1
Table 1: Representative error rates for cloud-accessible quantum computing
systems.

Current measurement mitigation techniques utilize the full $2^{N}$-dimensional
probability space, and thus do not scale beyond a handful of qubits. A
truncation scheme was developed in [23], however it did so at the loss of
measurement information, and still required explicit construction of the full
A-matrix. Creating a scalable mitigation strategy requires reducing the
dimensionality of the linear system in Eq. (1) without the need for computing
$A$ itself. Fortunately, present day cloud-accessible quantum computing
systems have measurement error rates of a few-percent or less, see Table 1,
indicating it is possible to view the measurement error process as a small
correction to the ideal probability distribution; measurement errors
redistribute small fractions of probability from a given bit-string primarily
to those that are a short Hamming distance away [e.g. see App. (A)]. To good
approximation, the solution is contained within $\vec{p}_{\rm noisy}$, and we
can mitigate errors in a re-normalized subspace defined by these bit-strings.
Worst case, this subspace dimension is equal to the number of times the input
circuit is sampled. For cloud-accessible quantum computers, the number of
times a circuit can be sampled is limited, typically $8192$ times on IBM
Quantum systems, and therefore the dimensionality of the corresponding reduced
assignment matrix $\tilde{A}$ can be markedly smaller than the full A-matrix
for $N$-qubits. In practice $\tilde{A}$ is often small enough such that the
solution is amenable to standard LU-factorization, returning a vector of
quasi-probabilities for use in a reduced version of Eq. (2).

However, for situations where explicitly forming $\tilde{A}$ is still
prohibitive due to large numbers of unique samples, it is possible to use
preconditioned matrix-free iterative linear solution methods. Such methods
have also been explored in numerical solutions of large-scale steady-state
density matrices [36]. In practice this gives quick convergence, typically in
$\mathcal{O}(1)$ steps, and is competitive with direct solution run times
while requiring orders of magnitude less memory. Although the methods
introduced here return quasi-probabilities, it is possible to find the nearest
probability distribution, in terms of $L2$-norm, in linear time [37].

In this paper, we detail this efficient mitigation method beginning with Sec.
(II) that motives the subspace reduction procedure, and describes how it is
performed. Section (III) shows how preconditioned matrix-free methods can be
utilized for a performant and memory efficient solution technique. In Sec.
(IV) we show that one can obtain bounds on the variance of the computed
expectation values in a similarly efficient manner with an overhead of
$\mathcal{O}(1)$ in terms of additional run time. We demonstrate our technique
in Sec. (V), showing the validity of our method, and mitigating readout errors
out to numbers of qubits that would be intractable on even the largest of
supercomputers using previous methods.

##  II Subspace reduction

Figure 1: a) Simulated discrete probability distribution for 30 samples
showing the ideal distribution for a Greenberger–Horne–Zeilinger (GHZ) state
subject to gate errors only as well as noisy data affected by both gate and
measurement errors. Observed bit-strings are in bold. The same random number
seed was used for both simulations such that the difference between ideal and
noisy distributions is solely due to measurement errors. b) Elements of the
assignment matrix $A$ that are used for constructing the reduced assignment
matrix $\tilde{A}$ that corresponds to the noisy distribution from part (a).
Columns of $\tilde{A}$ must be renormalized such that they again sum to one,
Eq. (5).

We aim to construct $\tilde{A}$ without necessarily forming the full
$A$-matrix. To this end we look to compute elements $A_{\tt row,\tt col}$,
directly from the bit-string values of the input noisy counts, and a small set
of calibration data matrices. For concreteness, we study the case of tensored
measurement errors as IBM Quantum systems are calibrated for high-fidelity
quantum non-demolition (QND) measurements where uncorrelated errors are
nominally dominant. Other quantum hardware vendors also report the same [26].
The full-dimensional tensored A-matrix $A^{(T)}$ over $N$-qubits can be
constructed from $N$ $2\times 2$ calibration matrices:
$A^{(T)}=S_{N-1}\otimes\dots S_{1}\otimes S_{0}$ where $S_{k}$ is the
calibration matrix for the $k$-th qubit with the form

$S_{k}=\begin{bmatrix}P^{(k)}_{0,0}&P^{(k)}_{0,1}\\\
P^{(k)}_{1,0}&P^{(k)}_{1,1}\end{bmatrix},$ (3)

where $P^{k}_{i,j}$ is the probability of the $k$-th qubit being in state
$j\in\\{0,1\\}$ and measured in state $i\in\\{0,1\\}$. Here we use the
convention that qubit $0$ corresponds to the least-significant bit. For two
bit-strings: $\texttt{row},\texttt{col}\in\\{0,1\\}^{N}$, the matrix element
$A^{(T)}_{\tt row,\tt col}$ can be computed using:

$A^{(T)}_{\texttt{row},\texttt{col}}=\prod_{k=0}^{N-1}P^{(k)}_{\texttt{row}[N-1-k],\texttt{col}[N-1-k]}.$
(4)

Therefore it is possible to compute individual matrix elements directly from
bit-string values and a number of calibration matrices that scales at most
linearly with the number of qubits. Accommodating correlated errors in our
method simply amounts to finding an equivalent expression to Eq. (4). We give
an example in App. (B), where we intentionally induce correlated errors into
the measurement process. As a corollary of grabbing elements individually, it
is possible to select only those elements within a given Hamming distance,
$d(\texttt{row},\texttt{col})\leq D$, where $D$ is the desired maximum
distance. This allows for varying the sparsity of $\tilde{A}$, and examining
the effect of low-distance approximations.

$\tilde{A}$ is defined to be the assignment matrix over only those bit-strings
observed in $\vec{p}_{\rm noisy}$. To understand the validity of this
reduction consider the simulation presented in Fig. (1a) where we compare the
distributions for $\vec{p}_{\rm ideal}$ and $\vec{p}_{\rm noisy}$. In effect,
measurement errors take blocks of probability, as defined by the finite number
of samples, from $\vec{p}_{\rm ideal}$ and redistribute them amongst other
bit-strings to get $\vec{p}_{\rm noisy}$. This redistribution is predominantly
to those bit-strings that are close in Hamming distance. Importantly, as seen
in Fig. (1a), when measurement errors are small, and a sufficient number of
circuit samples has been performed, it is unlikely that a given bit-string in
$\vec{p}_{\rm ideal}$ is completely redistributed; the bit-strings in
$\vec{p}_{\rm ideal}$ are also contained in $\vec{p}_{\rm noisy}$. Therefore,
mitigating measurement errors requires only elements $A_{\tt row,\tt col}$
that correspond to transitions between bit-strings in $\vec{p}_{\rm noisy}$.
An example is shown in Fig. (1b). Because we are grabbing only select elements
of the full $A$-matrix, the columns of $\tilde{A}$ must be renormalized such
that they once again sum to one. That is to say given any two bit-strings
$\texttt{row},\texttt{col}\in\vec{p}_{\rm noisy}$, the reduced matrix element
$\tilde{A}_{\texttt{row},\texttt{col}}$ is given by

$\tilde{A}_{\texttt{row},\texttt{col}}=\left\\{\begin{array}[]{cc}\cfrac{A_{\texttt{row},\texttt{col}}}{\sum_{\begin{subarray}{c}k\in\vec{p}_{\rm
noisy}\\\ d(k,\texttt{col})\leq D\end{subarray}}A_{k,\texttt{col}}}&\ \
\scriptstyle d(\texttt{row},\texttt{col})\leq D\\\ \\\ 0&\ \ \scriptstyle
d(\texttt{row},\texttt{col})>D\end{array}\right.$ (5)

where $D$ is the desired Hamming distance 111Keeping all elements is
equivalent to setting $D$ equal to the number of measured qubits, whereas
$D=0$ yields the identity matrix.. Performing a finite number of circuit
executions, even when measurement errors are weak, may result in completely
redistributing probability away from some small magnitude elements in
$\vec{p}_{\rm ideal}$ such that those elements are not present in
$\vec{p}_{\rm noisy}$; the solution vector will be missing these elements.
However as we will show, for typical numbers of circuit samples this effect is
minimal.

##  III Matrix-free solution

Although $\tilde{A}$ is much smaller than the original, when sampling circuits
with wide probability distributions many times, or executing on systems with
appreciable error rates, it is possible that $\tilde{A}$ itself may become too
costly to explicitly construct. Fortunately, being able to grab elements of
$A$ individually, Eq. (4), allows us to take advantage of matrix-free
iterative techniques [39]. The time to solution for iterative methods greatly
depends on the properties of $\tilde{A}$. $A$-matrices obtained from present
day cloud-accessible platforms nominally have strict diagonal dominance
$|A_{\tt row,\tt row}|>\sum_{\tt col\neq\tt row}|A_{\tt row,\tt
col}|~{}\forall~{}\tt row$, and are readily solved by simple iterative methods
such as Jacobi iteration [39]. However this condition does not hold in general
for systems with large error rates. Moreover, as the number of qubits grows so
does the number of possible error channels (i.e. number of states at low
Hamming distance) and it becomes harder to satisfy this stringent condition.
Thus general purpose methods such as generalized minimal residual (GMRES) [40]
or biconjugate gradient stabilized (BiCGSTAB) [41] methods must be used.
Importantly these Krylov subspace methods require computing only the product
$\tilde{A}\vec{p}_{\rm noisy}$, but not $\tilde{A}$ itself [42, 39]. However
having strict diagonal dominance, or close to it, suggests that we can
increase the rate of convergence by using a simple Jacobi preconditioner
$P^{-1}$ to solve

$P^{-1}\tilde{A}\vec{x}=P^{-1}\vec{p}_{\rm noisy},$ (6)

where $P^{-1}$ is a diagonal matrix with $P^{-1}_{i,i}=1/\tilde{A}_{i,i}$
[39]. In practice, Eq. (6) gives rapid convergence requiring only
$\mathcal{O}(1)$ iterations for an absolute tolerance value of $10^{-5}$ while
simultaneously dramatically reducing the memory requirements for mitigation.

##  IV Uncertainty estimates

Mitigating measurement errors does not come for free. Rather it results in an
increase in the uncertainty of repeated measurement outcomes that must be
compensated for by increasing the number of times the circuit is sampled. This
mitigation overhead $\mathcal{M}$ is determined by the one-norm of the inverse
of the A-matrix $\mathcal{M}=||A^{-1}||^{2}_{1}$ [11], and gives an upper
bound on the standard deviation of an observable
$\sigma_{O}\leq\sqrt{\mathcal{M}/s}$, where $s$ is the number of samples. Not
wanting to construct $\tilde{A}^{-1}$, here we use the iterative Hager-Higham
algorithm [43, 44] for estimating $||\tilde{A}^{-1}||_{1}$ using only linear
systems of equations involving $\tilde{A}$ and $\tilde{A}^{T}$. When using
direct factorization, the LU decomposition of $\tilde{A}$ can be cached, and
thus the overall run-time is $\sim 2x$ longer than mitigation alone. However,
for iterative methods, the overhead is between 4-10x longer depending on how
many steps the Hager-Higham routine requires. Although this method gives a
lower-bound on $\sqrt{\mathcal{M}}$, in practice it is often exact or nearly
so (see related discussion in Ref. [44]). Because our truncation method
selects only those rows and columns from $\vec{p}_{\rm noisy}$, the one-norm
of $||\tilde{A}^{-1}||_{1}$, and thus mitigation overhead, is dependent on the
circuit being executed and the noise properties of the device on which it is
run.

##  V Demonstrations

Our method [45] is implemented with NumPy [46], SciPy [47], and Cython [48],
and makes use of Qiskit [28] for calibration circuit construction and
execution. All timing data is taken on a quad core Intel i3-10100 system with
$32~{}\rm Gb$ of memory with NumPy and SciPy compiled against OpenBlas [49].

We begin by showcasing the veracity of our method by comparing expectation
values for the circuit shown in Fig. (2a) computed with the full-space
tensored method from Qiskit Ignis [28] along with our method called M3 222This
stands for: Matrix-free Measurement Mitigation. varying the number of samples
taken per circuit. Here the circuits are executed on the 27 qubit IBM Quantum
Kolkata system, mapping the virtual circuit qubits to physical qubits
$[1,4,7,10,12,13,14,11,8,5,3,2]$. The calibration data and raw input samples
are identical for both mitigation methods. In Fig. (2b) we see that, despite
using _at most_ $371$ bit-strings ($9\%$ of the full-dimensionality), M3
closely matches the full-dimensional $A^{(T)}$ results over the entire range
of samples, and agrees remarkably well with the Ignis results for the numbers
of executions typically employed in practice, $>1000$. An example
$\tilde{A}^{(T)}$ generated by M3 is presented in App. (A.3). The sub-space
reduction results in a run-time performance improvement as well, with Qiskit
Ignis mitigation taking $\sim 3~{}\rm s$ per circuit, where as M3 took at most
$\sim 7~{}\rm ms$. Additionally, we see that the uncertainty bound given by
$\mathcal{M}$ closely matches the experimental values, and verifies the use of
this technique in reporting faithful error bounds. The inset of Fig. (2b) also
shows that, for a tolerance of $\leq 10^{-5}$, $\tilde{A}^{(T)}$ is well-
approximated by keeping only terms out to $D=3$.

Figure 2: a) 12-qubit circuit with an exact distribution of 43 unique bit-
strings, and an expectation value of $\sim 0.446$. Measurements are omitted
for brevity. b) Expectation values before and after mitigating measurement
errors for the circuit in (a) executed 100 times on the IBM Quantum Kolkata
system while varying the number of samples per circuit. Error bars show one
standard deviation, while the shaded region gives error bounds determined by
the computed mitigation overhead and number of samples. The same calibration
data is used for both Ignis and M3 and was computed using $8192$ samples per
circuit. Inset figure shows the absolute error when truncating
$\tilde{A}^{(T)}$ to a given Hamming distance for the data collected at 8192
samples. Figure 3: a) Comparison of expectation values for GHZ states out to
42 qubits on the IBM Quantum Brooklyn system using M3, and Qiskit Ignis
tensored and bounded least-squares methods. M3 direct, iterative, and a $D=3$
Hamming approximation all yield the same values to a tolerance of $\sim
10^{-5}$ and the associated shaded region shows the error bound from
$\mathcal{M}$. Numbers alongside the M3 data show the number of bit-strings in
$\vec{p}_{\rm noisy}$ for each number of qubits. Inset shows the computed
mitigation overhead $\mathcal{M}$. All calibration data and GHZ circuit
execution was performed using $8192$ shots per-circuit. b) Timing information
(best of three runs) for the mitigation methods presented in (a). Annotated
numerical values indicate the number of iterations needed for a tolerance
$\leq 10^{-5}$. The full and $D=3$ iterative solutions used the same number of
iterations in all cases.

We now demonstrate the scalability of our M3 method by mitigating GHZ states
out to 42 qubits on the 65 qubit IBM Quantum Brooklyn system. Details of this
experiment are in App. (C). In Fig. (3) we compare M3 along with the Qiskit
Ignis tensored and bounded least-squares methods 333We have modified the
Qiskit least-squares method to use $\vec{p}_{\rm noisy}$ as the starting
vector as opposed to a random vector. This gives a 3x or more improvement in
runtime.. Only M3 allows for mitigating errors beyond 14 qubits due to
algorithmic breakdown or extreme run-times, for the tensored and least-squares
methods respectively, and shows the importance of performing measurement
mitigation for large-scale experiments. The overall expectation values drop as
the circuit depth increases where gate errors and decoherence, effects
measurement mitigation cannot resolve, start to dominate.

The mitigation overhead, inset of Fig. (3a), shows exponential scaling at
small numbers of qubits after which the overhead begins to plateau. The
exponential scaling arises as the diagonal elements of $\tilde{A}$, formed
from the product of $N$ probabilities, Eq. (4), are effectively being inverted
when computing $\tilde{A}^{-1}$, with additional contributions coming from
elements close in Hamming distance. Provided that the bit-strings in
$\vec{p}_{\rm noisy}$ sample sufficient portions of these short Hamming
distance elements, the re-normalization used in obtaining $\tilde{A}$ is
small, and one recovers the exponential scaling shown for the full-dimensional
A-matrix [11]. However if this is not the case then re-normalization increases
the magnitude of the elements in $\tilde{A}$ (in particular the diagonal
elements) suppressing the exponential growth in $\mathcal{M}$.

Figure (3b) details the timing across the different mitigation methods. We see
that M3 greatly improves the computed expectation values while taking at most
$1.2$ seconds to compute at $42$ qubits. When computing the mitigation
overhead, the total time increased to $2.4$ and $4.5$ seconds for the direct
and iterative solutions at $42$ qubits, respectively (not shown); an
extraordinary improvement upon the exponential run times observed for the
Qiskit methods. As with the example in Fig. (2), a $D=3$ Hamming approximation
well captures the full mitigation process to the desired tolerance, and
performs best at large numbers of bit-strings where the overhead from
computing the Hamming distance between elements starts to become less than the
cost of additional floating-point multiplications. For the M3 results, the
rate at which the run times increase begins to slow down at larger numbers of
qubits following a similar slow down in the number of unique bit-strings
observed, Fig. (3a), when additional qubits are added.

Finally, we note that at $42$ qubits, storing a full $2^{N}$-vector of single-
precision floating point values for $\vec{p}_{\rm noisy}$ requires $16~{}\rm
TiB$ of memory; well beyond the limits of our computer on which the mitigation
is implemented, but is amenable to storage on a supercomputer. Juxtapose that
with storing a sparse representation of $A$ using, for example, compressed
sparse column (CSC) format out to $D=3$. This requires $\sim 580~{}\rm PiB$ of
memory, which is $120$x more than that available in the Fugaku supercomputer
[52] [see App. (C.2) for details]. In contrast, the M3 iterative method uses
$\sim 1~{}\rm MiB$ of storage, highlighting the benefit of the techniques
presented here for mitigating measurements at scales amenable to
demonstrations of quantum advantage.

###### Acknowledgements.

We thank Doug McClure, David McKay, and Matthew Treinish for helpful
discussions. Figures (2-7) are produced using Matplotlib [53].

##  Appendix A Example A-matrices

###  A.1 Complete A-matrix

An example complete A-matrix computed by running $2^{N}$ circuits, one for
each computational basis state, on the IBM Quantum Kolkata system is given in
Fig. (4). Because of finite sampling, the matrix is nominally sparse, and only
those elements close in Hamming distance have appreciable transition
probabilities. The matrix has strict diagonal dominance and thus is guaranteed
to be invertible.

Figure 4: Full A-matrix for qubits $0\rightarrow 5$ on the IBM Quantum Kolkata
system. A circuit for each of the 32 bit-strings was executed 8192 times to
fill in the columns.

###  A.2 Tensored A-matrix

The A-matrix corresponding to tensored measurement errors, $A^{(T)}$ is
constructed by taking the tensor product of single-qubit calibration matrices
given by Eq. (3) in the main text. Unlike the complete A-matrix, $A^{(T)}$
contains only non-zero elements unless one or more qubits has no reported
measurement error for $P^{(k)}_{0,1}$ and/or $P^{(k)}_{1,0}$; $A^{(T)}$
expects every element of $\vec{p}_{\rm noisy}$ to have a non-zero entry. Like
Fig. (4), the matrix in Fig. (5) is strictly diagonally dominant, and
indicates that transitions between elements close in Hamming distance are more
likely. Data for Fig. (5) was taken immediately after that shown in Fig. (4).

Figure 5: Tensored A-matrix $A^{(T)}$ for qubits $0\rightarrow 5$ on IBM
Quantum Kolkata.

###  A.3 Example 12 qubit truncated A-matrix

Figure (6) shows one of the $100$ truncated $\tilde{A}^{(T)}$ used in the M3
mitigation performed in Fig. (2) at $8192$ counts. This matrix is also
strictly diagonally dominant. For each circuit execution, the number of
elements in $\tilde{A}^{(T)}$ may vary, as can their associated amplitudes.

Figure 6: Sample $\tilde{A}$ from the M3 results in Fig. (2b) taken at 8192
samples. The matrix contains $341$ bit-strings out of a possible $4096$. Rows
and columns are sorted in lexicographical order in terms of their bit-string
representation.

##  Appendix B Correlated errors

Although we have focused on tensored mitigation techniques, our method is
equally capable of handling correlated errors provided that matrix elements
$\tilde{A}_{\tt row,\tt col}$ can be obtained by their bit-string values as
done in Eq. (4). IBM Quantum systems operating normally are dominated by
uncorrelated measurement error, and are thus well mitigated by the tensored
A-matrix methods presented in the main text. It is possible however to
purposely induce correlated errors into the readout process and explore
correlated mitigation strategies.

To take into account pairwise correlated errors, we can modify Eq. (4),
grabbing elements using

$A_{\tt row,\tt
col}^{(C)}=\frac{1}{\binom{N}{2}}\sum_{k=1}^{N-1}\sum_{l=0}^{k-1}C_{kl}[a^{kl},b^{kl}]\prod^{N-1}_{\begin{subarray}{c}m=0\\\
m\neq k,l\end{subarray}}S_{m}[q_{m},q^{\prime}_{m}],$ (7)

where $\texttt{row}=q_{N-1}q_{N-2}..q_{0}$ and
$\texttt{col}=q^{\prime}_{N-1}q^{\prime}_{N-2}..q^{\prime}_{0}$. Here, the
$S_{m}$ are defined as in Eq. (3), and $C_{kl}$ are a $4x4$ stochastic matrix
(local noise matrix) between qubits $k$ and $l$ where $a^{kl}=2q_{k}+q_{l}$,
$b^{kl}=2q^{\prime}_{k}+q^{\prime}_{l}$, respectively. The elements of
$C_{kl}$ are obtained following Section $V$ of Ref. [11].

IBM Quantum Kolkata, a Falcon 5.11 series system, has readout output
multiplexing ratios ranging from $3:1$ to $5:1$, with readout frequencies in a
shared output typically separated by $50$-$60~{}\rm MHz$. This separation is
much larger than the average cavity linewidth $\kappa$ ($5.6~{}\rm MHz$) and
dispersive shift $\chi$ ($1.6~{}\rm MHz$), which when combined enable short
$330~{}\rm ns$ readout pulses for all qubits. Using default readout pulse
amplitudes, calibrated to optimize fidelity while maintaining QND readout,
Fig. (7a) shows expectation values produced by executing an 8 qubit GHZ
circuit 100 times using qubits $[8,5,3,2,1,4,7,6]$ that span two multiplex
readout groupings. Importantly we perform mitigation using the complete
$A$-matrix, $\tilde{A}^{(C)}$, as well as $\tilde{A}^{(T)}$, with results in
agreement with uncorrelated errors largely dominating the readout process.

Intentionally increasing readout pulse amplitudes by approximately 2x from the
optimized values results in correlated readout errors. With these larger
readout amplitudes, while there is no appreciable change in average readout
fidelity ($<0.2\%$) compared to the default setting, there is a substantial
uptick in non-QNDness. In Fig. (7b) we re-run our 8 qubit GHZ experiment from
Fig. (7a) under these new conditions. While the correlated M3 method using Eq.
(7) well captures the the correlated readout errors, as evident by agreement
with the complete A-matrix, the tensored mitigator strongly over-corrects due
to the misalignment of probabilities in $\tilde{A}^{(T)}$ with those actually
present in the system. Although Eq. (7) works for both uncorrelated and
pairwise correlated errors, each matrix element requires
$\mathcal{O}\left(N^{2}\right)$ floating-point evaluations, as opposed to $N$
in Eq. (4).

Figure 7: a) Histogram of expectation values collected from executing an 8
qubit GHZ circuit 100 times on the IBM Quantum F5.11 Kolkata system under
normal operating conditions and using qubits $[8,5,3,2,1,4,7,6]$ and $8192$
samples per circuit. Raw results are mitigated with the complete ($A$),
correlated M3 ($\tilde{A}^{(C)}$), and tensored M3 ($\tilde{A}^{(T)}$). b)
Repeated experimental results in the presence of correlated readout errors
generated by degraded readout.

##  Appendix C 42 qubit GHZ demonstration

###  C.1 Experimental details

Experiments are run on the 65 qubit IBM Quantum Brooklyn system. GHZ states
were prepared starting with a Hadamard gate on qubit $11$ and entangling
additional qubits as shown in Fig. (8). After entangling the first 6 qubits,
this pattern allows for increasing the GHZ state by 4 qubits per layer. The
average assignment and CNOT error rates across the qubits used is $2.15\%$ and
$1.01\%$, respectively.

Figure 8: Qubits used in generating GHZ states on the IBM Quantum Brooklyn
system. Gates entangled earlier are color coded darker. The initial two qubit
GHZ state is between qubits 11 and 4, while the final 42 qubit GHZ state is
formed by adding qubits 9, 37, 40, and 45 to the previous iteration.

###  C.2 Memory requirements for storing full 42 qubit A-matrix to $D=3$

Including elements up to a Hamming distance of three requires

$\binom{42}{0}+\binom{42}{1}+\binom{42}{2}+\binom{42}{3}=12384$ (8)

elements in each of the $2^{42}$ columns. Storing these values using single-
precision floating-point numbers, 4-bytes per entry, requires $193.5~{}\rm
PiB$ of memory. In addition, we must also specify the row and column indices
for these values. In compressed sparse column (CSC) format we need $2^{42}+1$
elements, the difference of which specify the number of non-zero elements in
each of the columns. Lastly, we also need the row index for each non-zero
matrix element. At 42 qubits, the size of these indices cannot be stored using
32-bit integers, and we must use 64-bits per entry. Storing these values
requires an additional $387~{}\rm PiB$ of memory. The total memory required is
therefore $580.5~{}\rm PiB$ and is $\sim 120$x larger than the $4.85~{}\rm
PiB$ of memory on the Fugaku supercomputer [52].

## References

  * Bravyi  _et al._ [2018] S. Bravyi, D. Gosset, and R. König, Quantum advantage with shallow circuits, Science 362, 308 (2018).
  * Bravyi  _et al._ [2020] S. Bravyi, D. Gosset, R. König, and M. Tomamichel, Quantum advantage with noisy shallow circuits, Nat. Physics 16, 1040 (2020).
  * Maslov  _et al._ [2021] D. Maslov, J.-S. Kim, S. Bravyi, T. J. Yoder, and S. Sheldon, Quantum advantage for computations with limited space, Nat. Physics 10.1038/s41567-021-01271-7 (2021).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error Mitigation for Short-Depth Quantum Circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Phys. Rev. X 8, 031027 (2018).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Extending the computational reach of a noisy superconducting quantum processor, Nature 567, 491 (2019).
  * McArdle  _et al._ [2019] S. McArdle, X. Yuan, and S. C. Benjamin, Error-Mitigated Digital Quantum Simulation, Phys. Rev. Lett. 122, 180501 (2019).
  * Giurgica-Tiron  _et al._ [2020] 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)_ (2020) p. 306.
  * Sun  _et al._ [2021] J. Sun, X. Yuan, T. Tsunoda, V. Vedral, and S. C. Benjamin, Mitigating Realistic Noise in Practical Noisy Intermediate-Scale Quantum Devices, Phys. Rev. Applied 15, 034026 (2021).
  * Kim  _et al._ [2021] Y. Kim, C. J. Wood, T. J. Yoder, S. T. Merkel, J. M. Gambetta, K. Temme, and A. Kandala, Scalable error mitigation for noisy quantum circuits produces competitive expectation values, arXiv:2108.09197 (2021).
  * Bravyi  _et al._ [2021] S. Bravyi, S. Sheldon, A. Kandala, D. C. McKay, and J. M. Gambetta, Mitigating measurement errors in multi-qubit experiments, Phys. Rev. A 103, 042605 (2021).
  * Geller [2020] M. R. Geller, Rigorous measurement error correction, Quantum Science and Technology 5, 03LT01 (2020).
  * Geller and Sun [2020] M. R. Geller and M. Sun, Efficient correction of multiqubit measurement errors, arXiv:2001.09980 (2020).
  * Hamilton  _et al._ [2020] E. Hamilton, Kathleen, T. Kharazi, T. Morris, A. J. McCaskey, R. S. Bennink, and C. Pooser, Raphael, Scalable quantum processor noise characterization, arXiv:2006.01805 https://arxiv.org/abs/2006.01805 (2020).
  * van den Berg  _et al._ [2020] E. van den Berg, K. Minev, Zlatko, and K. Temme, Model-free readout-error mitigation for quantum expectation values, arXiv:2012.09738 (2020).
  * Maciejewski  _et al._ [2020] F. B. Maciejewski, Z. Zimborás, and M. Oszmaniec, Mitigation of readout noise in near-term quantum devices by classical post-processing based on detector tomography, Quantum 4, 257 (2020).
  * Nachman  _et al._ [2020] B. Nachman, M. Urbanek, W. A. de Jong, and C. W. Bauer, Unfolding quantum computer readout noise, npj Quantum Inf 6, 84 (2020).
  * Hicks  _et al._ [2021] R. Hicks, C. W. Bauer, and B. Nachman, Readout rebalancing for near-term quantum computers, Phys. Rev. A 103, 022407 (2021).
  * Wang  _et al._ [2021] K. Wang, Y.-A. Chen, and X. Wang, Measurement Error Mitigation via Truncated Neumann Series, arXiv:2103.13856 (2021).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Gong  _et al._ [2019] M. Gong, M.-C. Chen, Y. Zheng, S. wang, C. Zha, H. Deng, Z. Yan, H. Rong, Y. Wu, S. Li, F. Chen, Y. Zhao, F. Liang, J. Lin, Y. Xu, C. Guo, L. Sun, A. D. Castellano, H. Wang, C. Peng, C.-Y. Lu, X. Zhu, and J.-W. Pan, Genuine 12-Qubit Entanglement on a Superconducting Quantum Processor, Phys. Rev. Lett. 122, 110501 (2019).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Wei  _et al._ [2020] K. X. Wei, I. Lauer, S. Srinivasan, N. Sundaresan, D. T. McClure, D. Toyli, D. C. McKay, J. M. Gambetta, and S. Sheldon, Verifying multipartite entangled Greenberger-Horne-Zeilinger states via multiple quantum coherences, Phys. Rev. A 101, 032343 (2020).
  * Mooney  _et al._ [2021a] G. J. Mooney, G. A. L. White, C. D. Hill, and L. C. L. Hollenberg, Generation and verification of 27-qubit Greenberger-Horne-Zeilinger states in a superconducting quantum computer, arXiv:2101.08946 (2021a).
  * Mooney  _et al._ [2021b] G. J. Mooney, G. A. L. White, C. D. Hill, and L. C. L. Hollenberg, Whole-device entanglement in a 65-qubit superconducting quantum computer, arXiv:2102.11521 (2021b).
  * Satzinger  _et al._ [2021] K. J. Satzinger, Y. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, C. Gidney, I. Aleiner, F. Arute, K. Arya, J. Atalaya, R. Babbush, J. C. Bardin, R. Barends, J. Basso, A. Bengtsson, A. Bilmes, M. Broughton, B. B. Buckley, D. A. Buell, B. Bukett, N. Bushnell, B. Chairo, R. Collins, W. Courtney, S. Bemura, A. R. Derk, D. Eppens, C. Erickson, E. Farhi, L. Foaro, A. G. Fowler, B. Foxen, M. Giustina, A. Greene, J. A. Gross, M. P. Harrington, S. D. Harrington, J. Hilton, S. Hong, T. Huan g, W. J. Huggins, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, T. Khattar, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, A. Locharla, E. Lucero, O. Martin, J. R. McClean, M. McEwen, K. C. Miao, M. . Mohseni, S. Montazeri, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, T. E. O’Brien, B. Villalonga, T. C. White, Z. Yao, P. Yeh, J. Yoo, A. Zalcman, H. Neven, S. Boixo, A. Megrant, Y. Chen, J. Kelly, V. Smelyanskiy, A. Kitaev, M. Knap, F. Pollmann, and P. Roushan, Realizing topologically ordered states on a quantum processor, arXiv:2104.01180 (2021).
  * Glick  _et al._ [2021] R. Glick, Jennifer, P. Gujarati, Tanvi, A. D. Córcoles, Y. Kim, A. Kandala, J. M. Gambetta, and K. Temme, Covariant quantum kernels for data with group structure, arXiv:2105.03406 (2021).
  * [28] Qiskit 0.26, https://qiskit.org.
  * Pashayan  _et al._ [2015] H. Pashayan, J. J. Wallman, and S. D. Bartlett, Estimating Outcome Probabilities of Quantum Circuits Using Quasiprobabilities, Phys. Rev. Lett. 115, 070501 (2015).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, and A. Aspuru-Guzik, A variational eigenvalue solver on a photonic quantum processor., Nat. Commun. 5, 4213 (2014).
  * ali [2020] quantumcomputer.ac.cn (Accessed August 27, 2020).
  * goo [2021] quantumai.google/hardware (Accessed June 01, 2021).
  * iqx [2021] quantum-computing.ibm.com/systems (Accessed August 23, 2021).
  * aws [2021] aws.amazon.com/braket (Accessed June 01, 2021).
  * ins [2020] www.quantum-inspire.com (Accessed August 27, 2020).
  * Nation  _et al._ [2015] P. D. Nation, J. R. Johansson, M. P. Blencowe, and A. J. Rimberg, Iterative solutions to the steady-state density matrix for optomechanical systems, Phys. Rev. E 91, 013307 (2015).
  * Smolin  _et al._ [2012] J. A. Smolin, J. M. Gambetta, and G. Smith, Efficient Method for Computing the Maximum-Likelihood Quantum State from Measurements with Additive Gaussian Noise, Phys. Rev. Lett. 108, 070502 (2012).
  * Note [1] Keeping all elements is equivalent to setting $D$ equal to the number of measured qubits, whereas $D=0$ yields the identity matrix.
  * Saad [2003] Y. Saad,  _Iterative Methods for Sparse Linear Systems_ , edited by 2nd (Society for Industrial and Applied Mathematics, 2003).
  * Saad and Schultz [1986] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. and Stat. Comput. 7, 856 (1986).
  * van der Vorst [1992] H. A. van der Vorst, Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems, SIAM J. Sci. and Stat. Comput. 13, 631 (1992).
  * Barrett  _et al._ [1994] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. Van er Vorst,  _Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods_ , 2nd ed. (SIAM, 1994).
  * Hager [1984] W. W. Hager, Condition estimates, SIAM J. Sci. Star. Comput. 5, 311 (1984).
  * Higham [1988] N. J. Higham, FORTRAN Codes for Estimating the One-Norm of a Real or Complex Matrix with Applications to Condition Estimation, ACM Trans. Math. Softw. 14, 381 (1988).
  * [45] https://github.com/qiskit-partners/mthree.
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, _et al._ , SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261 (2020).
  * Behnel  _et al._ [2011] S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D. S. Seljebotn, and K. Smith, Cython: The Best of Both Worlds, Computing in Science & Engineering 13, 31 (2011).
  * [49] http://www.openblas.net.
  * Note [2] This stands for: Matrix-free Measurement Mitigation.
  * Note [3] We have modified the Qiskit least-squares method to use $\vec{p}_{\rm noisy}$ as the starting vector as opposed to a random vector. This gives a 3x or more improvement in runtime.
  * fug [2021] https://www.r-ccs.riken.jp/en/fugaku/about/ (Accessed July 05, 2021).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing in Science & Engineering 9, 90 (2007).


4}+{\pi^{2}\over 6}+{\cal O}(\epsilon)\right]\;.$ (224)

In the theory of dilute Bose gases, the following integrals appear braaten

$\displaystyle I_{m,n}(M^{2})$ $\displaystyle=$
$\displaystyle\int_{p}{p^{2m}\over p^{n}\left(p^{2}+M^{2}\right)^{{n\over
2}}}\;.$ (225)

They satisfy the recursion relation

$\displaystyle{dI_{m,n}\over dM^{2}}$ $\displaystyle=$ $\displaystyle-{1\over
2}nI_{m+1,n+2}(M^{2})\;,$ (226)

which follows directly from the definition Eq. (225). Evaluating the integrals
in dimensional regularization, we find

$\displaystyle I_{m,n}(M^{2})$ $\displaystyle=$ $\displaystyle
e^{\gamma_{E}\epsilon}{M^{3+2m-2n}\over(4\pi)^{{3\over 2}}}\left({\Lambda\over
M}\right)^{2\epsilon}{\Gamma({d-n\over 2}+m)\Gamma(n-m-{d\over
2})\over\Gamma({n\over 2})\Gamma({d\over 2})}\;.$ (227)

We specifically need

$\displaystyle I_{0,-1}(M^{2})$ $\displaystyle=$ $\displaystyle{16\over
15}{M^{5}\over(4\pi)^{2}}\left[1+{\cal O}(\epsilon)\right]\;,$ (228)
$\displaystyle I_{1,1}(M^{2})$ $\displaystyle=$ $\displaystyle{16M^{3}\over
3(4\pi)^{2}}\left[1+{\cal O}(\epsilon)\right]\;.$ (229)

The integrals are finite in the limit $d\rightarrow 3$ reflecting that the
ultraviolet divergences are powerlike.

## References

  * (1) K. Rajagopal and F. Wilczek, At the frontier of particle physics, Vol. 3 (World Scientific, Singapore, p 2061) (2001).
  * (2) M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, Rev. Mod. Phys. 80, 1455 (2008).
  * (3) K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 014001 (2011).
  * (4) S. Weinberg, Physica A 96, 327 (1979).
  * (5) J. Gasser and H. Leutwyler, Ann. Phys. 158, (142) (1984).
  * (6) J. Gasser and H. Leutwyler, Nucl. Phys. B 250, 465 (1985).
  * (7) J. Bijnens, G. Colangelo and G. Ecker, Ann. Phys. 280, 100 (2000).
  * (8) D. T. Son and M. A. Stephanov, Phys. Rev. Lett. 86, 592 (2001); Phys. Atom. Nucl. 64, 834 (2001).
  * (9) K. Splittorff, D. T. Son, M. A. Stephanov, Phys. Rev. D 64, 016003 (2001).
  * (10) J. B. Kogut and D. Toublan,
  * (11) J. B. Kogut and D. K. Sinclair, Phys. Rev. D 66, 014508 (2002).
  * (12) J. B. Kogut and D. K. Sinclair, Phys. Rev D 66 034505 (2002).
  * (13) J. B. Kogut and D. K. Sinclair, Phys. Rev D 70 094501 (2004).
  * (14) D. K. Sinclair and J. B. Kogut, PosLat 2006 147 (2006).
  * (15) D. K. Sinclair, J. B. Kogut, PoSLAT 2007 225 (2007).
  * (16) B. B. Brandt and G. Endrődi, PoS LATTICE 2016, 039 (2016).
  * (17) B. B. Brandt, G. Endrődi, and S. Schmalzbauer, EPJ Web Conf. 175, 07020 (2018).
  * (18) B. B. Brandt, G. Endrődi, and S. Schmalzbauer, Phys. Rev. D 97, 054514 (2018).
  * (19) B. B. Brandt and G. Endrődi Phys. Rev. D 99, 014518 (2019)
  * (20) B.B. Brandt, F. Cuteri, G. Endrődi, and S. Schmalzbauer, Particles 3, 80 (2020), e-Print: 1912.07451 [hep-lat].
  * (21) B. B. Brandt, Francesca Cuteri, and G. Endrődi, PoS LATTICE2021, 232 (2022), e-Print: 2112.11113 [hep-lat]
  * (22) B. B. Brandt, Francesca Cuteri, and G. Endrődi, JHEP 07, 055 (2023).
  * (23) T. D. Cohen, Phys. Rev. Lett. 91, 222001 (2003).
  * (24) S. Carignano, A. Mammarella, and M. Mannarelli, Phys. Rev. D 93, 051503 (2016).
  * (25) S. Carignano, L. Lepori, A. Mammarella, M. Mannarelli, and G. Pagliaroli Eur. Phys. J. A 53, 35 (2017).
  * (26) T. Xia, L. He, and P. Zhuang, Phys. Rev. D 88, 056013 (2013);
  * (27) S. S. Avancini, A. Bandyopadhyay, D. C. Duarte, and R. L. S. Farias Phys. Rev. D 100, 116002 (2019).
  * (28) A. Ayala, A. Bandyopadhyay, R. L. S. Farias, L. A. Hernandez, and J. L. Hernandez, Phys. Rev. D 107, 074027 (2023).
  * (29) R. Chiba and T. Kojo, arXiv:2304.13920 [hep-ph].
  * (30) A. Ayala, B. S. Lopes, R. L. S. Faria, L. C. Parra, e-Print: 2310.13130 [hep-ph].
  * (31) M. S. Grønli and T. Brauner, Eur. Phys. J. C 82, 354 (2022).
  * (32) M. Mannarelli, Particles 2, 411 (2019).
  * (33) G. Ecker, J. Gasser, A. Pich, and E. de Rafael, Nucl. Phys. B 321, 311 (1989).
  * (34) R. Urech, Nucl. Phys. B 433, 234 (1995).
  * (35) R. Dashen, Phys. Rev. 183, 1245 (1969).
  * (36) U.-G. Meißner, G. Müller and S. Steininger, Phys. Lett. B 406, 154 (1997).
  * (37) J. Gasser, A. Rusetsky, and I. Scimemi, Eur. Phys. J. C 32, 97 (2003).
  * (38) M. Knecht and R. Urech, Nucl. Phys. B 519, 329 (1998).
  * (39) H. Leutwyler, Annals Phys. 235, 165 (1994).
  * (40) A. Pich, Rept. Prog. Phys. 58, 563 (1995).
  * (41) B. Gripaios, e-Print:1506.05039 [hep-ph]
  * (42) R. Penco, e-Print: 2006.16285 [hep-th].
  * (43) J. I. Kapusta, Phys. Rev. D 24, 426 (1981).
  * (44) H. E. Haber and H. A. Weldon Phys. Rev. D 25, 502 (1982).
  * (45) S. Scherer, e-Print: 0210398 [hep-ph].
  * (46) J. Bijnens, G. Colangelo, and G. Ecker, Ann. Phys. 280, 100 (2000).
  * (47) J. Bijnens, G. Colangelo, and G. Ecker, JHEP 02 020 (1999).
  * (48) M. K. Johnsrud, Master thesis, NTNU 2022. Thesis and code available at https://github.com/martkjoh/master.
  * (49) P. Adhikari, J. O. Andersen, and P. Kneshcke, Eur. Phys. J. C 79, 874 (2019).
  * (50) A. Nicolis and F. Piazza, JHEP 06, 025 (2012).
  * (51) A. Nicolis, R. Penco, F, Piazza, and R. Rattazzi, JHEP 06 155 (2015).
  * (52) P. A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * (53) J. Gasser and H. Leutwyler, Phys. Rep. 87, 77 (1982).
  * (54) A. Joyce, A. Nicolis, A. Podo, and L. Santoni, JHEP 09 066 (2022).
  * (55) I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products. Eighth edition (2014).
  * (56) J. Bijnens and G. Ecker, Ann. Rev. of Nuclear and Particle Science, 64, 149 (2014).
  * (57) P. Adhikari, J. O. Andersen, and M. A. Mojahed, Eur. Phys. J. C 81, 449 (2021).
  * (58) D.T. Son, e-Print: 0204199 [hep-ph].
  * (59) C. Manuel, A. Dobado, and F. J. Llanes-Estrada, JHEP 09, 076 (2005).
  * (60) C. Manuel and F. J. Llanes-Estrada, JCAP 08 001 (2007).
  * (61) M. A. Escobedo and C. Manuel, Phys. Rev. A 82, 023614 (2010).
  * (62) S. T. Beliaev, Sov. J. Phys. 7, 289 (1958).
  * (63) J. O. Andersen, Q. Yu, and H. Zhou, e-Print: 2306.14472 [hep-ph].
  * (64) L. He, Phys. Rev. D 82, 096003 (2010).
  * (65) N. N. Bogoliubov, J. Phys. (USSR) 11, 23 (1947).
  * (66) E. Braaten and A. Nieto Phys. Rev. B 56, 14745 (1997).
  * (67) T. D. Lee and C. N. Yang, Phys. Rev. 105, 1119 (1957).
  * (68) T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, 1135 (1957).
  * (69) T. T. Wu, Phys. Rev. 115, 1390 (1959).
  * (70) N. M. Hugenholz and D. Pines, Phys. Rev. 116, 489 (1959).
  * (71) K. Sawada, Phys. Rev. 116, 1344 (1959).
  * (72) E. Braaten and A. Nieto, Eur. Phys. J. B 11, 143 (1999).
  * (73) E. Braaten, H.-W. Hammer, and S. Hermans, Phys. Rev. A 63, 063609 (2001).
  * (74) A. Nicolis, A. Podo, and L. Santoni, JHEP 09, 200 (2023)
  * (75) P. Adhikari, J. O. Andersen, and M. A. Mojahed, Eur.Phys. J. C 81, 173 (2021).