Standard Model Physics and the Digital Quantum Revolution: Thoughts about the Interface

A necessary and sufficient condition for identifying quantum entanglement can be expressed in terms of the separability of the density matrix describing the union of two or more Hilbert spaces $\mathcal{H}=\cup_{i}\mathcal{H}_{i}$. For the specific case of $\mathcal{H}=\mathcal{H}_{A}\cup\mathcal{H}_{B}$, the two Hilbert spaces may be e.g., spatial regions of a field, separate momentum scales, valence and sea quantum number contributions, or a spatial region of a nucleon and its complement with a size set by the probing momentum transfer. If the density matrix can be written as a classical mixture of separable density matrices, $\rho_{A\cup B}=\sum_{i}p_{i}\ \rho_{A,i}\otimes\rho_{B,i}$ with $p_{i}\geq 0$, then the system is separable and not entangled. If the density matrix cannot be written in this way, the system is entangled. The existence of entangled states, flouting attempts to describe their correlations classically, is culpable for the exponential growth of Hilbert spaces (and thus classical computational resources) required to describe QMB systems. To date, every theoretical attempt to evade this exponential state space required by non-classical correlations has been ruled out by experiment, the most famous example being that of nature’s consistent response, through the violation of Bell’s inequalities, that a local hidden variable description is inconsistent with observation.

While experiment indicates that a complete evasion of the exponential growth of Hilbert space is not likely to provide a successful description of nature, it is possible to systematically delay the onset of this burdensome scaling. The language of one concrete and powerful example Bañuls et al. (2018) contains at its heart the pure-state Schmidt decomposition,

the expansion of a pure quantum state in terms of the separate eigenvectors of its $A$- and $B$-space reduced density matrices. For pure states, the two density matrix spectra, $\vec{\lambda}$, will have equal non-zero eigenvalues of number limited by the dimensionality of the smaller Hilbert space, $\chi_{m}$. Only for $\chi_{m}$ = 1, in which local pure states remain, is the state separable across the chosen bipartition. Entanglement-motivated approximations to the state in Eq. (1) can be achieved by introducing a truncation $\chi\leq\chi_{m}$, with a viable and affable entanglement measure identified as $\chi$’s logarithm. Despite the opacity of entanglement structure when an $n$-qubit state is written in the computational basis, $|\Psi\rangle=\sum\psi_{\vec{b}}\ |b_{1}\rangle\otimes\cdots\otimes|b_{n}\rangle$ with $b_{i}\in\{0,1\}$, an entanglement-driven hierarchical structure of the wavefunction can be exposed through a series of sequential Schmidt decompositions  Vidal (2003, 2004),

Here, the $\Gamma$’s are tensors of maximal dimension $\chi$ characterizing the single-qubit Hilbert space and the $\lambda$’s are locally-contracted vectors with a dimensionality of $\chi$ at the bipartition. These $\lambda$-vectors, serving as entanglement conduits between local qubit Hilbert spaces, may be truncated, leading to a family of systematically improvable approximations to $|\psi\rangle$ with reduced entanglement. Rather than the $2^{n}$ complex numbers capturing the exact wavefunction, the truncated description requires a memory that scales $\sim\chi^{2}n$, linearly in the number of qubits but exponentially in the entanglement measure mentioned above. Efficiently representing low-entanglement quantum states by truncating their local Schmidt decompositions underpins the robust technology of tensor networks Bañuls et al. (2018), continuing to provide leading capabilities in low spatial dimensions and for short time evolutions where local bipartite entanglement is naturally constrained.

Though the clear separability criterion above seems at first simple and innocuous, determining whether a particular density matrix can be transformed into a separable mixture is, in general, a difficult and open area of research. In particular, there is no known necessary and sufficient criterion for separability that is also efficiently computable as the dimensionality and/or number of Hilbert spaces becomes large. Unfortunately, large in the previous sentence happens to be around three qubits, a pair of three-state qutrits, or the presence of a four-state qudit. It may be no surprise that describing the complex non-local correlations that quantum systems are capable of expressing is not easily achieved through a single observable or criterion. Rather, analogous to the Bayesian approach to data analysis through classical probabilities as a “dialogue with the data” Sivia , the exploration and characterization of entanglement in quantum mechanical systems often requires a collection of entanglement criteria, each one perhaps being necessary or sufficient, whose combination of fragmented information may be synthesized into more complete logical conclusions.

For example, a quantum generalization of the Shannon entropy leads to the entanglement entropy, $\mathcal{S}=-Tr[\rho_{A}\log\rho_{A}]$, used to explore the entanglement of bipartite states that are globally pure. The entanglement entropy is calculated in terms of the reduced density matrix, $\rho_{A}=\text{Tr}_{B}\rho_{A\cup B}$, and is symmetric with respect to the system $A,B$ that is chosen to be traced over. Because entanglement entropy heralds the presence of entanglement by alerting us to mixedness in reduced density matrices (pure states exhibit $\mathcal{S}=0$), the interpretation of this entropy becomes ineffectual when the global state is already mixed. For the language of mixed states, one may turn to conditional entropy and mutual information that combine otherwise-ambiguous marginal distributions into more informative relative entropies. The important caveats continue as such mutual information criteria are sensitive to classical correlations in addition to quantum correlations; it is possible to find a non-zero mutual information for a classically mixed unentangled state, and thus non-zero mutual information is necessary though insufficient for detecting entanglement. To address this, one may look yet further to the family of positive transformations, those that retain positive eigenvalues of the density matrix. By acting locally (subsystem) with such a transformation, any emergence of a negative eigenvalue must herald inseparability. A common positive transformation to consider, due to its relatively amiable computability, is the transpose Horodecki et al. (1996); Simon (2000); Vidal and Werner (2002); Plenio (2005); Horodecki (1997); Horodecki et al. (1998). Physically, the partial transpose implements a local momentum inversion that will produce a valid, positive density matrix if the two regions distinguished by the transpose are unentangled. Though the presence of negative eigenvalues, non-zero negativity, is unperturbed by classical correlations, negativity is sufficient for detecting inseparability, but it is not necessary Horodecki (1997); Horodecki et al. (1998); Horodecki and Lewenstein (2000); Werner, R. F. and Wolf, M. M. (2001). Aptly named, the multitude of entanglement or separability witnesses Terhal (2000), will individually fail to identify the properties they are designed to observe in a subset of Hilbert space; their strength, however, often lies in collaboration among multiple witnesses and deductive reasoning. Following Murphy’s law, and in the spirit of EFTs, it is usually a safe default to assume: if states of particular entanglement properties can exist (e.g., positive under partial transposition but inseparable), they do.

In addition to entanglement witnesses that conclusively herald but inconclusively exclude entanglement or separability, is a category of measures regarding entanglement more fundamentally as a resource for quantum information protocols. Examples of such measures are the entanglement of formation Wootters (1998); Hill and Wootters (1997) or the distillable entanglement Bennett et al. (1996a, b), quantifying the asymptotic ratio of maximally entangled resource states necessary to produce an ensemble of the state or the ratio capable of being distilled back from a given ensemble through local operations, respectively. While operational measures provide welcomed clear physical interpretations, they are commonly computationally prohibitive. It is thus common to combine heralding measures with operational measures, e.g., the logarithmic negativity is an upper bound to the distillable entanglement Vidal and Werner (2002); Audenaert et al. (2003). Consistent with the above creation of entanglement measures in a collaborative array connected to e.g., operational communications protocols, it is likely that developing the dialogue for exploring the entanglement of QMB and QFT systems will entail the creation of an extended array of entanglement measures, inspired by the natural symmetries and entanglement fluctuating interactions of the subatomic world.

Though the landscape for quantifying the structure of entanglement in QMB and QFT systems can seem daunting and precarious at times, this connection has already begun and will continue to provide novel perspectives of entanglement…

The first proposal for a universal quantum simulation algorithm in Ref. Lloyd (1996) used the Lie-Trotter Trotter (1959) approximation for the time evolution operator, which requires a number of operations scaling as $\mathcal{O}(T^{2}/\epsilon)$ (per measurement) in order to ensure a maximum error $\epsilon$ in a simulation of time $T$. Higher order Suzuki-type integrators Suzuki (1991) enjoy a better scaling $\mathcal{O}(5^{k}T^{1+\frac{1}{k}}/\epsilon^{\frac{1}{k}})$ with $k\in\mathbb{N}$ the approximation order. Due to the exponentially increasing prefactor, these approximations are typically used only with $k=1,2$ (e.g., Ref. Childs et al. (2018) for a more detailed discussion). More recently, alternative approaches to simulating real-time evolution with a reduced dependence on $\epsilon$ have been proposed, enabled by increased qubit overhead. Important examples are techniques based on the linear combination of unitaries (LCU) Childs and Wiebe (2012) and the Taylor expansion Berry et al. (2015a, b), which achieve an exponential improvement in the error scaling as $\mathcal{O}(T\frac{\log^{2}(1/\epsilon)}{\log\log(1/\epsilon)})$, and also techniques based on Qubitization and Quantum Signal Processing (QSP) Low and Chuang (2017, 2019), which are able to achieve the optimal scaling $\mathcal{O}(c_{T}T+c_{\epsilon}\log(1/\epsilon))$. Better than $\mathcal{O}(T)$ scaling is forbidden by the no-fast-forward theorem unless some of the structural properties of the Hamiltonian are used in the design of the simulation scheme Berry et al. (2007); Atia and Aharonov (2017), i.e. codesign. It is important to point out that we do not expect to be able to implement generic unitary transformations Knill (1995) and that being able to construct a qubit or qudit Hamiltonian describing a physical system is a necessary but not sufficient condition for the existence of an efficient simulation strategy (see, e.g., Ref. Childs and Kothari (2010)).

Besides second quantization schemes, an alternative mapping can be obtained using first quantization language where the particle’s statistics are included by explicit symmetrization Abrams and Lloyd (1997); Aspuru-Guzik et al. (2005). The main advantage of working in this basis presents itself in situations where the number of particles $\eta$ in the system is conserved and discretization over a large number of fermionic modes $N$ is required—in second quantization the required number of qubits is at least $n=N$ (using e.g., the JW or BK mappings described in Sec. V.1.1 below), while in first quantization this can be reduced exponentially to $n=\mathcal{O}(\eta\log(N))$ instead Babbush et al. (2019) when $N\gg\eta\gg 1$. This reduction in Hilbert space complexity is typically accompanied by an increase in time complexity and, for some important situations, the two descriptions have comparable simulation performance (see, e.g., Ref. Su et al. (2020)). Analogous symmetry or conserved quantity projections are also under development in simulations of gauge field theories (see Section V.1.3), and are found to increase classical preprocessing requirements.

More recently, growing interest in reducing the number of quantum gates required for realistic simulations has led to an increasing appreciation of the potential benefits of including stochastic components to simulation algorithms. Notable ideas in this direction are, for example, randomizing the order of exponentials in the Trotter expansion of the evolution operator, resulting in effective higher-order integrators Childs et al. (2019), and utilizing the fully stochastic compilation of Trotter evolution with the quantum stochastic drift protocol (QDRIFT) Campbell (2019) and its extensions Chen et al. (2020); Faehrmann et al. (2021); Berry et al. (2020). Further developments toward realistic simulations have leveraged EFT renormalization group (RG) ideas to prioritize the low energy subspace in dynamical simulations Şahinoğlu and Somma (2020).

Due to the finite gate fidelities of real quantum devices, an inverse relationship exists between the number of qubits in a time-evolved system and the number of Trotter steps feasible within the “gate fidelity coherence time”. As discussed by Martonosi and Roetteler Martonosi and Roetteler (2019), different algorithms for simulations with different target objectives will perform optimally on different architectures. In particular, the scaling of the number of gates with increasing system size for a given algorithm will determine performance on near-term devices with given error rates. The important connection between gate errors and number of qubits for pursuing computation has led to the exploration of more holistic metrics than qubit number or error rate alone e.g., the quantum volume Moll et al. (2018). This metric characterizes the performance of “square” circuits e.g., if a 5 qubit depth 5 circuit provides reliable results (so defined from a random quantum circuit), but a 6-qubit, depth 6 circuit does not, then the quantum volume is $2^{5}=32$.

In light of the possibility of using error correction schemes to encode quantum information in a redundant way allowing for a direct reduction of the effective error rate afflicting a given simulation, it’s important to realize that quantum algorithms requiring seemingly very different logical quantum resources could be implemented successfully on the same general purpose quantum machine. In order to illustrate this point, we show in the right hand panel of Figure 1 three indicative configurations of computational problems (denoted with labels $A$,$B$ and $C$) requiring different logical resources. The gray region indicates instead the regime in logical resources we expect to be available in the near term NISQ era: a few hundred to a thousand logical qubits and a few thousand logical operations. By means of encoding, the problem classes defined at the logical level can be mapped to wide regions in terms of physical resources. Provided the error rate in the physical device is smaller than the threshold for error correction (indicated with a dashed blue line on the left panel of Figure 1), it is possible to trade a larger number of physical qubits for lower error rates. As a concrete example, the configurations of problems denoted as $A$ in Figure 1 could be implemented successfully on both a small quantum device with $\mathcal{O}(10^{2})$ physical qubits and very small error rates $\mathcal{O}(10^{-12})$ or equivalently on a large quantum device with $\mathcal{O}(10^{5})$ physical qubits and much larger error rates $\mathcal{O}(10^{-4})$. Due to the formidable engineering challenges in reducing the limiting noise level to arbitrarily low levels, a commonly adopted strategy in the pursuit of universal quantum computing devices e.g., as discussed by technology companies Kelly (2018), is to first achieve error rates a few order of magnitude below the error correction threshold, and then develop larger qubit arrays that will enable scalable fault tolerance. A representative path for this hardware development is shown as the dotted black line in Figure 1.

While estimates of classical computing resources are well established in many areas, such activities for quantum computing are just beginning. Research teams have undertaken extended periods of resource estimation and coordinated algorithm development tied to present-day simulation capabilities e.g., Refs. Wecker et al. (2014); Reiher et al. (2017); Babbush et al. (2018); Lee et al. (2020a), to reduce naive quantum resource requirements from, in some cases, billions of device-years to a handful of device-days. One of the common standards for performance measures and resource estimation is the Hubbard model and other spin systems with “simple” interactions, e.g., Ref. Childs et al. (2018). At this point is important to highlight the fact that problem instances depicted in both panels of Figure 1 are necessarily narrowly defined problems with clearly stated observables and target precision goals, examples of problems belonging to the 3 configurations depicted are:

• •

  configuration A: computing the ground-state fermion condensate in the Lattice Schwinger model with $N=16-64$ lattice sites and a electric field cutoff $\Lambda=2-8$ using an optimized scheme without amplitude estimation Shaw et al. (2020)

• •

  configuration B: estimating the frequency dependent, inclusive response function with resolutions $\Delta\omega=10-100$ MeV, at one value of the momentum transfer for a system of $A=\mathcal{O}(50)$ nucleons described using Lattice EFT Lee (2009) with a minimal basis defined on a $10^{3}$ lattice Roggero et al. (2020a)

• •

  configuration C: minimal instance of a problem in financial derivative pricing showing a quantum advantage Chakrabarti et al. (2021). The estimate is that a gate depth of $54\times 10^{6}$ on 8K logical qubits will be required.

With rapidly increasing focus on developing quantum algorithms to address scientific applications, and on dedicated and general purpose hardware, benchmarks and estimates of quantum resources requirements are expected to evolve rapidly. In some instances, the estimated resources will increase as our knowledge improves and naïvety recedes, but reductions in resource requirements for quantum advantages in scientific applications are anticipated as methodologies become more sophisticated and devices advance.

The flexibility enabled by error correction mappings inspires the viability of engineering universal quantum devices able to tackle a variety of problems with widely different computational requirements at the logical level. For our three problem configurations, for example, a quantum device characterized by the yellow pentagon in the left panel on Figure 1 would be able to execute instances from all three configurations by using various levels of error correction. On the other hand, we expect near term applications in the NISQ era to use little or no error correction and the gray region of logical resources in the right panel of Figure 1 will be mapped, almost unaltered, onto the physical resources. One complementary route to early success in obtaining a quantum advantage for a problem of scientific interest is the adoption of ideas of codesign to develop special purpose architectures tailored to specific problems ⁶⁶6 In the development of high-performance classical computing in the US, the QCDSP and QCDOC systems Christ (2011) were specifically codesigned for lattice QCD, with on-chip memory and low-latency communication fabric. These systems constituted the precursor to IBM’s successful Blue-Gene series of HPC systems. There were parallel developments in Europe, with the APE machines Marinari (1986). . Codesign of special purpose machines has the potential of drastically impacting the mapping between logical computational resources and physical ones displayed in Figure 1. For instance, it could be possible to design a quantum device with the physical properties denoted by the yellow star in Figure 1, which would be able to perform simulations of problems belonging to configuration $A$ more efficiently but unable to tackle problems belonging to other configurations.

The arrival of heterogeneous classical architectures led to an extended period of time during which scientific applications were ported from CPU-only to CPU-GPU hybrids, and now, most of the performance resides on the GPUs and not CPUs. We anticipate the same sort of transition in the introduction of QPUs into heterogeneous computing Pogorelov et al. (2021). Even assuming that there is a polynomial speed up to parts of the calculations, there is benefit in identifying parts of calculations that would be accelerated by a QPU. In particular, computations of highly entangled subsystems could be performed on QPUs, while the classically efficient elements of the computational workflow, and re-assembly of the results from the QPUs, could be performed with classical computers. Analogous to the widespread adoption of GPU frameworks in scientific computing aided by the standardization of the hardware interface, standardization at both the low-level API and abstracted hardware-specific compilers will provide portability and high-level approaches to quantum algorithmic design. A variety of different framework and standards are currently being explored e.g., Cirq team and collaborators (2020), Qiskit Aleksandrowicz et al. (2019), Qiskit Pulse Alexander et al. (2020), OpenQASM Cross et al. (2017), xacc McCaskey et al. (2020), Q# Microsoft Q# Language Design Team .

Systematic uncertainty quantification generally requires an ensemble of simulations with strategically chosen input parameters such that extrapolation to asymptotic configurations beyond the capability of a single simulation configuration (a hero run) can be achieved. This ensemble maps out a surface in parameter space that enables extrapolations and interpolations to reliably estimate observables of interest and their uncertainties e.g., Ref. Li and Benjamin (2017); Temme et al. (2017); Kandala et al. (2019), injecting another layer of analysis to resource estimation. As is the case for classical calculations and experimental design, designing an optimal manifold of simulations to achieve the desired precision in observables requires planning and simulations of the simulations ⁷⁷7 While hero runs are not the pure focus of simulation, they may be the relevant mode of operation for other quantum applications..

Part of designing the simulation workflow involves identifying and optimally interleaving such calibrations with the physics simulations. While it is obvious that the nature of the calibration is hardware dependent, they will also depend upon the structure of the quantum circuits being executed. Simultaneous adoption of a suite of error mitigation protocols Li and Benjamin (2017); Temme et al. (2017); Endo et al. (2018); Kandala et al. (2019); Dumitrescu et al. (2018) has been used to assess the reliability of the extraction of physical quantities and their sensitivity to the specific noise properties of the device being utilized, e.g. Ref. Roggero et al. (2020a). As NISQ-devices currently resemble experiments more than classical computers, device calibrations that are temporally correlated (integrated) with physics calculations are valuable to in vivo monitor and subsequently identify physics data accumulated during periods in which the device is performing within specified parameters. An example of such “cuts” is provided in Ref. Klco and Savage (2020a), where calibration circuits that measured the in-medium fidelity of a Hadamard gate were included to provide contextualized performance information in initializing symmetric wavefunctions.

Another important aspect of workflow depends upon the (user defined) partitioning of the simulation into classical and quantum components, dependent on the performance of the QPU relative to the CPU+GPU, and the efficiency of communication between them. The current state of the NISQ ecosystem mandates that scientific applications must include a diverse range of integrated controls and diagnostics of the quantum hardware, both accessible and inaccessible to the users. In addition to making available the first cloud-accessible quantum computers, IBM Q Experience IBM provides regular calibration data for each of their superconducting quantum systems, built-in functionality for some error mitigation protocols, and access to pulse-shaping for gate optimization through their Qiskit Aleksandrowicz et al. (2019) programming environment—demonstrating a successful path (and potential template) for co-design collaborations with universities and national laboratories toward scientific applications.

No discussion of quantum resources would be complete without considering complexity classes in the formal sense of computer science, see e.g., Refs. Watrous (2008); Bernstein and Vazirani (1993); Cleve ; com (2013). The scaling of computational resources required to determine a given observable with increasing system size for asymptotically large systems and resources defines the complexity class within which the observable resides. For instance, the P class (PTIME) contains all problems that require computational time scaling as a polynomial of systems size (on a deterministic Turing machine), and is loosely defined by the set of problems that can be solved efficiently. In some cases, the assignment of problems to complexity classes depends on the performance of known algorithms, though not all classifications provide constructive guidance. Because classical complexity classes, defined with respect to Turing machines (deterministic or non-deterministic), do not capture the scaling of problems addressed with quantum devices, additional quantum complexity classes have been introduced. For example, BQP is the complexity class containing problems with solutions of bounded-error requiring polynomial-scaling time using quantum resources. Since a quantum circuit can simulate a classical one, both P, and it’s probabilistic generalization BPP, are contained in BQP. Jordan, Krovi, Lee and Preskill Jordan et al. (2018) have shown that scattering within interacting scalar QFT with external classical sources lies within BQP-complete, indicating that all problems in BQP can be mapped to scattering in scalar QFT with polynomial-scaling time to solution, and further that the scattering problem itself is in BQP and thus efficiently simulatable quantum mechanically.

While it would have been convenient if all scientifically interesting problems were in BPP or BQP, the unfortunate fact ⁸⁸8This fact may be rather fortunate if this complexity is essential to support life. is that they are likely not. Many important problems lie in NP (nondeterministic polynomial time), requiring beyond polynomial resources to solve (assuming P $\neq$ NP) but only P resources to verify solutions. The quantum generalization of NP with polynomially intractable problems is denoted by QMA, and problems belonging to this class are not expected to be solved efficiently even with a quantum computer (unless BQP$=$QMA). Similar to the NP-completeness of finding the ground state energy of a classical Ising spin glass due to an exponentially vanishing gap to the first excited state (see e.g., Ref. Barahona (1982); Istrail (2000); Troyer and Wiese (2005)), the problem of determining the ground state energy of a $k$-body Hamiltonian is QMA-complete for $k\geq 2$ Kempe et al. (2005). A $k$-body Hamiltonian (often referred to as a $k$-local Hamiltonian in the computer science literature, though no sense of spatial locality is intended ⁹⁹9 The most interesting aspects of $k$-local_cs (as defined by computer scientists) systems to physicists are likely the $k$-nonlocal_phys attributes induced by entanglement.) is defined as an operator that admits a decomposition into a polynomial number of terms acting individually at most on $k$-qubits. Interestingly, the $2$-body Hamiltonian problem was shown to be QMA-complete Kempe et al. (2005) also when restricted to spatially local interactions on a 2D grid. These results show that it will not be possible to solve efficiently for ground-state energies of even simple Hamiltonians in general, similar to the NP-hardness of a general solution to the sign problem in Quantum Monte Carlo methods Troyer and Wiese (2005).

It is important to point out that membership to a particular complexity class is a statement about the worst-case instance of a particular problem, and does not necessarily reflect the average- or best-case. A great example of this dichotomy between average case complexity and worst-case complexity is the Minesweeper game: in this game, players are given a 2D lattice with $L^{2}$ sites containing $M$ mines. When a site without a mine is uncovered, the site is assigned a numerical value corresponding to the number of neighboring mines. The goal is to uncover all sites not containing mines using inference from an initial set of uncovered sites. This problem has been shown to be NP-complete Kaye (2000) (or more precisely co-NP-complete Scott et al. (2011)) and therefore extremely hard to win in the worst case scenario. For those who have played Minesweeper, this result may clash with practical experience, in which the game can often be solved very quickly. An important structural property of the Minesweeper game is the mine density $\rho=M/L^{2}$, and it was recently argued that a complexity phase transition occurs once the mine density exceeds a critical value $\rho_{c}$ Dempsey and Guinn (2020). At small mine density, the game is typically easy and inference can be carried out locally on the lattice. As the critical mine density is approached, successful inference of the mine’s location requires consideration of lattice patches approaching the full volume of the lattice, causing solutions to become expensive to find. This “frustration effect” shows up in certain regimes of a class of NP-complete problems, for which a prototypical example is the phase transition in random k-SAT problems Monasson et al. (1999). This discussion draws obvious analogies with chiral symmetry breaking, confinement, and the QCD phase diagram, where transitions are accompanied by the delocalization of natural degrees of freedom capturing the microscopic properties of quarks and gluons.

To effectively connect quantum simulations of SM observables to experiment, a complete quantification of uncertainties is required. Lattice QCD calculations have established well-defined uncertainty quantification methodologies (for a review, see Ref. Beane et al. (2015)), in addition to demonstrating the benefits of independent teams of researchers and code developments, and comparable access to independent computing environments. Lattice QCD is a low-energy EFT of QCD that can faithfully reproduce low-energy QCD observables. Such simulations involve discretizing volumes of spacetime, introducing a lattice spacing, and bounding a spatial extent in each of the four spacetime directions. Therefore, observables will necessarily deviate from their QCD value. In addition, as simulations only involve dimensionless quantities, an overall length scale is required to be determined, along with input values for the quark masses. This is typically achieved by making comparisons with a small number of precisely known experimental quantities, such as hadron masses. Such tunings are imperfect and iterative, and introduce further deviations from QCD values. EFTs, such as $p$-regime chiral perturbation theory ($\chi$PT), play an essential role in post-computational processing of simulation results to make reliable predictions of QCD with fully quantified uncertainties (see, e.g., Ref. Golterman for mesons and Ref. Savage (2006) for nuclei and multi-nucleon scattering systems). For small enough lattice spacings, compared with the QCD scale and external kinematics, the Symanzik action Symanzik (1983a, b) includes the structure of the gauge-invariant operators that may be introduced into the simulation (dependent upon the nature of the field discretization) to compensate for the finite lattice spacing in each direction, or alternately dictate the nature of the power-series in the lattice-spacing(s). While the quark-mass dependence of a number of low-energy observables is known from $\chi$PT, or HQET, calculations are typically performed in the neighbourhood of the physical point (or the mass points of interest) and interpolations are performed to mitigate tuning errors. On top of systematic errors, there are statistical errors and the range of behaviors of observables with sampling. The extraction of energies and observables requires sufficient statistics (sampling of gauge configurations) to recover correlation functions with sufficient precision to be able to fit relevant quantities with controlled estimates of mean values, statistical uncertainties and systematic errors.

For quantum simulations, there are further systematic errors associated with the quantum hardware and the qubit representation of the simulated system. Mapping bosons onto a finite Hilbert space introduces digitization errors and field truncation errors, to be quantified and/or mitigated. There are errors associated with state-preparation (similar to source structure in Euclidean space, which may inform quantum simulations), and measurement, and both are at early stages of being understood and quantified. Other sources of error that must be quantified are device noise and field pixelation (operator smearing) on quantum correlations (i.e., the structure of entanglement between separated regions of a latticized field, e.g. Ref. Klco and Savage (2021a, b)). The systematic and statistical uncertainties associated with a quantum simulation (determined by the device, algorithms and theoretical framework) can be written as,

where $\{\epsilon_{\rm th.},\epsilon_{a},\epsilon_{L},\epsilon_{\theta},\epsilon_{M},\epsilon_{\rm Meas.},\epsilon_{\rm ops.}\}$, associated with the theory, lattice spacing, volumes, parameter tuning, scale setting, measurements and operator structures, are familiar (but some will be handled differently) from classical simulations. The $\epsilon_{\rm digital},\epsilon_{\Lambda},\epsilon_{\rm Trotter},\epsilon_{\rm Noise}$ associated with field digitization and mappings, truncations in field space, Trotterized time evolution, and device noise, are further sources of uncertainties that need to be quantified.

The statistical behavior of an observable, and how it is constructed from ensemble measurements, is another important feature of algorithmic design. With direct access to the wavefunction of a system not practical (tomography scaling exponentially with system size), optimizing observable estimators is essential for efficient use of quantum simulation (as it is classically). Cancellations between contributions to expectation values can create sign problems at the measurement stage of a calculation, exacerbating the number of circuit executions required to achieve a given precision e.g., Ref. Roggero and Baroni (2020). Because the statistical structure of the measurements is mapping (basis) dependent, it is possible that circuit design and algorithms with less favorable asymptotic scaling may provide compensating enhanced statistical convergence.

The complexity classes of SM observables will determine the observables that can be computed with arbitrary precision using combinations of classical and quantum resources. However, if they lie outside P or BQP, it does not mean that significant progress, and potentially all of the required progress, cannot be made in addressing important problems, gaining significant physical insights and even proposing new experiments. Explicitly, some observables formally in P may scale with a sufficiently high-order polynomial that they are unreachable with the requisite precision for foreseeable computers, while some in NP or QMA may scale benignly for small or modest system sizes, enabling computations with “sufficient” precision to have near-term or long-term impact. The impact of any given simulation is limited by its precision, not only of the computation, but also of the theory input into the simulation. Scientific simulations have progressed organically with simulation providing results of precision $\epsilon^{(t)}$ (where the superscript $t$ denotes total uncertainty at the time of the simulation). Once $\epsilon^{(t)}$ is specified, there is freedom to define approximate (perturbatively-close) theoretical inputs and tolerances of the computation, to furnish a result that lies within $\epsilon^{(t)}$ of the true result. So, for bounded-error simulations, with a target tolerance $\epsilon^{(t)}$, as is the case with quantum devices, the input theory (model or EFT) can deviate from the target theory, provided that the induced errors are $\leq\epsilon^{(t)}$. This builds upon the recognition that only a subset of possible quantum observables are required to be simulated e.g., low-energy scattering amplitudes and not amplitudes up to the Planck scale, which leads to the employment of EFTs and phenomenological models ¹⁰¹⁰10 An implicit distinction in motivations for developing simulations of arbitrary Hamiltonians and domain science applications. .

Leveraging lower-complexity leading order interactions has enabled classical simulations to make significant progress in addressing problems formally beyond the reach of classical computing. The SM has many examples of QMB and QFT systems that are amenable to such perturbative expansions centered on lower-complexity leading orders:

• •

  Density Functional Theory (DFT) is a highly successful framework for nuclear structure and reactions (see e.g., Ref. Colò (2020)). The foundational result upon which DFT is based is the Hoenenberg-Kohn theorem Hohenberg and Kohn (1964), which states that the problem of determining the ground state energy of QMB or QFT systems can be replaced by an equivalent optimization problem over one-body density matrices $\rho$ using an appropriate functional $F[\rho]$ as a cost function. Central to the accuracy of the method is the availability of accurate approximations to the functional $F$ since, once this is available, the optimization can be usually carried out efficiently in P. At first glance, this seems at odds with the observation that estimating ground state energies of even simple $2$-body Hamiltonians is QMA-complete. Indeed, it is possible to show that, even though a DFT calculation is usually tractable once the functional $F$ is known, finding an accurate functional is a QMA-complete problem Schuch and Verstraete (2009). However, with the guidance of physical understanding and experiment, approximate functionals of sufficient quality to reproduce observables like binding energies within a given accuracy can often be found—nature hacking complexity. The QMA classification in this case is thus a statement about the intrinsic difficulty of designing functionals capable of arbitrary precision in general.

• •

  Low-energy nuclear forces approximately manifest Wigner’s emergent SU(4) spin-flavor symmetry Wigner (1937a, b) (which also emerges from large-N_c considerations Kaplan and Savage (1996), vanishing entanglement power considerations Beane et al. (2019), and from lattice QCD calculations Wagman et al. (2017)). Quantum Monte Carlo simulations of multi-nucleon system ground states, e.g., using the pionless EFT, generally exhibit a sign problem. Because the sign problem vanishes in the limit of SU(4) symmetry, successful hybrid schemes exist involving “pre-conditioning” systems with the SU(4) symmetric interactions followed by implementation of the full interaction, mitigate the impact of the sign problem Lee (2016). This process identifies a leading order part of the QMA calculation that is in BPP, and perturbatively estimates the full result to fixed precision.

• •

  Nuclear structure has been tremendously impacted by Green’s Function Monte Carlo (GFMC) calculations of light nuclei (see Ref. Carlson et al. (2015) for a review). Despite requiring exponentially increasing resources with baryon number, for an extended period of time, GFMC was able to provide accurate results for system sizes not reachable with polynomially-scaling algorithms, such as Coupled Cluster. This balance of practical computational power has almost entirely reversed in recent years due to the adoption of low-resolution interactions (derived from chiral-EFT Weinberg (1990, 1991); Ordonez et al. (1996); Kaplan et al. (1996, 1998a, 1998b); Epelbaum et al. (1998); Machleidt and Entem (2011)) and RG techniques Kaplan et al. (1998a, b); Bogner et al. (2003) (see, e.g., Ref. Hergert (2020) for a recent review).

• •

  Current lattice QCD calculations indicate that:

  • –

    Computing the mass of the pion to a given precision can be accomplished efficiently at scale, consistent with residing in BPP

  • –

    The exponentially degrading signal-to-noise at late times in nucleon and multi-baryon correlation functions Parisi (1984); Lepage (1989); Beane et al. (2008) (a sign problem) currently suggests Troyer and Wiese (2005) that the Quantum Monte Carlo computation of the binding energy of nuclei lies outside of P and BPP.

  Through its connection to the Hubbard model, finding the ground state of the pionless EFT is among the hardest problems in QMA Schuch and Verstraete (2009); O’Gorman et al. (2021), and therefore, it is likely that lattice QCD calculations of nuclei lie in QMA or beyond in PSPACE. However, calculations of energies and properties of smaller nuclei, which have been demonstrated Beane et al. (2013), permit the bounded-error extraction of counterterms in nuclear EFTs Barnea et al. (2015); Kirscher et al. (2015); Contessi et al. (2017); Bansal et al. (2018). These results may be then used to compute the energies, properties, and low-energy dynamics of somewhat larger nuclei Savage (2006) (for reviews on this subject, e.g., Ref. Beane et al. (2008, 2011); Savage (2016); Drischler et al. (2019); Davoudi et al. (2021a); Aoki and Doi (2020)).

• •

  Heavy-quark systems (exhibiting heavy-quark symmetry for $m_{Q}\rightarrow\infty$) are comprised of one or more heavy quarks, $Q$, and light degrees of freedom. Leading order wavefunctions of $Q\overline{l}$-mesons are a tensor product, $|Q\overline{l}\rangle\sim|Q\rangle\otimes|\overline{l}\rangle$. The higher-order $1/m_{Q}$ contributions to observables can subsequently be included in perturbation theory from this wavefunction. Powerful theorems, such as the Adellmolo-Gatto theorem Ademollo and Gatto (1964) and Luke’s theorem Luke (1990), ensure that leading order results in the symmetry breaking parameters are known without model dependence, and that symmetry breaking in certain matrix elements is suppressed by two powers of the breaking parameter. Further, any model dependence or parameters that need to be computed elsewhere enter at higher orders in the perturbative expansion(s), and associated uncertainties are thus parametrically smaller than leading order.

• •

  The parton model characterizes hadronic structure and processes at high momentum transfer. The scale invariance of QCD in this regime, due to asymptotic freedom at short distances, allows many QCD processes to factorize at leading order in pQCD, providing a foundation for the systematic inclusion of higher order corrections. In such situations, classical computing has been proven effective in computing higher-order matrix elements with the required precision to advance experimental programs.

For quantum simulations, this hacking of bounded-error complexity classes through perturbative changes to the theory is only now beginning. For example, precision first principles calculations of fragmentation lie beyond the reach of classical computing. It is also expected to lie outside of BQP, due to final state complexities, but perturbatively addressable with simulations that lie within BQP. The use of Soft-Collinear EFT (SCET) Bauer et al. (2000, 2001); Bauer and Stewart (2001) may enable simulations within desired error bounds (though this remains to be demonstrated) Bauer et al. (2021a). For such hacks to be successful, problems should be broken up into a large part that resides within BQP and can be performed efficiently on a quantum device, and a perturbative part (containing the terms that place the problem outside of BQP).

Remarkable progress has been made toward understanding nature through the SM even for problems that reside in complexity classes formally beyond the capability of classical or quantum computers. This progress has been enabled, to a large degree, by points of enhanced symmetry exhibiting reduced complexity. In such instances, nature has been “kind” to us by providing Hamiltonians that can be separated into “large” and “small” parts for the observables of interest. This separation has been enabled by the identification of approximate symmetries of the systems, with the “large” part(s) invariant under the symmetries, and the “small” part(s) breaking the symmetries. Generally, and in hindsight, these symmetries are identified by eliminating entanglement and identifying perturbatively close tensor-product states. However, examples of expanding about leading-order entangled states are known Kaplan et al. (1998a, b); van Kolck (1999), resulting from a non-trivial fixed-point in RG flow. For classical computing, the challenge has been to find symmetries that place leading order approximations into P or BPP with a controlled perturbative expansion that sufficiently converges. For quantum computing, the leading order symmetries need not eliminate entanglement. Instead, incorporation of quantum devices will allow the advantage of addressing important QMA-complete SM applications from the standpoint of leading order approximations residing in BQP. A nice example of this strategy can be found in low-energy nuclear structure where nucleon-nucleon correlations can be self-consistently included into a more effective single nucleon basis Robin et al. (2021), perhaps distilling the leading-order consequences of entanglement from a QMA-complete problem into one that is in BQP. Broadly, identification of the leading order entangled systems of enhanced symmetry presents an exciting challenge.

At this early stage of development of quantum simulation, asymptotic scaling does not conclusively determine the reach of present day or near-term calculations, including beyond the NISQ-era into the period where error-correction becomes routine. While a better understanding of the formal asymptotic scaling of classical and quantum resources required to address challenges facing the SM is important, it will not be the sole consideration in identifying goals for near-term quantum simulation, as important progress may be achievable, possibly even with demonstrations of quantum advantages, for problems that lie outside BQP.

The mapping one chooses to represent a quantum system of interest with quantum computational degrees of freedom will dramatically impact the quantum resources required. In particular, the representation will affect the cost of all steps in a quantum simulation from the entanglement structures created in initial state preparation, to the interactions needed to perform time evolutions, to the distribution of information desired to be probed in the final measurement procedure. It is expected that the basis-dependent property of entanglement may guide the development of intuition in mapping fields onto quantum degrees of freedom (see Section II.2.6). Furthermore, exploring the repercussions of various representations informs theoretical flexibility to codesign simulations.

A basic question to consider in quantum simulation design (that has begun to be addressed in the application to fields) is Where is the Line, referring to the flexible partitioning of computational responsibility between the quantum device and the large-scale classical resources to which it will be symbiotically married. One category of explorations shifting this line has been in state preparation e.g., through the use of quantum devices to optimize the design of interpolating operators for classical LGTs Avkhadiev et al. (2020), through the use of classical snapshots of a volume set by the correlation length scale to inform state preparation on quantum volumes capturing many correlation length scales Klco and Savage (2020b, c), or through purely quantum techniques inspired by the structure of classical tensor network design Moosavian et al. (2019). A second category of explorations shifting the line is in the basic representation of the field. The influential papers of Jordan, Lee, and Preskill (JLP) illustrated an efficient approach for expressing real-time dynamics of scattering events quantum mechanically Jordan et al. (2012, 2014). With polynomial resources in particle number, energy, precision, and volume, a microscopic description of interacting scalar field theory represented by local qubit degrees of freedom was shown to produce final state probability distributions that reflect, to systematically controlled precision, the exact probability distributions that would be experimentally observed. Further inspiring the use of quantum devices for this purpose, recent quantification of computational complexity in tensor network simulation of inelastic false-vacuum bubble scattering illustrated a growth of entanglement at high energies and multiple scatterings Milsted et al. (2020), clearly connecting quantum effects to the voracious consumption of classical computational resources in scattering processes. As algorithms for the quantum implementation of microscopic time evolution operators continue to progress, a complementary perspective from EFT descriptions is under development Provasoli et al. (2019); Bauer et al. (2021b, a), utilizing quantum devices to address non-perturbative physics isolated in EFT descriptions of gauge field dynamics e.g., of parton showers, high multiplicity final states, or soft functions of SCET. As a reliable methodology for systematically isolating energy regimes, the demonstrated synergy between EFTs and microscopic quantum simulations of fields is expected to reduce quantum resource requirements by multiple orders of magnitude and offer a clear route for integration of quantum devices into existing computational frameworks of experimental impact. Conveniently, the dynamical techniques developing for the microscopic simulation of gauge theories are closely related (often identical) to those leveraged for the calculation of real-time non-perturbative matrix elements within EFT descriptions. As such, explorations considering different quantum-classical partition lines may be explored in parallel, to mutual benefit, without concern of obsolescence dependent on the trajectory of quantum hardware development.

As a purist motivation for focusing on the quantum simulation of fields, it is expected that the quantum mechanical complexity at play in atomic scale and larger systems will be commensurate with that in subatomic and smaller systems. Thus, precision control of the former may enable efficient calculation of emergent properties of the latter, in particular with resources that no longer scale exponentially with volume. A concrete example of this commensurate scaling can be seen in the identification of vacuum-to-vacuum matrix elements in a self-interacting scalar field with classical external sources to be a BQP-complete problem Jordan et al. (2018). In this direction, the relevance of fields for quantum simulation is not limited to relativistic fields describing the dynamics of the building blocks of the SM, but also applies to non-relativistic systems, such as EFT descriptions of low-energy physics. From a less-purist perspective, the variety of qubit arrays that are being developed in laboratories around the world e.g., 2-dim arrays of implanted atoms in instrumented silicon substrates, ultracold atoms in optical lattices, ion lattices, superconducting circuit ladders and arrays of Josephson junctions, are the closest physical systems to latticized quantum fields that have been created to date. As such, the BQP-completeness of the scalar field Jordan et al. (2018) potentially indicates a natural language for quantum computation—any problem (physics or beyond) that can be efficiently calculated on a quantum device can be efficiently mapped to a scattering problem in scalar field theory, which itself is efficient to implement on a quantum device Jordan et al. (2012, 2014). One step further lies the connection with quantum error detection/correction, where the quantum information used to compute is embedded in topological excitations of a gauge field or a bulk holographic dual Almheiri et al. (2015); Pastawski et al. (2015), both providing a redundancy to enable robustness to local environmental perturbations. So, whether one is interested in fundamental simulations of nature or large-scale quantum computations more broadly, the QFT language provides a natural framework.

Typically considered to be in the realm of condensed matter, quantum spin liquids represent an important connection between research areas, bringing together QFTs, spin models and the distribution and protection of quantum information (for a review, see Ref. Savary and Balents (2016)). While there are a number of candidate systems to provide logical qubits and quantum memories, as reviewed in Ref. Campbell et al. (2017), Kitaev’s Toric Code Kitaev (2003) (and subsequent surface codes exploring alternate boundary conditions Kitaev (1997b); Bravyi and Kitaev (1998); Dennis et al. (2002); Fowler et al. (2012)) represents perhaps the first and clearest connection between these areas of research. In its simplest form, a 2-dim system of spin-${\frac{1}{2}}$ spins defining the links of a square lattice is subject to a Hamiltonian with vertex operators (stars), plaquette operators, and periodic boundary conditions. The ground state of the system is topologically ordered and consequently highly entangled, supporting two logical qubits defined by the action of Wilson-loop operators in the two directions. These topological states are robust against errors (local perturbations) below some threshold, which can be corrected by measuring $n-4$ independent stabilizers (commuting star and plaquette operators) and using decoders to determine an effective set of operations to restore the system. Reducing the number of independent stabilizers allows an increased dimension of the code space at the expense of error robustness. The quasi-particle picture is a convenient way to understand the behavior of these systems and errors. Localized errors create quasi-particles in the system and error correction corresponds to annihilating quasi particles, returning the systems back into the code space (ideally without inadvertently producing a logical error). Excitations above these ground states are (non-local) pairs of electric charges and pairs of magnetic vortices. The anyonic nature of these excitations and their braiding points to the high-degree of underlying entanglement of the ground state(s), and as such, local operations are unable to resolve which state of the (logical) ground-state space the system is in. An auxillary spin can be included at each lattice vertex to transform this system into a $Z_{2}$ LGT, of similar form as the Kogut-Susskind Hamiltonian for SU(N) LGTs (for a review, see Ref. Kitaev and Laumann (2009)), and generalizing in detail to SU(N) gauge theories Mathur and Sreeraj (2016). The Toric Code is also established in 3-dim and 4-dim using hypercubic (and other) Euclidean space lattice geometries. A major challenge to the execution of circuits on logical qubits is the required action of gates across the physical qubits defining the logical qubits. This development is mature for the Toric Code, see e.g., Refs. Dennis et al. (2002); Jochym-O’Connor and Yoder (2021), and also for other encodings, such as color codes Bombin and Martin-Delgado (2006).

Interestingly from an SM physics perspective, higher-dimension SU(N) LGTs, due to the phenomena of confinement manifest in non-Abelian gauge theories, have utility in defining logical qubits due to their innate mechanisms to localize qubit errors within “hadrons” Bombín (2015); Brown et al. (2016a) (a similar effect can be induced by Anderson localization Wootton et al. (2010)). Unlike the Toric code, where electric charges and magnetic vortices are unconfined and error correction requires multiple rounds of stabilizer measurements, confinement keeps SU(N) Gauss-law-violating vertex charges induced by qubit errors from propagating “too far” from each other, restricting them to approximately within the confinement volume. A possible exception to this is in instances where multiple color-singlet “hadrons” can be formed from the SU(N) charges and wrap “around the world”, an effect that is suppressed at small lattice spacings and is assessed to be topologically benign due to confinement. As a result of this connection, the techniques and technology that continue to be developed for classical LGTs are likely to be of direct impact upon logical qubits and error-correcting protocols, e.g., Ref. Andrist et al. (2011); Brown et al. (2016b).

Another lattice system that has been developed to define logical qubits is the 2-dim honeycomb lattice with spin-$\frac{1}{2}$ spins at each vertex and with link-direction-dependent couplings between spins Kitaev (2006). This system can be written in terms of Majorana fermions with both gapped and gapless phases. Imposing a magnetic field on the gapless phase induces a gap and non-trivial topological structures with chiral edge states carrying Chern numbers. Edge states, and their deep connection to topological structure are commonly used to simulate chiral fermions in LGT calculations—domain wall fermions employ a 5th dimension to isolate a fermion of definite chirality at the domain wall on the 4-dim boundary Kaplan (1992, ).

While efficient non-dynamical approaches to state preparation can be determined for systems with sufficient structure, the ability to simulate the real-time dynamics of a quantum system is commonly a useful resource in studying static properties, e.g., the ground-state energy density. There are three main classes of approaches to attain this goal: adiabatic state preparation, the use of energy filters, and variational approaches. Using adiabatic evolution, it is possible to drive a quantum state from the ground-state of an initial Hamiltonian $H_{i}$ to the ground-state of a final Hamiltonian $H_{f}$ of interest Farhi et al. (2000); Albash and Lidar (2018). Given a procedure to initialize a quantum-computer in the ground-state of $H_{i}$, an appropriate adiabatic path can be designed to produce a final state after time evolution with a large overlap onto the target ground state. The scale for the evolution time $T$ is set by the smallest energy gap $\Delta E$ between the instantaneous ground-state and first excited state throughout the adiabatic path. Several strategies have been proposed to optimize the efficiency of this scheme, ranging from explicit coupling to an external “bath” of degrees of freedom Kaplan et al. (2017); Lee et al. (2020b); Metcalf et al. (2020), gap amplification techniques Somma and Boixo (2013), and specially designed adiabatic paths that leverage previous knowledge of the excitation spectrum to avoid closing gaps during the adiabatic evolution (e.g., Refs. Wecker et al. (2015); Du et al. (2021)). This method has been applied to lattice field theory simulations Jordan et al. (2012, 2014) and applications to the nuclear many-body problem have recently started to be developed Du et al. (2021).

Quantum Phase Estimation (QPE) Kitaev (1995); Abrams and Lloyd (1997, 1999) or it’s more recent variants Svore et al. (2014); Wiebe and Granade (2016), is a candidate technique for computing ground-state properties with a quantum computer at scale. This versatile approach can be used to measure the eigenvalues of the Hamiltonian operator $H$ by exploiting the ability to simulate real-time evolution under the operator $U(t)=\exp(-itH)$ to implement an energy filter of resolution $\delta H=\mathcal{O}(1/T)$, with $T$ the longest time interval required for the algorithm. For sufficiently long propagation times $T$ the probability to measure the ground-state energy when performing QPE on an initial trial state $|\Psi_{T}\rangle$ is proportional to it’s overlap $|\langle\Phi_{0}|\Psi_{T}\rangle|^{2}$ with the ground-state Ovrum and Hjorth-Jensen (2007). It is therefore crucial to be able to initialize the quantum device in a high-overlap trial state $|\Psi_{T}\rangle$, while keeping the complexity of this initialization low for ease in simulating the real-time dynamics of QPE within the limited coherence time of the devices. These trial states are typically obtained using the variational approaches described below, but both adiabatic state preparation or a preliminary coarse energy filter could be employed. Similar to the calculation of dynamical response functions, polynomial expansion methods using Qubitization/QSP can be employed to define alternative unitary operators to use in place of $U(t)$ Poulin et al. (2018) or directly to implement more efficient energy filters for ground-state projection Lin and Tong (2020). Recent proposals for ground-state preparation algorithms like the Rodeo algorithm Choi et al. (2021) or simulated imaginary-time propagation Turro et al. (2021) also belong to the class of energy filters, though their resource scalings are currently less understood.

Variational approaches address the ground-state preparation problem through an (in general non-convex) optimization. This is achieved by constructing quantum states using a fixed quantum circuit with free parameters and using an appropriate variational principle to find a good approximation of the ground-state by optimization. Typically the expectation value of the energy is used as a cost function in the optimization e.g., the Variational Quantum Eigensolver (VQE) Peruzzo et al. (2014); McClean et al. (2016); Sim et al. (2019), which directly optimizes over the parameters of a general unitary operator. Alternatively, the Quantum Approximate Optimization Algorithm (QAOA) Farhi et al. (2014) mimicks the adiabatic state preparation described above, but optimizes a parametrization of the adiabatic path with the goal of reducing the dependence on energy gaps in the spectrum. Other cost functions have also been considered e.g., the action generating imaginary time evolution McArdle et al. (2019) or few-body correlations Motta et al. (2019). The two main advantages of variational approaches are the relatively small number of quantum operations required (at least for simple parametrizations) and the inherent robustness to miscalibration of the quantum device and to noise. An intuitive way to understand this last feature is to realize that the various noise sources will effectively shift the location of the global minimum for a given parametrization, but this effect may be partially compensated by the optimization procedure O’Malley et al. (2016). Variational algorithms are heuristic in nature. On one hand, it would be classically hard to simulate the measurement statistics given a sufficiently complex circuit structure. On the other hand however, finding the optimal parameters using non-convex optimization is in general NP-hard and we do not expect to be able to reach the globally optimal set of parameters for all instances of the problem Bittel and Kliesch (2021). The NP-hardness of these variational algorithms is a natural occurrence not necessarily specific to quantum algorithms e.g., the corresponding NP-hardness of finding the optimal mean-field solution using the Hartree-Fock method Schuch and Verstraete (2009), and it does not necessarily imply a limitation in the applicability of variational algorithms in practice (see Sec. IV), NP-hardness is a statement about the behavior in worst-case instances and is not necessarily indicative of the average complexity, especially if additional structure is present in the problem.

Given their versatility, limited use of quantum resources and resilience to device noise, Variational Quantum Algorithms (VQA) are a flourishing area of research (see Ref. Cerezo et al. (2020) for a recent review). VQA will likely be important as a stepping stone towards preparing high-accuracy trial states to be used as initial conditions for more sophisticated quantum algorithms. This is similar to the use of Variational Monte Carlo as a precursor to more accurate, and expensive, GFMC calculations of nuclei (see Ref. Carlson et al. (2015) for a review). In addition to such compounding applications, another possible use of VQA is as components for Quantum Machine Learning Biamonte et al. (2017) tasks (possibly coupled to classical neural networks e.g., Ref. Filipek et al. (2021)), such as characterizing the output states generated by a quantum algorithm or even directly as a way to perform inference on classical data (for a recent review of applications to SM problems, see Ref. Guan et al. (2021)). Presently, it is not clear if an exponential advantage could be obtained for the latter task (see, e.g., Refs. Aaronson (2015); Wiebe (2020)).

A dominant component of the cost in variational techniques comes from the number of repetitions required to estimate the cost function (and/or it’s gradient) used for VQE, and considerable effort has been devoted in the recent past to design approaches to reduce this cost. These include optimization of commuting sets of operators that can be measured simultaneously and techniques using information from short-time dynamical evolution Roggero and Baroni (2020); Di Matteo et al. (2021); Verteletskyi et al. (2020); O’Brien et al. (2019); Sim et al. (2019); O’Brien et al. (2021); Funcke et al. (2021); Huggins et al. (2021); Guzman and Lacroix (2021). One last class of techniques developed recently for efficient estimation of general expectation values uses random sampling before measurement to drastically reduce the number of measurements required for convergence Brydges et al. (2019); Elben et al. (2019). For example, by informing a ’classical shadow’ Aaronson (2018); Huang et al. (2020) using the results of random measurements, the expectation value of $M$ distinct observables can be estimated using only $\mathcal{O}(log(M))$ measurements on a quantum device.

Extracting reliable dynamical information from QFTs and QMB systems is an extremely challenging task due to sign problems that generically plague real-time path integrals. This difficulty is particularly impactful in the SM since the vast majority of experimental information is obtained through scattering processes. Progress continues to be made in developing more efficient implementations of time evolution for QFTs and QMB systems, as outlined in Section III, but there is a further level of sophistication that is required in designing time-dependent calculations performed on a device to optimally address the physics problems at hand.

A common approach to evaluate scattering cross sections in the linear response regime using classical computing is to first compute Euclidean-time correlation functions, readily available from imaginary-time path integrals, and then attempt a numerical inversion to real-time. This process is ill-posed, in the sense that small errors in the Euclidean correlators can lead to large errors in the result of the inversion, and various heuristic approaches have been proposed to extract dynamical properties with minimal bias (see e.g., Refs. Magierski et al. (2009); Burnier and Rothkopf (2013); Lovato et al. (2016)). A similar approach popular in few-body nuclear physics is to attempt a numerical inversion of integral transforms that are better behaved than the Laplace transform e.g., the Lorentz Barnea et al. (2010) or the Sumudu integral transforms Roggero et al. (2013), with the goal of simplifying the inversion procedure. As we will see further below, similar ideas will play an important role in proposed quantum algorithms for real-time dynamics.

Using the efficient simulation strategies presented in the previous section, quantum computers are expected to be able to simulate efficiently the real-time dynamics of systems with local interactions. As perhaps the most straightforward application of these ideas, direct simulations of Hamiltonian dynamics starting from a reference state have recently appeared e.g., the probability of pair production and chiral dynamics in the Schwinger model Martinez et al. (2016); Klco et al. (2018); Kharzeev and Kikuchi (2020), the three-body contact density in a simple model for the triton Roggero et al. (2020a), the time-dependent Coulomb excitation of the deuteron colliding with an heavy-ion Du et al. (2020), the spin dynamics of a pair of nucleons interacting through a tensor interaction Holland et al. (2020) and the entanglement and flavor evolution in collective neutrino oscillations Yeter-Aydeniz et al. (2021); Hall et al. (2021).

Particularly important for comparison with experiments is the calculation of scattering cross sections in the linear response regime, where the interaction between system and probe is well approximated by a single vertex. Apart from kinematical factors, the full energy dependence of the cross-section is contained in the the response function $R(\omega)=FT\left[\langle O(0)O(t)\rangle\right]$. In this expression $FT$ denotes the Fourier transform and $\langle O(0)O(t)\rangle$ is the real-time two-point function of the vertex operator $O$ evaluated on the ground state of the target. A direct calculation of this two point function is possible in both digital quantum devices Somma et al. (2002) and more general quantum simulators Baez et al. (2020). In practice, the Fourier transform is performed using the two-point function signal up to a maximum time $T$, leading to a resolution in the frequency domain $\Delta\omega=\mathcal{O}(1/T)$. Since the dominant cost of the digital algorithm is the implementation of the real-time evolution operator for a total time $T$, the final cost of these types of approaches scale as $\mathcal{O}(1/\Delta\omega)$. This agrees well with the no-fast-forward theorem mentioned above, and doesn’t exclude the possibility of an important speed up if structural information of the system’s Hamiltonian is used in the design of the algorithm. It is important to remember that the Minkowski-space correlation functions extracted from any finite volume quantum computation need to be properly extrapolated to the continuum. Contrary to their Euclidean time counterparts, e.g. Refs. Lüscher (1991, 1986a, 1986b), these real-time correlators have been shown to pose additional challenges in this extrapolation for certain kinematical regimes, and practical strategies for mitigating the problem have been proposed Briceño et al. (2021).

A more direct approach to scattering matrix elements, inspired by the integral transform techniques used in few-body nuclear theory, can be formulated to work directly in the frequency domain avoiding the approximation of the Fourier transform in both the UV and IR Roggero and Carlson (2019). The main idea is to apply an approximate projector into a frequency band, in the same spirit as the Lorentz Integral Transform method of Ref. Barnea et al. (2010), followed by a direct measurement of the population in that band to reveal the response function strength. Besides a better control of the frequency-space errors, this type of approach opens the way for characterization of semi-exclusive information about the final state of the reaction since, after the application of the energy filter, the register of qubits representing the target is mapped to the set of final states. This is an important step towards a more complete characterization of semi-exclusive cross sections e.g., those required to extract neutrino parameters from long baseline experiments like DUNE, and first estimates of the quantum resources required for this goal have appeared Roggero et al. (2020a), suggesting further improvements are needed for near term applications on real devices. The scheme can be generalized by considering different frequency projectors constructed either using the real-time unitary operator (see e.g., Refs. Somma (2019); Lu et al. (2021)) or more general expansions in orthogonal polynomials Roggero (2020); Rall (2020), which allow for superior scaling with the target precision.

The close connection between spin models that can be used to define logical qubits and LGTs (see Section V.2) might suggest that we dedicate much of this section to discussing results related to the Toric code, quantum spin liquids Savary and Balents (2016) and related systems. Instead, we focus on hardware implementations explicitly targeting the SM. The first digital quantum simulation of real-time dynamics in a LGT Martinez et al. (2016) not only demonstrated the progress that had been achieved in the control and manipulation of quantum hardware, but provided tangible inspiration that the vision of digital quantum simulation for microscopic descriptions of nature was becoming a reality. Complementary to larger U(1) simulations on analog quantum devices and classical simulations, e.g., Refs. Bañuls et al. (2013); Buyens et al. (2017); Bañuls et al. (2016); Funcke et al. (2020); Yang et al. (2020), this digital experiment Martinez et al. (2016) was performed on four trapped calcium ions representing the fermionic degrees of freedom at two spatial sites of latticized 1+1 dim quantum electrodynamics (QED), also known as the lattice Schwinger model. First proposed in Ref. Hauke et al. (2013), this work showed the dynamical generation of $e^{+}e^{-}$ pairs emerging from the trivial vacuum under time evolution, as shown in Figure 2 ¹¹¹¹11 In each figure panel of this section, icons from Ref. Klco and Savage (2020a) are assigned to indicate the type of compute device that generated the results, i.e. a classical computer without a noise model, a classical computer using a noise model, or a quantum device. These icons and their descriptions can be found at https://iqus.uw.edu/resources/icons/. .

Encouragingly, even with the $\sim 50$ quantum gates necessary for each step in the time evolution, control and coherence was demonstrated up to four Trotter steps with a survival probability in the physical, zero-charge subspace (six of the $2^{4}$ fermionic configurations) of greater than 70%, allowing for error-mitigating post-selection while maintaining sufficient statistics. As a theory in one spatial dimension, the Schwinger model is a system where a number of dynamical attributes can be addressed with precision using modern analytical and numerical tools, such as tensor networks, e.g., Refs. Bañuls et al. (2013); Buyens et al. (2017, 2016); Pichler et al. (2016); Butt et al. (2020); Bañuls et al. (2020). However, continuing to perform quantum simulations of this and similar theories aids the benchmarking of quantum devices, and the development of algorithmic techniques.

While 1+1 dim QED is not of direct interest to fundamental physics, it shares attributes with QCD, such as charge screening, a fermion condensate and a spectrum of composite particles (“hadrons”), including two-body and three-body bound states (“nuclei”). Of particular note, it has non-trivial topological structure through a $\theta$-term, proportional to the electric field at non-zero coupling, and also a chiral chemical potential Kharzeev and Kikuchi (2020). Early studies of the impact of time-dependent $\theta(t)$ and $\theta$-quenches have been undertaken using classical simulations of quantum devices Kharzeev and Kikuchi (2020); Zache et al. (2019) (see also related work with quantum link models Huang et al. (2019)), shown in Figure 2. However, in contrast to 3-dim gauge theories, the U(1) gauge field of the Schwinger model, like any gauge field in 1+1 dim, is not dynamical; constrained by Gauss’s law at each lattice site, it is dependent only on the boundary conditions, which conveniently allows an entirely fermionic representation. Since its first implementation, an array of simulations and proposals for simulations (e.g., Refs. Zohar et al. (2017); Muschik et al. (2017); González-Cuadra et al. (2017); Klco et al. (2018); Bañuls et al. (2020); Luo et al. (2020); Aidelsburger et al. (2021)) of U(1) gauge theories using different types of quantum devices have been established, see Figure 2. Though the long-range interactions arising from the non-dynamical gauge field are amenable to trapped ion implementation Hauke et al. (2013); Martinez et al. (2016); Davoudi et al. (2020), the utilization of superconductors inspired the exploration of an alternate approach of classical pre-processing for global and local symmetry projections Klco et al. (2018). This early example served to emphasize the extent to which the design of a digital quantum simulation of LGTs must consider the physics of the underlying quantum device for near-term progress.

First considerations have been made in determining the Schwinger Model’s asymptotic scaling Shaw et al. (2020), and further detailed studies continue to provide important information, e.g., Refs. Kaplan and Stryker (2020); Stryker (2019); Magnifico et al. (2020); Davoudi et al. (2020); Stryker (2021); Halimeh et al. (2020a, b); Van Damme et al. (2020); Davoudi et al. (2021b), including detecting and correcting violations of Gauss’s law and mitigating the effects of evolution into gauge-variant parts of the hardware Hilbert space. While these and future studies reveal important aspects that will be encountered going forward, the fact that the gauge field is not dynamical limits the implications for simulations of the SM. Formalisms for simulating QED in higher dimensions have been put in place for a number of systems, e.g., Refs. Zohar and Reznik (2011); Zohar et al. (2012); Tagliacozzo et al. (2013); Zohar et al. (2013c); Wiese (2013); Marcos et al. (2014); Kuno et al. (2015); Bazavov et al. (2015); Kasper et al. (2016); Brennen et al. (2016); Kuno et al. (2016); Zohar et al. (2017); Kasper et al. (2017); González-Cuadra et al. (2017); Ott et al. (2020); Paulson et al. (2020a); Kan et al. (2021); Aidelsburger et al. (2021), but remain to be executed on a quantum device.

The last several years has seen a number of exciting developments in establishing frameworks for simulations of non-Abelian gauge theories, starting with the pioneering work of Byrnes and Yamamoto Byrnes and Yamamoto (2006). That work builds upon the Kogut-Susskind Hamiltonian Kogut and Susskind (1975) for SU(N) Yang-Mills and first suggestions for implementation by Zohar, Cirac and Reznik Zohar et al. (2013b), which in turn built upon an extensive literature on Hamiltonian formulations of Yang-Mills gauge theories and QCD, e.g., Refs. Wilson et al. (1994); Heinzl (1995); Bronzan (1985); Ligterink et al. (2000a, b). While there are a number of formulations of non-Abelian gauge theories for hardware implementation (see Section V.1.3), the real-time dynamics (on small lattices) of only one construction has so far been simulated on quantum devices, using IBM’s quantum systems Klco et al. (2020); Ciavarella et al. (2021) and a D-wave annealing device A Rahman et al. (2021), as shown in Figure 2. The local gauge invariance at each lattice site has been exploited to implement the non-Abelian constraints explicitly by integrating over the gauge group Bañuls et al. (2017); Klco et al. (2020); Ciavarella et al. (2021). This reduces the dimensionality of the hardware Hilbert space to capture only the irrep of each link, and introduces controlled-plaquette operators. Simulations that have been performed on quantum devices have employed global bases for one and two plaquettes. With the knowledge gained from these small systems, performing simulations with more than two plaquettes in higher dimensions is a priority (for a recent review see Ref. Paulson et al. (2020b)).

Nuclear structure and reactions reside in an important low-energy sector of the SM, with precision targets expected to become achievable with future quantum simulation capabilities. As illustrated in Figure 3, the first ground-state calculation of a nuclear system was a VQE simulation of the deuteron using a harmonic oscillator basis Dumitrescu et al. (2018) with the superconducting quantum devices of IBM and Rigetti. Extensions of the approach soon followed with implementations of the same model on trapped-ion quantum computers Shehab et al. (2019), explorations of different mode-to-qubit encodings using the Gray code Di Matteo et al. (2021), and calculations of the energies of larger systems with $A\leq 4$ Lu et al. (2019). These latter calculations, in large Hilbert spaces containing up to 68 states, were made possible by the use of an all-optical quantum device called a Quantum Frequency Processor (QFP), which allows for an entire calculation to be formulated using a single $d$-dim qudit whose state is manipulated using standard quantum-optics elements. The depth achievable with these devices is unfortunately limited by photon loss and, more generally, the qudit encoding employed requires resources scaling exponentially with the number of orbitals. However, specialized photonic systems like the QFP may have a role in the global ecosystem of quantum technologies for SM physics. The use of variational approaches based on imaginary-time evolution like QITE was also studied recently in Ref. Yeter-Aydeniz et al. (2020). This class of techniques, at least in principle, could be helpful not just as a way of preparing approximations to ground-states (by working at low temperatures) but also open the possibility of studying finite temperature properties in systems like neutron star matter.

Particularly interesting is the possibility to study real-time properties of nuclear systems in thermal equilibrium e.g., transport coefficients in strongly correlated nuclear matter. In general, these calculations will require knowledge of Minkowski-space correlation functions, which are not easily available with classical methods. It would have been “nice” if Euclidean time correlations functions could be readily obtained from Monte Carlo simulations and analytically continued to real time (see e.g., Ref. Meyer (2011)). As discussed in Sec. V.4, this procedure is in general ill-posed for finite precision correlators (e.g., those extracted from MC sampling), but nevertheless in some situations it is possible to gain considerable insight into transport properties by analyzing features of the Euclidean correlation functions (see e.g., Ref. Roggero and Reddy (2016)). However, a large-scale, high-fidelity quantum computer is expected to be able to evaluate real-time correlators accurately. Depending on the particular problem at hand, it could also be advantageous to work directly in frequency space using well-behaved integral kernels e.g., the Fejer Roggero and Carlson (2019) or Gaussian Roggero (2020). Similar in spirit are approaches to extract Green’s functions from a quantum simulation Ciavarella (2020), which uses the expectation values $F(t)=\langle\Psi|e^{-iHt}|\Psi\rangle$ evaluated over a range of times to approximate $g(\omega,\eta)=\langle\Psi|\frac{1}{\omega-H+i\eta}|\Psi\rangle$, which in the limit $\eta\to 0$ gives the exact Green’s function. This Green’s function algorithm was investigated through a simple model of the $\phi\to 2\chi$ decay of a (fictitious) heavy scalar $\phi$ into two lighter scalar particles $\chi$ using one of IBM’s superconducting quantum devices Ciavarella (2020), as shown in Figure 3. Variational approaches with reduced computational cost have also been proposed recently Endo et al. (2020), possibly opening the way to exploration of more realistic simulation on near term devices.

Another advantage of frequency space methods is their ability to provide access to information about the possible final states of a nuclear reaction, thus opening the way to a better characterization of semi-exclusive scattering cross sections. This capability will be important in situations where reliable theoretical predictions for the properties of final states is required as input to achieve the physics goals of an experiment. An example of this is long-baseline neutrino oscillation detectors, e.g., DUNE, which are designed to constrain neutrino properties like masses, mixing-angles and the CP-violating phase in the lepton sector. Even though these parameters are captured by the inclusive scattering of neutrinos off nuclei in the detector (followed possibly by a particle shower), the strong energy dependence of the inclusive cross section, together with the wide range of neutrino energies in the beam, requires an explicit reconstruction of the incoming neutrino energy on an event-by-event basis (using information from the observed final states). For Argon based detectors like DUNE, it is possible to reconstruct the trajectories of charged particles, and in principle reconstruct the kinematic parameters of the reaction. In practice, however, additional theoretical input is required in order to properly take into account the possible presence of neutral particles in the final state. The dominant quantum resource required for the quantum simulation of these processes is connected with the implementation of the energy filters used to approximately project into intervals of the excitation energy. A first estimation of the quantum resources required to perform a simulation with possible impact on the physics goals of DUNE was carried out recently Roggero et al. (2020a), using leading-order pionless EFT on a lattice as a minimal choice for realistic nuclear interactions. These results suggest that, at least for the Fejer kernel, a minimally realistic simulation of scattering from ⁴⁰Ar will require $\mathcal{O}(10^{8})$ one- and two-qubit operations. Using recent improvements in both integral kernels Roggero (2020) and more efficient representations of the real-time evolution operator for similar models Su et al. (2020), this estimate could be reduced below $\mathcal{O}(10^{6})$, leaving realistic studies of nuclear scattering out of reach for NISQ devices.

Thinking beyond the universal gate set towards custom-designed operators, SRF-cavity quantum devices, under development at Lawrence Livermore National Laboratory (LLNL), Fermi National Accelerator Laboratory (FNAL) and elsewhere, are superconducting quantum devices capable of operating high dimensional Hilbert spaces per cavity. Time-ordered radio-frequency signals are introduced into the cavities to accomplish unitary transformations between the many transmon-coupled states supported in high-Q cavities, describing the time evolution of the QMB system of interest. Classical simulations performed to support the LLNL cavity devices show that entangled spin systems, relevant to low-energy nuclear physics, can be efficiently evolved in time with high fidelity Holland et al. (2020). Figure 4 shows the classically-simulated frequency-space evolution of the exact neutron-neutron spin persistence probability during time evolution.

Beyond the preparation of a good approximation $|\Psi_{0}\rangle$ to the ground state of the target nucleus, techniques working directly in frequency space Roggero and Carlson (2019) require one more non-trivial step: applying the vertex operator $\hat{O}$ describing the interaction between probe and target while maintaining a normalized state $|\Phi_{O}\rangle\propto\hat{O}|\Psi_{0}\rangle$. The main complication in this step is that the vertex operator $\hat{O}$ is Hermitian but not necessarily unitary. A detailed comparison of two approaches to carry out this step has been made Roggero et al. (2020b) using a schematic representation of the $np\leftrightarrow d\gamma$ reactions. The first method Roggero and Carlson (2019) proposed using real time evolution under the action of the vertex operator $\hat{O}$ controlled by the state of an additional ancilla qubit, while the second method Childs and Wiebe (2012) is based on the general linear combination of unitaries (LCU) scheme, developed predominantly for fault tolerant quantum devices (as in general it requires control logic and large ancilla registers). Results obtained in Ref. Roggero et al. (2020b) and shown in Figure 4, indicate that the computational cost of LCU can be substantially reduced by a careful choice of operator basis and by optimizing the implementation with respect to the limited connectivity available in the device. The net result is that the LCU-based state preparation is expected to outperform the more memory efficient time-dependent algorithm in both the fidelity and success probability in most cases.

Neutrinos initially produced in a flavor eigenstate can experience flavor oscillations, with frequencies determined by its energy and the mass-squared gaps between mass eigenstates. The flavor evolution of a neutrino “cloud” can be modified substantially through weak interaction with a background e.g., the MSW effect generated by interactions with a background of leptons Wolfenstein (1978); Mikheyev and Smirnov (1985). In astrophysical settings with large neutrino densities e.g., core-collapse supernovae and the early universe, neutral-current neutrino-neutrino interactions can become important and modify substantially the transport of flavor in these environments. Neutrino-neutrino interactions can lead to collective neutrino oscillations, and these modes can have important effects on the dynamics of extreme astrophysical environments and the early universe Fuller et al. (1987); Savage et al. (1991). Considering only the leading order neutrino-neutrino interactions in forward scattering, and neglecting momentum-changing interactions suppressed as $G_{F}^{2}$, the dynamics act only on the flavor degrees of freedom, which can be represented by a collection of $SU(N_{f})$ spins. In this framework, the neutrino momenta are fixed and only the spin (flavor) polarization evolves in time. In the simplified approximation with $N_{f}=2$ flavors, the model Hamiltonian is equivalent to an all-to-all Heisenberg model with an interaction strength proportional to the neutrino density Pehlivan et al. (2011).

In general, solving for the exact real-time dynamics of these models is out of reach to direct methods based on exact diagonalization. In highly symmetric homogeneous situations, the neutrino Hamiltonian becomes integrable and the Bethe ansatz can be used Pehlivan et al. (2011); Patwardhan et al. (2019). For more general conditions, mean field approximations are usually employed (see e.g., Refs. Duan et al. (2010); Chakraborty et al. (2016) for recent reviews). Due to the widespread use of mean-field simulation techniques, understanding the accuracy of this approximation in capturing the dominant effects in collective oscillations is therefore essential Bell et al. (2003); Friedland and Lunardini (2003); Cervia et al. (2019); Rrapaj (2020). Due to the complexity of the calculations, however, the systems sizes that can be explored with these direct approaches is limited to $\mathcal{O}(10)$, making it difficult to understand extrapolations to the thermodynamic limit.

One venue where quantum technologies can help in breaking this computational barrier is in simulations of the flavor dynamics. First explorations of this possibility with superconducting qubit devices on small systems with up to $N_{\nu}=4$ neutrino amplitudes have been performed Hall et al. (2021); Yeter-Aydeniz et al. (2021). One of the main bottlenecks for efficient implementations of these models is the need to implement the all-to-all pairwise interactions on devices with limited connectivity. A general procedure to side-step this problem has been proposed Hall et al. (2021) by means of a swap-network construction that can be implemented efficiently even on devices with linear nearest-neighbor connectivity. The detectable difference from the mean field approximation may result from non-trivial evolution of entanglement and different techniques to reliably extract information about different entanglement measures from a digital quantum simulation have been presented in Ref. Hall et al. (2021). A complementary approach is to exploit the long-range nature of interactions in trapped-ion systems Monroe et al. (2021) to carry out analog simulations of flavor evolution, by generalizing earlier work on dynamical phase transitions in Ising models Zhang et al. (2017).

First explorations of employing an MPS ansatz for the many-body neutrino wave functions have been carried out Roggero (2021a, b). These calculations show that, at least for conditions with a high degree of symmetry, the bipartite entanglement in a neutrino cloud starting as a product state grows slowly with system size as $\approx\log(N_{\nu})$. This allows for controllable simulations of more than a hundred neutrino amplitudes and the identification of a link between the appearance of instabilities in the flavor evolution and the presence of a dynamical phase transition in the underlying spin model Roggero (2021b). This connection has the potential of providing a general framework to understand the conditions leading to collective oscillation phenomena without resorting to a mean-field approximation. A particularly exciting prospect of using MPS-based techniques for studies on collective neutrino oscillations is the ability to detect the presence of configurations capable of generating substantial levels of entanglement by checking the convergence of the calculations. The identification of these special conditions, if they exist, will allow a focusing of quantum resources on the simulation of these instances, which will be out of the reach of tensor-network-based methods.