Quantum computers as universal quantum simulators: state-of-art and perspectives

TL;DR

This paper reviews the current state and future prospects of using near-term quantum computers, based on trapped ions and superconducting qubits, as universal quantum simulators for complex physical models.

Contribution

It provides a comprehensive overview of recent experimental achievements, theoretical background, and future challenges in digital quantum simulation with noisy intermediate-scale quantum devices.

Findings

Quantum devices can simulate spin models with tens of qubits.

Comparison of trapped ion and superconducting architectures.

Outlook on achieving quantum advantage in many-body simulations.

Abstract

The past few years have witnessed the concrete and fast spreading of quantum technologies for practical computation and simulation. In particular, quantum computing platforms based on either trapped ions or superconducting qubits have become available for simulations and benchmarking, with up to few tens of qubits that can be reliably initialized, controlled, and measured. The present review aims at giving a comprehensive outlook on the state of art capabilities offered from these near-term noisy devices as universal quantum simulators, i.e. programmable quantum computers potentially able to calculate the time evolution of many physical models. First, we give a pedagogic overview on the basic theoretical background pertaining digital quantum simulations, with a focus on hardware-dependent mapping of spin-type Hamiltonians into the corresponding quantum circuit model as a key initial step towards simulating more complex models. Then, we review the main experimental achievements obtained in the last decade regarding the digital quantum simulation of such spin models, mostly employing the two leading quantum architectures. We compare their performances and outline future challenges, also in view of prospective hybrid technologies, towards the ultimate goal of reaching the long sought quantum advantage for the simulation of complex many body models in the physical sciences.

Equations

$$\frac{d\vec{x}}{dt}=M\vec{x}$$

$$\vec{x}(t)=e^{Mt}\vec{x}(0)$$

$$\frac{d|\Psi ⟩}{dt}=-i\mathcal{H}|\Psi ⟩$$

$$|\Psi (t)⟩=\mathrm{U}(t)|\Psi (0)⟩$$

$$\mathcal{H}=\sum _l^L{\mathcal{H}}_l$$

$$\tilde{\mathrm{U}}=\prod _l{\mathrm{U}}_l(t)$$

$$e^{-i{\sum }_l{\mathcal{H}}_{l}t}=\underset{n\to \infty }{\lim }{\left(\prod _l{e}^{-i{\mathcal{H}}_{l}t/n}\right)}^n$$

$$\mathrm{U}(t)=e^{-i{\sum }_l{\mathcal{H}}_{l}t}={\left(\prod _l{e}^{-i{\mathcal{H}}_{l}t/n}\right)}^n+O\left(\frac{t^2}{n}\right)$$

$$\mathcal{H}\to H(\{{\sigma }_{\alpha }\})$$

$$H=\sum _l{H}_l$$

$$e^{-iHt}=\prod _l{e}^{-i{H}_{l}t}$$

$$e^{-iHt}\simeq {\left(\prod _l{e}^{-i{H}_{l}t/n}\right)}^n$$

$$\begin{array}{r}{\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right){\sigma }_y=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right){\sigma }_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\end{array}$$

$$[{\sigma }_{\alpha },{\sigma }_{\beta }]=2i{ϵ}_{\alpha \beta \gamma }{\sigma }_{\gamma },\{{\sigma }_{\alpha },{\sigma }_{\beta }\}=2{\delta }_{\alpha \beta }\mathbb{I}$$

$$H=\sum _{\begin{array}{c}i=1\\ \alpha =x,y,z\end{array}}^N{h}_{\alpha ,i}^{(1)}{\sigma }_{\alpha }^{(i)}+\sum _{\begin{array}{c}i,j=1\\ \alpha ,\beta =x,y,z\end{array}}^N{h}_{\alpha \beta ,ij}^{(2)}{\sigma }_{\alpha }^{(i)}{\sigma }_{\beta }^{(j)}$$

$$\mathrm{U}(\theta ,\varphi ,\lambda )=\left(\begin{array}{cc}\cos (\theta /2)& -e^{i\lambda }\sin (\theta /2)\\ e^{i\varphi }\sin (\theta /2)& e^{i(\lambda +\varphi )}\cos (\theta /2)\end{array}\right)$$

$$\mathrm{H}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)$$

$$\Phi (\delta )=\left(\begin{array}{cc}1& 0\\ 0& e^{i\delta }\end{array}\right)$$

$$\mathrm{U}(\theta ,\varphi ,\lambda )=e^{-i\theta /2}\Phi \left(\frac{\pi }{2}+\varphi \right)\mathrm{H}\Phi (\theta )\mathrm{H}\Phi \left(-\frac{\pi }{2}+\lambda \right)$$

$${\mathrm{R}}_{\alpha }(\theta )=\exp \left(-i\frac{\theta }{2}{\sigma }_{\alpha }\right)\alpha =x,y,z$$

$$\mathrm{U}(\theta ,\varphi ,\lambda )={\mathrm{R}}_z(\varphi ){\mathrm{R}}_x(\theta ){\mathrm{R}}_z(\lambda )$$

$$H_1^{(i)}=\sum _{\alpha =x,y,z}{h}_{\alpha ,i}^{(1)}{\sigma }_{\alpha }^{(i)}$$

$${\mathrm{U}}_1^{(i)}(t)=e^{-i{H}_1^{(i)}t}$$

$${\mathrm{U}}_{\alpha \beta }^{(i,j)}(t)=e^{-i{H}_{\alpha \beta }^{(i,j)}t}=e^{-i\delta {\sigma }_{\alpha }^{(i)}\otimes {\sigma }_{\beta }^{(j)}}$$

$$H=J\sum _{⟨i,j⟩}\left({\sigma }_x^{(i)}{\sigma }_x^{(j)}+{\sigma }_y^{(i)}{\sigma }_y^{(j)}+{\sigma }_z^{(i)}{\sigma }_z^{(j)}\right),$$

$$H=\sum _{⟨i,j⟩}\left(J_{xx}{\sigma }_x^{(i)}{\sigma }_x^{(j)}+J_{yy}{\sigma }_y^{(i)}{\sigma }_y^{(j)}+J_{zz}{\sigma }_z^{(i)}{\sigma }_z^{(j)}\right),$$

$$H=\sum _i{h}_i{\sigma }_x^{(i)}+\sum _{⟨i,j⟩}{J}_{zz}{\sigma }_z^{(i)}{\sigma }_z^{(j)}$$

$${\mathcal{S}}_1=\{{\mathrm{R}}_{\alpha }(\theta ),\mathrm{CNOT}\}$$

$$\mathrm{CNOT}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)$$

$$\mathrm{ZZ}(\delta )=e^{-i\delta {\sigma }_z\otimes {\sigma }_z}$$

Full text

I Introduction

When trying to accurately describe the dynamical behavior of physical systems made of several interacting fundamental constituents, and from these explain the complexity of natural aggregates following a bottom up approach, the well established classical laws of physics fail to give an accurate picture of reality, as it is now accepted and understood. In fact, quantum mechanics is arguably the most complete and successful theory we currently have to effectively describe the dynamics of the elementary constituents of our universe. A great deal of methods and simulation tools have been developed in the last century, such as quantum Monte-Carlo gubernatis_quantum_2016 , molecular dynamics haile_molecular_1992 , and tensor networks montangero_introduction_2018 to name a few examples, which allow solving some of the theoretical models formulated in quantum mechanical terms and correctly describe a large variety of quantum phenomena. The very concept of “simulation” has a broadly understood meaning in Science, Technology, Engineering, and Mathematics (STEM) applications. In fact, simulating any natural phenomenon is equivalent to artificially reproduce its properties and its dynamical evolution in time. This is primarily carried out through an accurate mathematical modeling, i.e. a mapping of the information we know about a system of interest onto a certain set of variables and equations, followed by an analytic or most often numerical solution. The resulting set of mathematical identities (or the computer with its numerical program aimed at solving them) can then be named a simulator. Such a simulator is used to study the behavior of the real system under fairly general conditions, to make predictions and to test new hypotheses, the only limitations being the validity of the initial modeling and the available computational power.

It is generally accepted that most of the models we currently deal with cannot be solved exactly with classical computing machines, such as modern supercomputers obeying the laws of classical physics. The main reason lies in the exponential scaling of time and memory resources needed to correctly capture the dynamics of the relevant physical variables with increasing system size. This is especially true when strong correlations between the system parties play a dominant role, which is the case in most interesting situations. In such cases, even the most elaborate but inevitably approximate classical simulation approaches so far developed fail in giving the correct answers. Hence, quantum simulators have long been proposed as a possible solution, building on the general idea that since Nature ultimately behaves quantum mechanically, only a computing machine obeying quantum mechanical laws would be able to accurately simulate it manin_computable_1980 ; benioff_computer_1980 ; feynman_simulating_1982 . A quantum simulator is a system under high control of the experimenter, which is able to reproduce the dynamical behavior of a given physical model, irrespective of the degree of internal correlations or entanglement between the model’s degrees of freedom. Following this route, a plethora of analog quantum simulators have been proposed fisher_boson_1989 ; jaksch_cold_1998 ; hartmann_strongly_2006 ; greentree_quantum_2006 ; angelakis_photon-blockade-induced_2007 ; carusotto_numerical_2008 ; carusotto_fermionized_2009 ; tian_circuit_2010 ; gerace_analog_2012 ; viehmann_observing_2013 ; du_superconducting_2015 ; reiner_emulating_2016 and developed greiner_quantum_2002 ; friedenauer_simulating_2008 ; gerritsma_quantum_2010 ; kim_quantum_2010 ; islam_onset_2011 ; nguyen_acoustic_2015 ; steinhauer_observation_2016 ; labuhn_tunable_2016 ; bernien_probing_2017 ; roushan_spectroscopic_2017 ; zhang_observation_2017 ; zhang_experimental_2018 , in which the physical properties of a targeted model are reproduced on a physical set-up under externally controlled conditions. On the other hand, digital quantum simulators are programmable and general purpose quantum devices, which promise a larger flexibility on the models to be solved lloyd_universal_1996 ; ortiz_quantum_2001 ; somma_simulating_2002 ; verstraete_quantum_2009 ; weimer_rydberg_2010 ; santini_molecular_2011 ; casanova_quantum_2012 ; mezzacapo_digital_2012 ; jordan_quantum_2012 ; raeisi_quantum-circuit_2012 ; hauke_quantum_2013 ; las_heras_digital_2014 ; chiesa_digital_2015 ; garcia-alvarez_digital_2017 ; jiang_quantum_2018 ; kivlichan_quantum_2018 . In this respect, digital quantum simulators are quantum computing machines not restricted to emulate the dynamics of targeted models, but satisfying DiVincenzo criteria divincenzo_physical_2000 for quantum computation. Here we will consider such digital quantum computers as universal quantum simulators (UQS) lloyd_universal_1996 , meaning that they are able, in principle, to reproduce with arbitrary precision the dynamics of any Hamiltonian model that can be suitably encoded on a given quantum register and translated into a sequence of gate operations, as schematically illustrated in Fig. 1. The time evolution of the physical model is mapped onto an effective model defined on the quantum hardware degrees of freedom, in which the time evolution can be programmed in digital steps through a sequence of unitary operations defined by a quantum circuit nielsen_quantum_2000 . This mapping will be the focus of the present review. As it will be recalled, the concept of universality requires that the model Hamiltonian be the sum of locally interacting terms, which in turn implies that the number of required operations does not grow exponentially with the system size. This, in practice, restricts the class of Hamiltonian models that can be usefully simulated (although to the ones that are most physically relevant).

On a more refined level, it is worth mentioning that hybrid digital-analog quantum simulators have also been proposed, aimed at combining the easier scalability of analog approaches with the intrinsic universality of digital quantum simulations mezzacapo_digital_2015 . Here, analog blocks allow for the direct simulation of the time dynamics on a large number of variables, thus reducing the number of digital operations and errors, while digital blocks are included to introduce a variety of possible interaction models. This paradigm is hailed as a promising route leading to universal digital-analog quantum computation.

Several excellent reviews have been published in the last few years, giving a broad account of quantum simulators, either general purpose buluta_quantum_2009 ; georgescu_quantum_2014 ; sanders_efficient_2013 ; preskill_quantum_2018 or more focused on specific categories and/or quantum hardware blatt_quantum_2012 ; bloch_quantum_2012 ; aspuru-guzik_photonic_2012 ; houck_-chip_2012 ; wendin_quantum_2017 ; lamata_digital-analog_2018 . Here we give a targeted overview on near term digital quantum computers as devices able to perform universal quantum simulations. This goes in line with the fast pace of advancement in different quantum computing technologies that have made programmable devices available, thus attracting widespread interest worldwide. In fact, current quantum processors promise to overcome the intrinsic limitations of simulating complex many body physics with classical computing machines, although it is still difficult to predict when this will happen. A targeted goal would be to reach the long-sought “quantum advantage”, i.e. a certified gain in either memory or temporal efficiency obtained for the solution of a quantum problem with respect to the equivalent simulation being performed on a classical supercomputer. Without entering into the subtleties related to a rigorous definition of quantum advantage, here we just consider that a quantum computer with fully operational $N=50$ qubits is able to store something like $8\times 2^{N}\sim 9\cdot 10^{15}$ bytes of information (i.e. 9 Pb, assuming 8 bytes to store a complex number in single-precision), which roughly corresponds to the random access memory of state-of-art supercomputers pednault_breaking_2017 ; boixo_characterizing_2018 . Should this threshold be met with an actual quantum simulation involving the whole quantum hardware in the current Noisy Intermediate Scale Quantum devices (NISQ) era preskill_quantum_2018 , it would represent the quantum advantage turning point. On a longer and still unpredictable timescale, even farther-reaching consequences are expected should fully fault tolerant and scalable quantum hardware become available schindler_experimental_2011 ; you_simulating_2013 ; barends_superconducting_2014 ; corcoles_demonstration_2015 , in which $N>100$ logical qubits have to be complemented with a much larger number of auxiliary quantum bits aimed at correcting errors.

Our aim is to give an overview of the field that could be useful to the beginning researcher or student, trying to keep a pedagogic approach over the elementary theoretical background throughout the manuscript, and then summarizing the main experimental achievements and prospective developments. Since most Hamiltonian models can be mapped onto spin-type ones, being able to efficiently simulate spin models on actual quantum computing devices is crucial, not only because they possess interesting many body dynamics themselves but also to open the door to the universal quantum simulation of a large class of quantum models (typically interacting fermionic particles) that are intractable by classical computation means troyer_computational_2005 . Paradigmatic examples are the Hubbard model in condensed matter casanova_quantum_2012 ; barends_superconducting_2014 , or the Schwinger model in lattice gauge field theory hauke_quantum_2013 ; martinez_real-time_2016 . In particular, we emphasize the role of specific quantities that are known to be difficult to compute but extremely important in the description of the dynamical properties of many body systems, such as quantum correlations.

In terms of actual quantum hardware, we will focus on reviewing the main experimental achievements obtained in the last decade, specifically dealing with the simulation accuracy of targeted spin Hamiltonians on different quantum platforms. While several alternatives are currently being pursued to realize actual non-error corrected quantum processors ladd_quantum_2010 , from photonic integrated circuits aspuru-guzik_photonic_2012 to spins in semiconductors awschalom_quantum_2013 , we concentrate upon the two leading architectures that have been dominating the scene: trapped ions manipulated through external microwave or optical fields monroe_scaling_2013 ; schindler_quantum_2013 ; bruzewicz_review-trap-ion_2019 , and superconducting circuits working at microwave frequencies clarke_superconducting_2008 ; schoelkopf_wiring_2008 ; devoret_superconducting_2013 ; Gu_review_SCcircuits_2017 . We anticipate that interesting results might already be within reach in NISQ processors, despite the relatively small number of useful operations and non-error corrected qubits currently available on such devices. On a parallel sight, while the main object is restricted to quantum simulations of physical models and STEM applications in general, actual quantum processors might eventually turn to solve complex problems in other fields as well. As examples, classification and scheduling tasks, stock market pricing woerner_quantum_2019 ; martin_towards_2019 and machine learning biamonte_quantum_2017 might benefit from speedup advantages over classical computers. The basics of quantum circuit programming reported in this review may be a useful starting point. Last but not least, these topics settle within the quantum technologies roadmap promoted at the European level through the recently funded Quantum Flagship acin_quantum_2018 .

II Theory of Digital Quantum Simulations

When the main object of a physical theory is to determine the evolution in time of a given system, most problems are formulated in terms of a set of differential equations. Their solution is at the heart of many simulation protocols nowadays, from molecular dynamics to aircraft design. A very common situation is, for example, a linear set of equations such as

[TABLE]

where $M$ is a matrix and $\vec{x}$ represents a vector of dynamical variables. Once an initial condition $\vec{x}(0)$ is given, the formal solution to the above equation is simply

[TABLE]

Implementing such a solution on a computer routine gives a useful tool to fully solve the system dynamics, provided that the size of the numerical problem is within reach of the available computational resources. In quantum mechanics, the paradigmatic example is the Schrödinger equation (here and in the following, we take $\hbar=1$)

[TABLE]

where $\mathcal{H}$ is known as the Hamiltonian operator. This complex-valued differential equation is solved by computing the unitary time-evolution operator $\mathrm{U}(t)=e^{-i\mathcal{H}t}$. Indeed, once the latter is known, any initial condition can be evolved linearly as

[TABLE]

Matrix exponentiation is a very common numerical task arising in many interesting simulation scenarios, and crucially in the field of quantum mechanical systems. On classical computers, this task turns out to be provably difficult in terms of the matrix size, most notably for quantum mechanical simulations, where the exponential increase of the size of the Hilbert space of a composite system with the number of sub-systems leads to an exponential demand of time and memory resources.

In 1982, Richard Feynman conjectured that using a controllable quantum mechanical system as a computing resource, instead of a classical object, would provide significant advantages in the simulation of quantum systems feynman_simulating_1982 . Indeed, just about fifteen years later, in 1996, Seth Lloyd proved that idea to be essentially correct lloyd_universal_1996 , with the sole limitation that the systems to be simulated only carry local interactions between their constituent subsystems. In these review, we will thus concentrate on system Hamiltonians of the form

[TABLE]

where $\mathcal{H}_{l}$ acts locally only on a portion of the total system. As a matter of fact, many Hamiltonian models of physical interest can be formulated as in Eq. (5).

II.1 The quantum computer as a universal quantum simulator

Given a Hamiltonian $\mathcal{H}$ that models the physical system under investigation, the problem of computing the corresponding time evolution operator $\mathrm{U}(t)=\exp\left(-i\mathcal{H}t\right)$ is equivalent to the task of implementing a well defined unitary matrix. A quantum computer endowed with a universal set of quantum gates is by definition able to perform any unitary transformation, albeit not necessarily in an efficient number of elementary operations nielsen_quantum_2000 . What Lloyd actually proved is that universal quantum computers can calculate $\mathrm{U}(t)$ efficiently (i.e. with polynomial time and memory resources in the size of the target system) when $\mathcal{H}$ is a sum of local terms. The proof is based on two fundamental facts: first, in the circuit model for UQS we can implement generic transformations by successively performing elementary unitary operations (quantum gates), and appending one unitary $\mathrm{U}_{A}$ after another $\mathrm{U}_{B}$ in the circuit results in a total unitary, which is mathematically the product $\mathrm{U}_{A}\mathrm{U}_{B}$ being applied to the state of the qubit register. Second, any unitary operation $\mathrm{U}$ acting on $N$ qubits can be implemented with $O(2^{2N})$ elementary operations (we recall that the dimension of the Hilbert space associated to $N$ qubits is $d=2^{N}$) barenco_elementary_1995 ; nielsen_quantum_2000 . Following Ref. lloyd_universal_1996, , we now suppose that we are given a Hamiltonian expressed as the sum of $L$ local terms, as in Eq. (5), where $L\propto p\cdot N$ such that $p$ measures some degree of locality (e.g., the number of nearest neighbors or second-to-nearest neighbors in a lattice), and $N$ is the total number of qubits required to encode the computation. Hence, the number of local terms, $L$, scales polynomial with $N$. In general, according to the rules above, computing directly $\mathrm{U}(t)=\exp\left(-i\mathcal{H}t\right)$ requires $O(2^{2N})$ operations, and it is therefore exponentially inefficient. However, let us call $m_{l}$ the dimension of the subsystem over which the action of $\mathcal{H}_{l}$ is restricted. Typically, we will have $m_{l}\ll 2^{N}$, since local terms only involve few-body interactions. In this case, the unitary $\mathrm{U}_{l}(t)=\exp\left(-i\mathcal{H}_{l}t\right)$ can be computed with $O(m_{l}^{2})$ operations. The overall product

[TABLE]

can therefore be obtained on a universal quantum computer by juxtaposing the circuit implementations of the single $\mathrm{U}_{l}(t)$ unitaries, and it would take at most $O(Lm_{max}^{2})$ elementary operations, where $m_{max}=\operatorname{max}_{l}{m_{l}}$. The final step of the reasoning lies in the following mathematical identity, which is known as the Suzuki-Trotter (ST) decomposition:

[TABLE]

Unless all the $\mathcal{H}_{l}$ operators commute, in which case the ST identity is exact already for $n=1$, the product of local unitaries will not be exactly equal to the total target unitary $\mathrm{U}(t)=\exp\left(-i\mathcal{H}t\right)$. However, it can be shown that $\forall n$

[TABLE]

which means that we can approximate arbitrarily well the desired unitary operator by repeating $n$ times the sequence of gates corresponding to the product of local terms for time slices $t/n$. All in all, we were able to break our original problem into smaller pieces, $e^{-i\mathcal{H}_{l}t/n}$, which can now be implemented efficiently using only a limited set of elementary gates and give the correct answer up to an arbitrarily small digital error $O(t^{2}/n)$. Indeed, for any $\epsilon>0$ and $t$, there exists a $n_{\epsilon}$ such that $\mathrm{U}(t)$ can be computed within an approximation $\epsilon$ in at most $n_{\epsilon}Lm_{max}^{2}$ operations. This is polynomial in $N$ whenever $L=\text{poly}(N)$, as for example in the case of nearest neighbors interactions.

II.2 Quantum simulations cookbook

From now on, we will assume to work with a universal quantum computer, described in the standard circuit model as a (quantum) digital device, i.e. qubit-based, obeying the algebra of Pauli matrices and operating with a universal set of quantum gates nielsen_quantum_2000 . The problem of quantum simulation can then be formulated and solved on such a machine by taking a few simple steps, which we are going to outline in the following.

First, define a model Hamiltonian of interest $\mathcal{H}$. This should contain all the dynamical information necessary to describe and characterize the physical quantum system under investigation. The most appropriate set of variables and operators will appear in the mathematical structure of $\mathcal{H}$.

Second, map the target Hamiltonian $\mathcal{H}$ onto its representation on the qubit Pauli algebra

[TABLE]

In simpler terms, this means finding a suitable encoding of the degrees of freedom of the target system into a number $N$ of qubits. The resulting mapped Hamiltonian $H$ will then be written in terms of Pauli matrices. Notice that this mapping is straightforward for physical systems consisting of collections of spin-$1/2$ objects, as they also obey Pauli algebra, but it is possible in principle for a large class of physical system, as it will be shown in the following with some specific example. The quantum simulation will be efficient whenever such $H$ is the sum of local terms. Notice that this is usually not a limitation in many practical cases, as most physical processes are inherently local in nature. However, local Hamiltonian models get mapped into non-local ones, as, e.g., when the well known Jordan-Wigner transformation is applied jordan_uber_1928 to encode fermionic degrees of freedom. These models might be efficiently simulated on quantum hardware implementing multi-qubit gates, as it will be detailed in the following.

Third, assuming the target Hamiltonian is mapped onto a sum of local contributions

[TABLE]

check whether $[H_{l},H_{l^{\prime}}]=0\,\,\forall l,l^{\prime}$. If that is the case, then

[TABLE]

with no digital error. Otherwise, choose the number of ST steps (sometimes referred to as Trotter steps), $n$, that is appropriate for the required degree of precision, in such a way that

[TABLE]

This application of the ST formula is sometimes called Trotterization in quantum simulations jargon.

Fourth, translate each local unitary $e^{-iH_{l}t}$ (or $e^{-iH_{l}t/n}$) into a sequence of quantum gates. This is always possible in at most $O(m_{l}^{2})$ operations and with any universal set of single- and two-qubits operations available on a general purpose quantum computer nielsen_quantum_2000 . The total quantum circuit encoding the time evolution will be the juxtaposition of all the sequences corresponding to the factors in the ST decomposition, repeated $n$ times.

Finally, add initial state preparation at the beginning of the circuit and an appropriate set of measurements at the end to recover expectation values of the relevant observable quantities on the evolved quantum state.

The points above represent a quite general set of instructions towards the design of a quantum simulation algorithm. In the following, we will give some explicit examples to show how this is done in practical cases. Of course, such techniques are not limited to actual simulations of real physical systems, but can become a tool for a larger class of computational tasks whenever the problem of interest can be encoded in a Hamiltonian quantum dynamics.

III Quantum Circuits

Among the steps that must be undertaken in order to practically design and realize a digital quantum simulation, the translation of unitary operators into elementary quantum gates is the one that is most typically hardware-dependent. It is also critical in terms of results and performance, particularly in the present era of noisy and intermediate-scale prototypes of quantum processors, where the interplay between hardware properties and target features is stronger.

Several universal sets of single- and two-qubit gates are known nielsen_quantum_2000 , all in principle equally valid as a primitive set to realize any quantum simulation. However, every real hardware platform usually comes with a native set of operations that, due to the physical characteristics of the device, are readily implemented in practice. The platform is in itself capable of implementing universal quantum computation, and is thus a potential UQS, if and only if the native set is a universal set in the usual quantum computing sense. If that is the case, any target unitary evolution can be translated in a combination of the native operations without unnecessary overhead. Processors based on different technological platforms may also exhibit distinct topological properties, i.e. different qubit-qubit inter-connectivity and limitations in gate directionality. While these do not pose hard limitations to the computational power of the platform, since they can always be compensated via, e.g., $\mathrm{SWAP}$ operations, they may results in some overhead in the total length of the simulations. Hence, in this NISQ era some platforms are more suitable for the simulation of certain physical models (e.g. trapped ions, featuring built-in all-to-all connectivity, can more easily simulate long-range interactions), thus making a fair comparison of performances less straightforward cross_validating_2018 . Of course, it should be reminded that, as a general rule, only systems described by local interaction terms are somehow guaranteed to be efficiently mapped on a quantum computing register.

III.1 Pauli algebra and spin Hamiltonians

The mathematical properties of qubits are those of spin-$1/2$ systems, thus obeying the algebraic properties of Pauli matrices. The latter can be written in the computational basis representation as

[TABLE]

and satisfy the following commutation and anti-commutation rules

[TABLE]

where $\alpha,\beta,\gamma\in\{x,y,z\}$, $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor, $\delta_{\alpha\beta}$ is the Kronecker delta and $\mathbb{I}$ is the identity matrix.

In order to be simulated on a qubit-based architecture, any target Hamiltonian, $\mathcal{H}$, has to be mapped into an equivalent Hamiltonian, $H$, of interacting spin-$1/2$ operators. As already mentioned, this step is straightforward for paradigmatic spin-$1/2$ Hamiltonians (e.g., implementing Heisenberg or Ising models), but effective mappings are known for a large variety of cases, ranging from spin $S>1/2$ santini_molecular_2011 ; chiesa_digital_2015 ; tacchino_electromechanical_2018 to fermionic and fermionic-bosonic systems ortiz_quantum_2001 ; santini_molecular_2011 ; casanova_quantum_2012 ; mezzacapo_digital_2012 ; chiesa_digital_2015 ; mezzacapo_digital_2015 ; barends_digital_2015 ; garcia-alvarez_digital_2017 ; kandala_hardware-efficient_2017 ; dumitrescu_cloud_2018 ; kivlichan_quantum_2018 , including lattice models related to gauge theories martinez_real-time_2016 ; klco_quantum-classical_2018 . The generator of time evolution in a $N$-qubit digital quantum simulation therefore takes the general form

[TABLE]

containing in general both single- and two-spin terms, to which any other many body term time evolution can, in principle, be reduced (see Sec. III.4). Whenever the overall structure of $H$ retains a local nature, as it is the case for many physically relevant examples, its translation into elementary gate operations can be done efficiently. In the following, we will provide a dictionary of useful decomposition rules in terms of different universal sets of gates. Most of them are derived from real use-case scenarios and can therefore be straightforwardly applied to well known physical models.

III.2 Single-qubit rotations

Once the target Hamiltonian $\mathcal{H}$ is reduced to its counterpart $H$ on $N$ spin-$1/2$ systems, a register of $N$ qubits can be used to encode and carry out the quantum simulation via the identification of each qubit with a single spin-$1/2$ element. All currently proposed and realized quantum computing platforms allow addressing single qubits with tailored control pulses to perform single qubit gates. The most general single qubit $\mathrm{SU}(2)$ operation has the form

[TABLE]

and can be obtained, for example, by combining well known single qubit quantum gates such as the Hadamard gate

[TABLE]

and the phase gate

[TABLE]

Indeed, the following identity holds:

[TABLE]

Rotations around the coordinate axes

[TABLE]

can be implemented, up to global phase factors, by choosing particular parameters in $\mathrm{U}(\theta,\phi,\lambda)$. For example, $\mathrm{R}_{z}(\lambda)=e^{-i\lambda/2}\Phi(\lambda)=\mathrm{U}(0,0,\lambda)$, $\mathrm{R}_{x}(\theta)=\mathrm{U}(\theta,-\pi/2,\pi/2)$ and $\mathrm{R}_{y}(\theta)=\mathrm{U}(\theta,0,0)$. Vice-versa, any platform capable of implementing single-qubit rotations around the coordinate axes can in principle realize an arbitrary $\mathrm{U}(\theta,\phi,\lambda)$ via the following identity

[TABLE]

In Eq. (15), any single-spin term

[TABLE]

essentially represents a magnetic field applied to the $i$-th qubit along the direction identified by the vector $\vec{h}=(h^{(1)}_{x},h^{(1)}_{y},h^{(1)}_{z})$. The induced time evolution

[TABLE]

is a precession around the $\vec{h}$ axis, with the corresponding action on a qubit being a rotation of the Bloch vector. This can always be expressed in the $\mathrm{U}(\theta,\phi,\lambda)$ form, and therefore as a combination of rotations around the coordinate axes or of Hadamard and phase gates. Other decompositions of general $\mathrm{SU}(2)$ transformations, as well as approximate results employing only a finite set of fixed-phase single qubit operations instead of continuous-valued ones, are also known. nielsen_quantum_2000 ; chow_universal_2012

III.3 Two-qubits gates

Two-spin interactions appearing in the general Pauli Hamiltonian, Eq. (15), are usually implemented in digital quantum simulation protocols as combinations of single- and two-qubits gates. The typical evolution operator has the form

[TABLE]

where $\delta$ is a dimensionless phase factor. These terms arise naturally in the simulation of many renown spin models such as the Heisenberg model

[TABLE]

the $\mathrm{XYZ}$ model

[TABLE]

which reduces to the so called $\mathrm{XY}$ model if $J_{zz}=0$, or the transverse field Ising model

[TABLE]

Here $\langle i,j\rangle$ denote nearest neighbors spin pairs.

The exact and most effective decomposition of $\mathrm{U}_{\alpha\beta}^{(i,j)}(t)$ terms into elementary quantum gates varies from platform to platform, depending on the available set of native operations. One common situation, typical of, e.g., superconducting qubit technology with cross-resonance interactions rigetti_fully_2010 ; chow_simple_2011 ; sheldon_procedure_2016 , is a native universal set

[TABLE]

containing single qubit rotations and the two-qubit $\mathrm{CNOT}$ entangling gate

[TABLE]

Let $\mathrm{ZZ}(\delta)$ be the unitary operation

[TABLE]

This can be realized using the elementary quantum gates belonging to $\mathcal{S}_{1}$ with the following quantum circuit:

[TABLE]

Other terms generated by $\sigma_{\alpha}\otimes\sigma_{\beta}$ can be obtained from the construction above by suitable changes of reference frames, implemented with single qubit rotations. Indeed, remembering the following identities

[TABLE]

it is straightforward to verify that

[TABLE]

and

[TABLE]

These gate sequences can be combined to simulate all of the paradigmatic spin models mentioned above. For example, it is straightforward to prove that for the two-qubit Heisenberg model we have

[TABLE]

where $\mathrm{AB}(\delta)=e^{-i\delta\sigma_{a}\otimes\sigma_{b}}$, and all of the terms on the right hand side commute with each other. More detailed examples will be given in Sec. III.7.

Another universal set, defined $\mathcal{S}_{2}$, that often arises in superconducting realizations and proposals of quantum simulators replaces the $\mathrm{CNOT}$ gate with a parametric $\mathrm{XX}+\mathrm{YY}$ interaction las_heras_digital_2014 ; salathe_digital_2015 ; mckay_universal_2016 ; tacchino_electromechanical_2018

[TABLE]

In this case, we can take as the fundamental building block $\mathrm{XX}(\delta)=e^{-i\delta\sigma_{x}\otimes\sigma_{x}}$, to which all other unitary evolution terms generated by $\sigma_{\alpha}\otimes\sigma_{\beta}$ can be reduced with single-qubit changes of reference frame. The $\mathrm{XX}(\delta)$ gate is realized in $\mathcal{S}_{2}$ as

[TABLE]

Finally, let us call $\mathcal{S}_{3}=\{\mathrm{R}_{\alpha}(\theta),\mathrm{C}\Phi(\delta)\}$ the universal set of quantum gates containing all single qubit rotations and the controlled phase gate

[TABLE]

The latter is natively implemented on superconducting platforms with state dependent frequency shifts carretta_quantum_2013 ; barends_superconducting_2014 ; chiesa_robustness_2014 ; barends_digital_2015 ; reagor_demonstration_2018 , and is closely related to the Ising interaction generated by $H_{\text{Ising}}\propto\sigma_{z}\otimes\sigma_{z}$ jones_robust_2003 . In view of the latter property, it is not surprising that the $\mathrm{ZZ}(\delta)$ building block can be obtained directly from a single $\mathrm{C}\Phi(\delta)$ just with single qubit corrections and apart from an overall phase:

[TABLE]

An equivalent construction with two $\mathrm{C}\Phi(\delta)$ is the following

[TABLE]

where rotations around $\alpha=x,y$ enable the range of negative and small angles in those real experimental setups where the achievable phases $\delta$ in a single $\mathrm{C}\Phi(\delta)$ gate might be limited due to hardware constraints barends_digital_2015 .

In quantum simulators based on trapped ions technology schindler_quantum_2013 ; schindler_quantum_2013-1 ; lanyon_universal_2011 , the fundamental set of operations, which we will call $\mathcal{S}_{4}$, typically includes individual single qubit $z$ rotations

[TABLE]

collective non-entangling operations

[TABLE]

where $\sigma_{\phi}=\cos\phi\sigma_{x}+\sin\phi\sigma_{y}$, and Mølmer-Sørensen collective entangling gates molmer_multiparticle_1999

[TABLE]

Any subset of qubits can in principle be addressed with the collective gates, while leaving the others untouched. On a 2-qubit quantum register, $\mathrm{T}_{4}(\delta,0)$ can for example be used to obtain $\mathrm{XX}(\delta)$. Of course, the naturally collective character of trapped ions quantum gates is best exploited for the quantum simulation of long range and multiple-body interactions.

It is worth pointing out that while the elementary decomposition of typical two-qubits interaction terms reported here can be used to perform the digital quantum simulation of generic spin Hamiltonians, this is not necessarily the optimal strategy in general. Indeed, further optimization of, e.g., combined two qubit operations can lead to an overall reduction of the total number of gates for particular target Hamiltonian models vidal_universal_2004 ; las_heras_digital_2014 ; klco_quantum-classical_2018 ; chiesa_quantum_2019 . Examples of these techniques applied to the Heisenberg model simulated with $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ universal sets are discussed in Sec. III.7 below.

III.4 Multiple-qubit interactions

The generalization of $\mathrm{U}_{\alpha\beta}^{(i,j)}(\delta)$ building blocks to $N$-qubit interactions leads to unitary evolution terms of the form

[TABLE]

These can be in principle always decomposed into single- and two-qubit operations. An example within the $\mathcal{S}_{1}$ universal set is the following:

[TABLE]

The pattern can be generalized to any $N>3$, and changes of reference frames can be applied to individual qubits as done for the $N=2$ case.

In trapped ions processors, whose universal set $\mathcal{S}_{4}$ natively contains many body interactions, the decomposition of $N$-body terms can usually be done very efficiently using Mølmer-Sørensen gates muller_simulating_2011 and the limits on $N$ are in principle dictated only by the scalability of the hardware set-up itself.

III.5 Suzuki-Trotter decomposition and digital error

When designing a quantum simulation which requires the non-trivial application of Suzuki-Trotter approximation formula, Eq. (7), the degree of acceptable digital error must be carefully assessed. This is critical for intermediate-scale non error-corrected quantum processors, where the increase in the number of gates which comes with the increase in the number $n$ of Trotter steps cannot proceed indefinitely without affecting the quality of the results. In practical cases, it is usually sufficient for the digital error to be just smaller than the hardware noise. If $\mathrm{O}_{1}$ and $\mathrm{O}_{2}$ are two operators such that $[\mathrm{O}_{1},\mathrm{O}_{2}]\neq 0$, the so called first-order Suzuki-Trotter formula gives hatano_finding_2005

[TABLE]

A better scaling of the digital error can be obtained at the cost of an additional factor per iteration using the second-order formula hatano_finding_2005

[TABLE]

In both cases, a ratio $r_{\epsilon}=\delta^{p}/n^{q}$ controls the digital error as a function of the target evolution phase and the number of Trotter steps. Two different strategies can therefore be envisioned.

On the one hand, one could aim at a fixed digital precision $\epsilon$ over the whole range of the dynamical simulation. This requires to increase the number of Trotter steps, and consequently the total length of the quantum circuit to be computed, keeping the ratio $r_{\epsilon}$ fixed. As an example, for the first-order formula in Eq. (46) we get

[TABLE]

Notice that while the number of digital steps increases, the phase evolution $\delta_{n}=\delta/n$ required in each step decreases as $1/n$, thus keeping the overall computation time on the physical hardware linear in the total phase provided that each digital step can be implemented with a coherent operation of duration $t\propto 1/\delta_{n}$ lloyd_universal_1996 .

On the other hand, when the maximum length of quantum circuits that can be faithfully realized is de facto limited, such as in state-of-the-art noisy quantum processors, it might be convenient to follow a different approach, namely to keep fixed the length of the quantum circuit (i.e. the number of steps $n$). This produces a phase-dependent digital error scaling e.g. with $\delta^{2}$ in the first-order case. The fixed computational complexity, and consequently the uniform effect of hardware noise over the whole simulation, comes at the cost of a limited range of phases (and therefore of physical times) in which the results of the simulation agree with the target model. Hybrid solutions are also possible, e.g. by selecting reasonable number of steps $n$ in different intervals of phases $\delta$, always with the primary goal of balancing the total error arising both from the hardware noise and software-level approximations. In Fig. 2 we compare the two different approaches (fixed $\epsilon$ or fixed $n$) by showing how the fidelity $\left|\langle\psi_{0}|\psi_{n}(\delta)\rangle\right|$ of the digitally evolved state $|\psi_{n}(\delta)\rangle=\left(\exp\left(\mathrm{O}_{1}\frac{\delta}{n}\right)\exp\left(\mathrm{O}_{2}\frac{\delta}{n}\right)\right)^{n}|\psi_{0}\rangle$ with respect to the exact evolution $|\psi_{ex}\rangle=\exp\left(\left(\mathrm{O}_{1}+\mathrm{O}_{2}\right)\delta\right)|\psi_{0}\rangle$ decreases at long evolution times $t\propto\delta$ when $n$ is fixed or increases only linearly with the phase. In this simple 2-qubit case, we choose $\mathrm{O}_{1}=-i(\sigma_{x}^{(1)}+\sigma_{x}^{(2)})$ and $\mathrm{O}_{2}=-i\sigma_{z}^{(1)}\sigma_{z}^{(2)}$, corresponding to interaction terms typical of an Ising model in a transverse field.

III.6 Extracting physical observables

At the end of a quantum simulation, the final state $|\psi(t)\rangle$ of the quantum register is measured to retrieve information about the physical properties of the system under study. With an appropriate mapping of the generic observable of interest, $\mathbb{O}$, onto a combination spin-$1/2$ operators, the expectation value $\langle\mathbb{O}(t)\rangle=\langle\psi(t)|\mathbb{O}|\psi(t)\rangle$ can be reconstructed by a readout procedure combining, e.g., appropriate unitary operations $\mathrm{U}_{\text{meas}}$ and measurements in the computational basis. The reason why $\mathrm{U}_{\text{meas}}$ might be needed is that the eigenstates of $\mathbb{O}$ are, in general, different from the computational basis states: for example, if $\mathbb{O}=\sigma_{x}$ for a single qubit, the readout of $\langle\sigma_{x}(t)\rangle$ can be done by performing a Hadamard gate (i.e., mapping $\sigma_{x}\mapsto\sigma_{z}$) followed by a standard measurement in the computational basis. Joint qubit measurements are also possible, in general, as a way of characterizing the output quantum state salathe_digital_2015 .

More refined strategies allow the extraction of complex physical quantities and to optimize the efficiency of the measurement process. Here we will review in particular ancilla-assisted observation of dynamical correlation functions and of the spectrum of an Hermitian operator. This topic is discussed in detail in Ref. somma_simulating_2002, .

Given a $N$-qubit state $|\psi\rangle$, a Hamiltonian $H$ generating time evolution and two unitary operators $\mathrm{V}$ and $\mathrm{W}$, we define the dynamical correlation $\mathcal{C}_{\mathrm{VW}}(t)$ function as the quantity

[TABLE]

The quantum circuit in Fig. 3 describes how to compute $\mathcal{C}_{\mathrm{VW}}(t)$ using a quantum register and an ancilla qubit $a$. Here we assume that the quantum register is already prepared in the desired state $|\psi\rangle$, e.g. the ground state of the target physical system, and that the ancilla starts in the quantum superposition $\sqrt{2}|+\rangle=|0\rangle+|1\rangle$. The joint initial state of the quantum register $R$ and the ancilla is therefore $|\phi\rangle_{aR}=|+\rangle_{a}|\psi\rangle_{R}$. The first step is a $\mathrm{W}$ unitary performed on $R$ and controlled by the ancilla:

[TABLE]

A quantum circuit implementing the digital simulation of the time evolution $\mathrm{U}(t)=e^{iHt}$ is then applied to the quantum register to evolve the state $|\psi\rangle$, thus leading to

[TABLE]

Finally, a $\mathrm{V}$ unitary is applied to $R$, controlled by the state $|0\rangle$ of the ancilla (this can be obtained by adding $\mathrm{X}\equiv\sigma_{x}$ quantum gates on $a$ before and after the standard controlled operation). The output state is:

[TABLE]

A measure of the observable $\sigma_{x}$ on the ancilla gives

[TABLE]

In a similar way, $\operatorname{Im}\left[\mathcal{C}_{\mathrm{VW}}(t)\right]$ can be obtained by measuring $\langle\sigma^{(a)}_{y}\rangle$, in a second run of the algorithm. In total

[TABLE]

where $2\sigma_{+}=\sigma_{x}+i\sigma_{y}$. The same scheme can be applied to equal-time correlations by removing the unitary evolution or by moving it at the beginning of the circuit to evolve some initial state. It is worth noting explicitly that the useful information at the end of the proposed procedure is accessible through the ancilla $a$ alone, while the larger quantum register $R$ needs not to be measured at the end. The algorithm can also be generalized efficiently to the extraction of $n$-point time-correlation functions pedernales_efficient_2014 and of the expectation value of any operator which can be expressed as $\mathbb{O}=\sum_{j}c_{j}\mathrm{V}_{j}^{\dagger}\mathrm{W}_{j}$ where $\mathrm{V}_{j},\mathrm{W}_{j}$ are unitary operators somma_simulating_2002 .

With the addition of a classical Fast Fourier Transform (FFT), the strategy described above for time correlation functions can be used to extract the spectrum of a Hermitian operator $\mathrm{Q}$. The most relevant example in physical problems is certainly $\mathrm{Q}=H$, for some Hamiltonian of interest $H$. The hybrid quantum-classical approach, first proposed in Ref. somma_simulating_2002, and then further developed and applied (see e.g. Ref. chiesa_quantum_2019, ), requires the quantum register $R$ to be initialized in a state $|\psi\rangle$ with some overlap with the eigenstates $|Q_{l}\rangle$ of $\mathrm{Q}$

[TABLE]

Since by hypothesis the target operator is Hermitian, its exponential $\mathrm{U}_{\mathrm{Q}}(\theta)=e^{-i\mathrm{Q}\theta}$ is a unitary operator. This can be realized on the quantum register in exactly the same way as any standard time-evolution operator $\mathrm{U}_{H}(\theta)=e^{-iH\theta}$. We can then compute the expectation value $\langle\psi|\mathrm{U}_{\mathrm{Q}}(\theta)|\psi\rangle$ with the ancilla-based protocol described in the prevoius paragraph, setting e.g. $|\psi\rangle_{R}=|\psi\rangle$, $t=0$ (i.e. removing the time evolution $\mathrm{U}(t)$ part in Fig. 3), $\mathrm{W}=\mathrm{U}_{\mathrm{Q}}(\theta)$ and $\mathrm{V}=\mathbb{I}$. In general, the result will be of the form

[TABLE]

where $q_{l}$ are the eigenvalues of $\mathrm{Q}$. Applying FFT to the variable $\theta$ then yields

[TABLE]

III.7 Examples

The Hamiltonian for a 2-qubit isotropic Heisenberg model is

[TABLE]

The induced time evolution then reads

[TABLE]

where $\delta=Jt$ and the second equality, which is essentially the ST formula for $n=1$, follows from $[\sigma_{\alpha}^{(1)}\sigma_{\alpha}^{(2)},\sigma_{\beta}^{(1)}\sigma_{\beta}^{(2)}]=0\,\forall\alpha,\beta$. Recalling Eq. (35) and the results in Sec. III.3, a 6-$\mathrm{CNOT}$ decomposition for arbitrary $\delta$ can be given using the universal set $\mathcal{S}_{1}$, as shown in Fig. 4a. An equivalent and more efficient circuit in terms of number of two-qubit operations can be designed, according to the results discussed in Ref. vidal_universal_2004 , if we consider the time evolution operator globally as a single two-qubit transformation, see Fig 4b. Within $\mathcal{S}_{2}$, besides juxtaposing gate sequences of the form shown in Eq. (37), an optimal decomposition, using again only three 2-qubit gates instead of six, is reported in Fig. 4c, based on the identity las_heras_digital_2014 ; ferrando-soria_switchable_2016

[TABLE]

where $H_{\alpha\alpha\beta\beta}=\sigma_{\alpha}^{(1)}\sigma_{\alpha}^{(2)}+\sigma_{\beta}^{(1)}\sigma_{\beta}^{(2)}$. In $\mathcal{S}_{3}$ a decomposition with three $\mathrm{C}\Phi(\delta)$ follows immediately from Eq. (39) and single qubit changes of reference frame. Finally, in $\mathcal{S}_{4}$ a possible realization of the Heisenberg interaction can be obtained for some digital resolution $\delta$ as

[TABLE]

where $\mathrm{A}=\mathrm{T}_{4}(\delta,0)$, $\mathrm{B}=\mathrm{T}_{4}(\delta,\pi/2)$ and $\mathrm{C}=\mathrm{T}_{3}(\pi/4,\pi/2)$. With any of the above elementary decomposition in quantum gates, the digital quantum simulation of the 2-qubits Heisenberg model can be performed and physical information can be extracted by using the methods discussed in Sec. III.6. In a numerical example reported in Fig. 5a we show the digital quantum simulation of the individual magnetization of the two spins, which can be extracted by measuring the observable $\sigma_{z}^{(i)}$ and using the definition $\langle s_{z}^{(i)}\rangle=(1/2)\langle\sigma_{z}^{(i)}\rangle$. No digital error is present in this case.

The decomposition of the 2-spin Heisenberg model into elementary quantum gates presented above can be used as a building block, in combination with single qubit rotations, to perform more complex digital quantum simulations. A 3-spin Heisenberg chain with open ends and $N_{b}=2$ bonds, put in an external field, has a Hamiltonian of the form

[TABLE]

where

[TABLE]

describes a magnetic field oriented along the $z$-direction and each of the spin-spin bonds corresponds to a term

[TABLE]

In general, the two bonds can be nonequivalent, i.e. $J_{12}\neq J_{23}$. Since $[H_{\text{Heis},2}^{12},H_{\text{Heis},2}^{23}]\neq 0$ (independently from the coupling constants $J_{ij}$), the quantum simulation must be carried out using the ST digital procedure, alternating the application of the results presented for the 2-spin case on the two bonds

[TABLE]

where $\delta_{ij}=J_{ij}t$. The part describing the magnetic field on equivalent spins (we set the gyromagnetic ratio $g_{1}=g_{2}=g_{3}=g$) corresponds to single qubit rotations around the $z$ axis. Since this part commutes with the rest, it can be performed at the beginning of the circuit without any phase discretization. In Fig. 5b we show how these results can be used to compute the time evolution of the occupation probability of an initial state $|\psi_{0}\rangle=|100\rangle$.

As a third example, we recall that the Hamiltonian of the transverse field Ising model (TIM), introduced in Eq. (27), in the two-qubit case can be written as

[TABLE]

where

[TABLE]

The quantum simulation of the transverse field Ising model corresponds to the following digital process

[TABLE]

Apart from straightforward single qubit rotations around the $x$ axis, the required quantum circuit contains only $\mathrm{ZZ}$ operations, which can easily be translated into elementary quantum gates as shown in Sec. III.3. The time evolution of the total magnetization of the spin dimer along $z$ can then be extracted by measuring the expectation values of $\sigma_{z}^{(i)}$, see Fig. 5c.

Finally, we also report the example of a 3-spin open Heisenberg chain with an application of the ancilla-based algorithm discussed in Sec. III.6 to the extraction of spin-spin dynamical correlations $\mathcal{C}_{ij}^{\alpha\beta}(t)=\langle s_{\alpha}^{(i)}(t)s^{(j)}_{\beta}\rangle=(1/4)\langle\sigma_{\alpha}^{(i)}(t)\sigma^{(j)}_{\beta}\rangle$ on the system ground state. The latter, for the model under study and for a sufficiently strong external field, $B$, is well approximated by $|\psi\rangle=\left|\downarrow\downarrow\downarrow\right\rangle$, which is then assumed as the initial state on the quantum register. The structure of the required quantum circuit is shown in Fig. 6 for the case of next-to-nearest neighbors cross correlations. Autocorrelations and nearest neighbors correlations can be computed in a similar way by changing the target qubit involved in the operations controlled by the ancilla. Numerical results based on $\mathcal{S}_{1}$ decompositions are presented in Fig. 7.

Despite the relatively small size of the systems presented in the previous examples, all the elements introduced in this section can be used as basic modules to extend the quantum simulation to an arbitrary number of spins with pairwise interactions. When scaling up any spin chain to a larger number of interacting elements, i.e., with $N_{b}>2$ and possibly to different inter-qubits connectivity, one may also consider that all the terms involving disjoint sets of locally interacting spins generate independent time-evolution terms commuting with each other. These can then be simulated in parallel by using the elementary quantum circuit building blocks defined above, thus reducing the overall complexity of the quantum simulation (i.e., depth of the quantum circuit). As an example, let’s consider a chain of 4 spins labelled from $1$ to $4$ with pairwise $1-2$, $2-3$, and $3-4$ interactions: the $1-2$ and $3-4$ blocks can be run in parallel santini_molecular_2011 ; chiesa_digital_2015 . It is also worth mentioning that, concerning the simulation of dynamical correlation functions, in any $N$-spin system there are $\mathcal{O}(N^{2})$ two-body sigma correlations of the form $\langle s_{\alpha}^{(i)}(t)s^{(j)}_{\beta}\rangle$. These quantities, which are often of great physical interest chiesa_quantum_2019 , can then be extracted efficiently, in principle, with the ancilla-based methods discussed in Sec. III.6, e.g. by repeating a polynomial number of times the calculation with slightly modified circuits for each spin pair.

As mentioned in Sec. II.2, many problems in digital quantum simulations can be reduced to an equivalent spin Hamiltonian, which can in turn be directly mapped onto a $N$-qubit quantum register. A significant example is provided by fermionic systems, which play a central role in fields such as quantum chemistry, solid state physics and material sciences and, at the same time, are typically hard to treat with classical computational methods. To understand how the mapping $\mathcal{H}\mapsto H(\{\sigma_{\alpha}\})$ is carried out in such a case, we consider an elementary two-sites Fermi-Hubbard model:

[TABLE]

where $h.c.$ means the Hermitian conjugate and the operator $c^{\dagger}_{i,s}$ ($c_{i,s}$) creates (annihilates) a fermion with spin $s\in\{\uparrow,\downarrow\}$ on site $i\in\{1,2\}$. In $\mathcal{H}$, the coefficients $V$ and $U$ are hopping and on-site repulsion energies. Fermionic operators obey the canonical anticommutation rules

[TABLE]

The mapping is obtained by applying the Jordan-Wigner transformation jordan_uber_1928 ; bari_classical_1973 ; ortiz_quantum_2001 ; somma_simulating_2002 , which in this case takes the explicit form barends_digital_2015 ; lamata_digital-analog_2018

[TABLE]

Here, as usual, we use $2\sigma_{+}=\sigma_{x}+i\sigma_{y}$ and the anticommutation relations are preserved due to the properties of Pauli matrices. The resulting Pauli Hamiltonian reads

[TABLE]

Since $H$ is now the sum of local single- and two-body Pauli terms, the corresponding time evolution operator can be digitally simulated on a 4-qubit quantum register, using the techniques and the building blocks described above for paradigmatic spin-$1/2$ interactions.

While a one-dimensional fermionic chain with on-site and nearest-neighbor interactions can be mapped onto a spin Hamiltonian including only one- and two-spin terms santini_molecular_2011 , the application of the Jordan-Wigner transformation in more than one dimension or with long-range couplings leads to multi-spin interactions in the resulting Hamiltonian. Indeed, in the most general case, the mapping of Eq. (71) takes the form $c^{\dagger}_{j}=\left(\prod_{l=1}^{j-1}-\sigma_{z}^{(l)}\right)\sigma_{+}^{(j)}$ somma_simulating_2002 , where we have introduced a single label $j=\{i,s\}$ for fermionic modes, addressing both site and spin variables. The resulting spin Hamiltonian contains several terms, depending on the number of occupied states existing between each pair of interacting fermions after lattice sites have been properly ordered chiesa_digital_2015 . The quantum simulation of these terms is demanding, in terms of quantum circuit depth, for any architecture implementing only two-qubit interaction terms. Hence, in the NISQ era, the practical quantum simulation of fermionic models can strongly benefit from native many body terms in the hardware Hamiltonian, as it is the case for the trapped ion quantum processors to be briefly recalled in the following.

IV Experimental achievements and prospective technologies

The last few years have represented a timeline of intense development towards the realization of quantum computing architectures. Among the plethora of possible platforms, two leading technologies have been emerging: trapped ions and superconducting quantum circuits. Here we give a brief overview of the main achievements reported to date in these two experimental set-ups. The focus will be on the perspective development of these two platforms as universal quantum simulators, rather than reviewing the technical details of the different setups, for which thorough reviews already exist. We notice that a comparison between the two leading platforms can be done based, for instance, on relevant figures of merit such as gate fidelities or the ratio between coherence and gating times, roughly representing the number of operations that the quantum hardware can reliably perform before the quantum information is degraded. However, it is important to stress that most of the remarkably high experimental fidelities reported in the literature have been achieved on two-qubit setups that, although in principle scalable, are optimized for a specific target and are often pushed to the experimental limits. Conversely, when several qubits are connected and operated together, many new challenges emerge, such as the need to selectively address only some of them or to stem cross-talk errors. As discussed below, this yields slower and more error-prone gates. Hence, a fair study on the actual performances of a quantum simulator should be done on an architecture including several interacting qubits, in order to assess the effective scalability of the given setup. A useful concept, in this respect, is that of “quantum volume”, recently introduced by IBM researchers with the aim of quantifying with a single paramenter the performance of a quantum computer, based on how efficiently a quantum algorithm can be run moll_quantum_2018 . The quantum volume takes into account both width (i.e., the number of qubits used) and depth (i.e., the number of reliable operations performed) of a given quantum circuit. Despite a few non trivial issues related, e.g., to individual qubit addressing and frequency crowding, adding qubits to the actual hardware mostly appears to be a technological challenge, while implementing a sequence of high-fidelity operations on several qubits could be more demanding. Indeed, this requires a considerable improvement of gate fidelities, suppression of both qubit decoherence and coherent errors due to imperfect qubit manipulations, and reduction of unwanted qubit-qubit interactions (cross-talk) whose harmful effect increases with the system size. Increasing the quantum volume truly represents the next challenge for current NISQ devices to evolve into useful quantum simulators.

In Fig. 8, we try to give a quantitative summary of the main experimental achievements reported in the recent literature. As mentioned above, it is generally difficult to directly compare experiments performed on different platforms, under different initial conditions, and reporting slightly different figures of merit. However, the graph aims at giving a visual idea of the scenario on the digital quantum simulation of spin Hamiltonians up to date, and the capabilities offered in different quantum hardware. Evidently, there is still a considerable correlation between the number of digital steps included in the simulation and the fidelity of the final state obtained. As a general comment, trapped ions quantum simulators allow deeper quantum circuits with better performance, i.e. a larger number of Trotter steps is possible. This is in line with recent studies comparing the two platforms when challenged with similar quantum algorithms on 5 qubits processors linke_experimental_2017 . One may notice that 5 Trotter steps are currently a limiting value for superconducting circuit quantum simulators, where the fidelity drops to values slightly above $60\%$, still far from an acceptable result, also in view of scalability. Finally, in terms of size of the simulated model we see that spin-models with up to 6 spins have been simulated on trapped ions processors lanyon_universal_2011 , while up to 4 spins on superconducting ones chiesa_quantum_2019 .

As a last comment, it is worth mentioning that great efforts have been lately devoted to design error-mitigation strategies to improve the quality of the final results from the data directly extracted from the quantum hardware. Some of them are problem specific, based, e.g., on the symmetries of the target Hamiltonian mcardle_error_2019 , or on general properties of the measured observables chiesa_quantum_2019 , which set constraints on the output of the quantum calculation. In particular, the stabilizer-like method proposed in Ref. mcardle_error_2019 could enable the detection of up to 60-80 % of depolarizing errors, thus largely improving the results of variational algorithms, as well as of quantum simulations. Theoretical proposals have also been developed, demonstrating that the accuracy in the expectation values of computed quantum observables can be improved by suitably interpolating to zero noise the results on a series of experiments at varying noise levels PRX_Benjamin_2017 ; PRX_Benjamin_2018 . This strategy has been recently applied to suppress incoherent errors in variational calculations on superconducting circuits kandala_error_2019 (see below). The combination of these two approaches could further improve the accuracy of the computed observables mcardle_error_2019 .

We will now give a more detailed description of each of the two aforementioned experimental platforms, and the corresponding key results on the digital quantum simulation of spin models. A word of warning should be given, however, before any further discussion takes place: a number of challenges still need to be faced before either one or the other may ultimately develop into a fully operational quantum simulation technology with clear advantages over classical simulations of quantum many body models. In fact, while scalability requires that the number of qubits on a given quantum hardware be arbitrarily increased without degrading the reliability of individual preparation and read-out, as well as quantum gating fidelities, there is no current quantum technology that is truly fulfilling this requirement. The key challenges to be faced in the forthcoming years are thus related to practically addressing these issues, even before fault tolerant quantum computation be developed.

IV.1 UQS with trapper ions

Atomic ions trapped in a suitable combination of static and time-varying electric fields have been representing one of the most promising routes towards realizing fully operational quantum processors since the late nineties blatt_entangled_2008 ; blatt_quantum_2012 ; monroe_scaling_2013 . As compared to neutral atoms, the trapping potential for charged atomic species can be much stronger, thus allowing to hold each single ion for several hours, and even days, with very long coherence times and a remarkably high degree of external control. This makes them a reliable solution for quantum information processing satisfying all of the required DiVincenzo criteria (for a comprehensive recent review on the state-of-art of this technology focused on quantum computing, we refer to Ref. bruzewicz_review-trap-ion_2019, ). Today, digital quantum processors with up to 11 programmable ions in a linear Paul trap have been made available wright_benchmarking_2019 , also through spin-off companies such as IonQ offering restricted cloud access (see https://ionq.com/, accessed 2019-09-01). Other companies are working on providing access to general purpose trapped ion quantum computing hardware with up to 20 qubits, such as Alpine Quantum Technologies (see, e.g., https://www.aqt.eu/, accessed 2019-09-01). In such realizations the trapping is generally created through radio frequency oscillating electric fields that generate a stable two-dimensional potential well for the charged atoms leibfried_review-trap-ion_2003 , such that ions are then confined by a further harmonic trap along a linear chain, in which the single particles are stably kept a few microns apart due to the mutual Coulomb repulsion monroe_scaling_2013 . While this is the configuration that currently allows the most performing type of apparatus as programmable NISQ devices, it is considered hardly scalable towards larger numbers of trapped ions within the same chain, as well as to higher dimensionality of the confining potential, for which new solutions will have to be realized. In particular, modularity bautista_multilayer-iontraps_2019 and 2D arrays bruzewicz_2Darray_ions_2016 are being considered as valuable approaches, for which we will likely see progress in the near future. It is worth noticing that this type of hardware requires ultra-high vacuum and possibly laser cooling, but not necessarily cryogenic apparatus to be operated.

Trapped ions may be employed as reliable, coherently manipulated qubits, in which the information is encoded into the internal quantum states of the charged atomic species. While we refer the interested reader to more specialized review works on trapped ion quantum technologies leibfried_review-trap-ion_2003 ; blatt_entangled_2008 ; schindler_quantum_2013 , we hereby summarize its main aspects for completeness, especially concerning the ongoing developments towards programmable NISQ devices. First, we remind that different types of trapped ion qubits can be defined, depending on the frequency spacing between the relevant energy eigenstates selected to encode the logical $|0\rangle$ and $|1\rangle$. In general, hyperfine qubits are encoded into a pair of hyperfine energy levels typically separated by GHz frequencies, while optical qubits are defined corresponding to quadrupole active transitions in the hundreds THz range; also Zeeman qubits can be defined by application of an external static magnetic field, opening a low frequency (few MHz) and tunable gap between magnetic levels. Pros and cons characterize each specific choice, and a dedicated experimental apparatus must necessarily be developed for one qubit type or the other. Among the different possibilities, particularly advanced appear the technologies based on 40Ca+ as optical qubits friis_observation_2018 ; hempel_quantum-chemistry_2018 , and 43Ca+ or 171Yb+ as hyperfine qubits ballance_fidelity-gates_2016 ; debnath_demonstration_2016 , respectively, although several other atomic species with a single outer electron can be successfully trapped monroe_scaling_2013 ; bruzewicz_review-trap-ion_2019 . In fact, irrespective of the specific qubit realization, initialization and readout of their quantum state is performed through coherent manipulation, e.g. by using external lasers of suitable frequency to optically pump the ion state into the desired one (preparation), or by detecting resonantly scattered radiation from an optical transition (readout). The qubits’ individual control (i.e., single-qubit operations) can be performed by directly coupling the $|0\rangle$ and $|1\rangle$ eigenenergy levels, for which the required tools inevitably depend on the qubits type (e.g., hyperfine qubits require microwave control or stimulated Raman coupling, while optical qubits internal states are directly coupled through a resonant laser) leibfried_review-trap-ion_2003 . Two-qubit entangling gates between ions trapped along the same chain are realized by exploiting the transverse normal vibrational modes of the whole ion string trapped in a harmonic potential cirac_quantum_1995 , which are used as a bus to transfer quantum information. Currently, the most performing multi-qubit operations rely on a controlled-phase type of gate that was originally proposed from Mølmer and Sørensen molmer_multiparticle_1999 , whose formal expression is given in Eq. (43). The details of the specific implementation of these gates depend on the type of qubit schindler_quantum_2013 ; shapira_robust_2018 ; webb_resilient_2018 . In general, they require only global control lasers to entangle multiple qubits, hence avoiding beams to be independently focused on each ion (which may anyway be necessary, e.g., for single-qubit control). In hyperfine qubits-based quantum hardware, a Mølmer-Sørensen gate with a global laser and suitably detuned individually addressing beams allow to realize an effective two-qubit XX gate, which then becomes the native two-qubit gate on that hardware debnath_demonstration_2016 . Elementary combinations of single-qubit rotations and such XX two-qubit gate finally lead to the CNOT gate nielsen_quantum_2000 . Independently of the specific implementation, an unequivocal advantage of any trapped ion quantum hardware is that each qubit can be connected (and entangled) to any other qubit in the chain, thus practically realizing all-to-all connectivity.

In terms of absolute performance specifically meant for prospective applications in quantum computing, the different experimental platforms have essentially shown quite comparable figures of merit. In particular, when working with isolated or pairs of trapped ions, single-qubit operations with fidelities in the order of $99.9999\%$ have been shown harty_high-fidelity_2014 , as well as two-qubit gates reaching fidelities in excess of $99.9\%$ even in different experimental setups gaebler_fidelity-gate-set_2016 ; ballance_fidelity-gates_2016 . Typical duration for single-qubit gates varies between 100 ns and few tens of $\mu$s, while two qubit gating times range in the $\mu$s to few hundred $\mu$s interval ballance_fidelity-gates_2016 ; schafer_fast_2018 ; bruzewicz_review-trap-ion_2019 ; readout is typically performed in hundreds of $\mu$s with fidelities in the $99.99\%$ range meyerson_readout_2008 . Generally speaking, a trade-off between speed and fidelity has to be found, where faster usually means more error prone. Considering also that the reported coherence times of trapped ion qubits vary between few hundred ms bermudez_progress_2017 and hundreds of seconds (i.e., up to several minutes) harty_high-fidelity_2014 ; wang_memory_2017 , depending on the type of qubit, the trapped ion quantum hardware is the one currently allowing to achieve the highest coherence vs. gating time ratio, ideally in the order of $10^{5}$ to $10^{6}$ (e.g., assuming two-qubit gate times in the order of 100-200 $\mu$s, and depending on the qubit type).

It is worth emphasizing, however, that the considerations above are mostly limited to few qubits quantum hardware: while analog quantum simulators of Ising chains have been shown with more than 50 trapped ions zhang_observation_2017 , digital quantum simulators are quite challenging to scale up. In fact, increasing the number of ions in the chain ultimately limits gates fidelity, mainly due to the difficulties of individually addressing single qubits while avoiding cross-talks between the different beams. Moreover, errors affecting the overall success of a quantum computation generally arise from two distinct mechanisms: decoherence, e.g. from undesired qubit-environment coupling such as spontaneous emission, frequency shift, or motional heating, and imperfect control fields, such as miscalibrated or noisy control field amplitude, frequency, or polarization, which typically result in quantum gate errors. Targeted strategies to mitigate these issues will become crucial to further improve the trapped ion quantum hardware bruzewicz_review-trap-ion_2019 . On a quantitative level, it has recently been reported that the fidelity of multi-qubit entangling Mølmer-Sørensen gates degrades from $99.6\%$ with 2 optical qubits down to $86\%$ when the same quantum hardware is loaded with 10 optical qubits erhard_characterizing_2019 . On the other hand, a quantum processor based on 11 hyperfine qubits has recently been shown to achieve an all-to-all connectivity with single- and two-qubit XX gate fidelities of $99.5\%$ and $97.5\%$ on average, respectively wright_benchmarking_2019 . It is interesting to notice that these results improve previous reports on an analogous 5-qubits quantum processor linke_experimental_2017 , as a clear signature of the ongoing development.

After reviewing the basic aspects of trapped ion quantum hardware, and the up-to-date figures of merit currently achieved in state-of-art experimental platforms, we hereby show a few examples of digital quantum simulation performed on one such a NISQ device. An example of results showing the first universal quantum simulations performed on NISQ processor based on optical qubits is given in Fig. 9, as reported back in 2011 in Ref. lanyon_universal_2011, from Lanyon et al. This represents a seminal work since it shows that the same quantum hardware can be experimentally reprogrammed to run the digital quantum simulation of different spin models and interaction terms, even for Hamiltonian terms not natively implemented on the simulator. Each spin-1/2 is directly mapped onto a single ionic qubit, and unitary operations ($C$, $D$, $E$, and $F$) are defined in terms of the universal set of gates in Eqs. 41-43, which is native on this hardware (see original reference for details). The time evolution is quantified by a dimensionless phase $\theta=Et/\hbar$, which is reported on the horizontal axis. The initial state is chosen as an eigenstate of $\sum_{i}\sigma_{x}^{(i)}$, and the population in each of the eigenstates is monitored as a function of $\theta$. The same work reported digital quantum simulation of up to 6 spins and multi-spin interaction terms, allowing to envision the potentialities of digital quantum simulators for fundamental physics studies. Hence, the work from Lanyon et al. sets a reference standard for all the following demonstrations of digital quantum simulation on NISQ processors. At time of writing, we are not aware of any digital quantum simulation of spin models performed on quantum processors based on hyperfine qubits such as the one presented, e.g., in Ref. wright_benchmarking_2019, , although it would be interesting to test them with some of the algorithms presented in this review, indeed.

The quantum hardware employed to obtain the results of Fig. 9 has more recently been applied to the quantum simulation of the real time dynamics underlying particle-antiparticle pair creation in lattice gauge field theories martinez_real-time_2016 , where up to 4 qubits were employed to run the related quantum circuit. There, digital quantum simulation is obtained after mapping the fermionic degrees of freedom into Pauli spin operators applying the Jordan-Wigner transformation, as outlined before. Given the size of the quantum register, a toy model successfully simulating the electron-positron spontaneous creation from vacuum fluctuations and the persistence of their entanglement was reported, which creates a bridge between digital quantum simulators and elementary particle physics.

Without restricting to the main object of the present review, it is worth mentioning recent results obtained on trapped ion quantum technology. First, a 20 qubits register has been shown to reliably allow for the creation of multi-qubit entangled states friis_observation_2018 , thus opening the door to quantum simulations of larger spin systems. The same quantum hardware has been used to show a hybrid quantum-classical approach to the simulation of the Schwinger model kokail_self-verifying_2019 . The latter belongs to the class of variational optimization algorithms commonly employed in quantum chemistry mcclean_theory_2016 , an approach now commonly defined Variational Quantum Eigensolver (VQE). In fact, the latter methods have been recently applied to accurately calculate the ground state energy of simple molecules hempel_quantum-chemistry_2018 . Conceptually similar approaches have also been applied to the quantum simulation of effective field theories in nuclear physics, such as calculating the deuteron nucleus binding energy with percent accuracy and using record depth quantum circuit on a ion trap quantum processor shehab_toward_2019 . A VQE algorithm has been recently described and applied to obtain the ground state energy of a large chemical complex, such as the water molecule, by using a hyperfine qubits-based quantum processor nam_water_2019 .

In summary, despite considerable progress and the maturity reached, many challenges still need to be addressed before practical quantum computers based on trapped ion technology may allow universal quantum simulations outperforming classical computations, in perspective. It currently appears a significant challenge to operate a quantum computer with a number of fully controlled qubits larger than a few tens, still preserving the required figures of merit in terms of gates fidelity and coherence properties, which might still need a few years of intense research and development.

IV.2 UQS with superconducting circuits

Superconducting quantum circuits have lately emerged as a practical quantum computing technology after a fast development in the last decade blais_cavity_2004 ; wallraff_strong_2004 ; koch_charge-insensitive_2007 ; gambetta_protocols_2007 ; majer_coupling_2007 ; gambetta_quantum_2008 ; mariantoni_implementing_2011 ; rigetti_superconducting_2012 ; barends_superconducting_2014 ; nersisyan_manufacturing_2019 . In fact, this platform has reached the level of reliability typical of trapped ions quantum processors in only few years of continuous improvement linke_experimental_2017 . It is worth mentioning that some of the worldwide leading high-tech companies that are currently targeting the realization of NISQ hardware are concentrating their efforts on this technology. As for the trapped ion case, restricted cloud access is also provided to superconducting quantum processors of up to 20 qubits from IBM (see https://www.research.ibm.com/ibm-q/, accessed 2019-09-01) and up to 16 qubits from Rigetti Computing (see https://www.rigetti.com/, accessed 2019-09-01), respectively. These devices work with cryogenic set-up in a 3He/4He dilution refrigerator with 10-15 mK base temperature, in which qubits can be efficiently encoded into the anharmonic energy spectrum of the lowest collective charge/current excitations in a micro-LC resonator, with a nanostructured Josephson junction playing the role of a nonlinear inducting element. With respect to quantum hardware based on trapped atomic species, these architectures certainly bring along the advantages of a solid state chip-based microelectronic technology, although working at low temperature and with intrinsic limitations in inter-qubit connectivity. Several families of superconducting qubits have been realized (i.e., phase qubits, rf-sQUID, flux qubits, charge qubits, transmons), characterized by different ratios between the characteristic parameters of their quantized Hamiltonian, namely charging energy of a Cooper pair, inductive and Josephson energy. These systems are multi-level quantum oscillators whose energy spectrum is made sufficiently anharmonic to enable selective addressing of only the two lowest energy levels, which then become the $|0\rangle$ and $|1\rangle$ qubit states. While additional levels may become a resource to implement two-qubit gates Nature_DiCarlo_2009 ; mariantoni_implementing_2011 or quantum-error correction codes hussain_coherent_2018 , they also represent a possible source of leakage errors that must necessarily be limited by using, e.g., long control pulses for qubit transitions. We refer to Ref. wendin_quantum_2017, for a detailed description and a comparison between the different types of superconducting qubits. Here we only note that the evolution of these devices led to the realization of the transmon koch_charge-insensitive_2007 , which has allowed to reach coherence times in the 100 $\mu$s range for single qubits rigetti_superconducting_2012 . This is a development of the Cooper pair box into a circuit less sensitive to charge noise but still characterized by a sufficiently anharmonic energy spectrum. The frequency of the transmon can be tuned by varying the Josephson energy using a SQUID. The 2D transmon is nowadays the elementary unit of several scalable architectures corcoles_demonstration_2015 , but other transmon-like qubits have been realized, such as the Xmon barends_superconducting_2014 , which has already shown remarkable fidelities even in a setup consisting of 9 qubits Nature_Kelly_2015 ; Nature_Barends_2016 .

Different qubits in these platforms are typically interconnected through superconducting transmission line resonators, and they can be individually addressed through other transmission lines wired at the edges of the chip board. The latter allow to perform single-qubit initialization, manipulation, and read out through microwave pulses. Current superconducting quantum circuits allow for single-qubit gate fidelities above $99.9\%$ barends_superconducting_2014 . Each qubit is dispersively coupled to a resonator to mediate an effective qubit-qubit interaction, employed to implement two-qubit gates. These are obtained either (i) by tuning the qubit transition frequency by a local magnetic field, or (ii) by using a cross-resonant (CR) drive, in which microwaves resonant with a target qubit are applied on another control qubit. While the first realizations of scheme (i) were obtained by bringing the two qubits into mutual resonance to get a virtual photon exchange (thus yielding an effective XY interaction salathe_digital_2015 ), the most promising implementations are currently based on tuning one qubit along a “fast adiabatic trajectory”barends_superconducting_2014 that moves the $|11\rangle$ component of the wave-function close to the avoided level crossing with state $|02\rangle$, leading to a state-dependent phase and hence to the implementation of the controlled-phase gate Nature_DiCarlo_2009 . This approach results in very fast (40 ns) and high-fidelity two-qubit gates, compared to relaxation and coherence times in the 20-40 $\mu$s range on average barends_superconducting_2014 . Conversely, scheme (ii) uses fixed-frequency qubits in order to avoid frequency crowding and implements a CNOT gate by activating the resonator-mediated interaction via a cross-resonant drive chow_simple_2011 . This requires more selective pulses and hence results in slower two-qubit gates, taking $\sim 200-300$ ns chow_simple_2011 ; gambetta_building_2017 with only slightly longer coherence times (50 $\mu$s on average quality devices) and average fidelity of $96\%$. These results make platforms based on scheme (i) more promising in view of concatenating several gates to perform quantum simulation algorithms. Furthermore, the possibility to tune the phase of $C\Phi$ gates makes architecture (i) more flexible than (ii), and it reduces the circuit depth for the simulation of many model Hamiltonians of physical interest. In this respect, a novel proposal to directly implement exchange-type gates with tunable amplitude and phase on fixed frequency qubits (thus making also this architecture much more flexible for quantum simulations) has been recently reported PRApplied_Ganzhorn_2019 .

Readout still represents one of the main limitations of superconducting platforms. Indeed, accurate readout based on transmission measurements of the frequency shift of auxiliary resonators is usually slow, taking few hundreds of ns, with average accuracy in excess of $96\%$ gambetta_building_2017 . However, high-fidelity (i.e., more than $99\%$) readout was demonstrated in 140 ns on a four-qubit setup PRL_Sank_2014 .

As already noted for trapped-ion architectures, we finally point out that while remarkable results have been obtained in samples specifically aimed at testing the basic operations on few qubits circuits, thus reaching fidelities no more limited by unitary errors PRA_Gambetta_2016 , the realization of a scalable platform able to implement fast and high fidelity single- and two-qubit gates, as well as efficient readout, still represents a key challenge. Indeed, recent studies have demonstrated that systematic coherent errors due to an imperfect implementation of the elementary gates are still one of the leading error sources when several qubits are operated together in a quantum simulation chiesa_quantum_2019 . This will also become evident from the examples of digital quantum simulations reported below.

The first digital quantum simulations of spin models on superconducting quantum hardware were experimentally reported in 2015 salathe_digital_2015 . Here, the evolution of the spin magnetization in Heisenberg and Ising models was systematically studied for 2-spin type Hamiltonians on a 4-qubits quantum processor, as a function of the number of ST steps. Superconducting processors with tunable frequency qubits (through external flux bias lines) naturally implement a XY-type interacting spin Hamiltonian, which can be used as the basis to digitally program a full Heisenberg or Ising type evolution through a circuit model, as outlined in the previous section and explained in detail in the original references las_heras_digital_2014 ; salathe_digital_2015 . An example is reported in Fig. 10, where the digital evolution is explicitly shown for the two spins projections along the magnetic field direction, $z$, with up to 3 digital time steps, for an initial state prepared in $\left|\uparrow\right\rangle\left(\left|\uparrow\right\rangle-i\left|\downarrow\right\rangle\right)/\sqrt{2}$ that evolves non trivially in time. A summary of the fidelity obtained from these quantum simulations on the same quantum hardware with up to 5 Trotter steps is also reported from the original reference salathe_digital_2015 . While the ideal fidelity of the simulated quantum state with respect to the exact evolution increases against the number of Trotter steps, the experimental one starts to decrease after about 2 or 3 digitized steps, depending on the phase. It is evident that 5 digitized steps in the simulated time evolution still presented limited fidelities, due to the short coherence times and to systematic circuit errors. Nevertheless, such results have set a milestone as a proof of concept demonstration of universal quantum simulations in superconducting quantum circuits.

Correlation functions represent some of the most useful and informative quantities to be calculated in quantum many body physics. A first attempt at simulating the digital time evolution of two-point correlations was already reported in salathe_digital_2015 . More recently, dynamical correlation functions have been experimentally simulated on the quantum processors made freely available by IBM, through their IBM Quantum Experience (see quantumexperience.ng.bluemix.net/qx/experience) chiesa_quantum_2019 ; smith_simulating-manybody_2019 . These processors are based on fixed frequency superconducting qubits allowing the experimenter to implement arbitary single-qubit rotations and CNOT two-qubit gates, along the lines of scheme (ii) described above (see qiskit.org). These elementary operations are then combined as in Fig. 4 to obtain complex expectation values for the correlation functions. The results reported in Ref. chiesa_quantum_2019, are encouraging in view of scalability: dynamical correlations were digitally simulated for various basic spin models, ranging from Ising to isotropic and anisotropic exchange Hamiltonians, both for spin dimers and trimers chiesa_quantum_2019 . The largest number of ST steps that could be reliably simulated on the IBM quantum processor was $n=4$ for spin dimers, and $n=2$ for trimers (due to the larger depth of the corresponding quantum circuit). An example of the digital quantum simulation of time-dependent two-body correlation functions, as defined in the previous section, is reported in Fig. 11 for a three spin-1/2 Heisenberg model in external magnetic field along $z$ for an initial state $\left|\downarrow\downarrow\downarrow\right\rangle$, as compared to the ideal digitized evolution with $n=2$ ST steps, showing truly remarkable agreement. A sketch of the available 5-qubit quantum hardware is also shown, with an outline of chip connectivity. In Ref. chiesa_quantum_2019, , the largest quantum simulation reported on the actual IBM quantum hardware actually employed 5 qubits (4 encoding the target system, plus 1 ancilla for correlations readout), showing good agreement with the expected behavior despite the noisy nature of the quantum processors. It is important to note that the accuracy of these results was largely improved by systematic post-processing corrections based on general properties of the extracted correlation functions, a problem-specific procedure which could be applied to many other quantum simulations, after analysis of the symmetries of the target Hamiltonian mcardle_error_2019 . Fitting of such digitally simulated correlations allows one to extract four-dimensional inelastic neutron scattering spectra, a crucial experimental tool to characterize magnetic molecules baker_spin_2012 ; garlatti_portraying_2017 . The speedup of a quantum processor in simulating the dynamical correlations needed to compute the inelastic neutron cross-section could allow for an efficient and real time interpretation of experiments on complex molecules, a task which is nowadays infeasible with classical computer simulations. This study also reports a systematic analysis of the errors propagating on the quantum hardware, highlighting that one of the main limitations to circuit depth is currently represented by systematic off-resonant driving errors mckay_efficient_2017 . This could be a direct consequence of the many transitions that must be addressed on a chip consisting of several qubits with only slightly different transition frequency (the so-called frequency crowding). In fact, recent quantum simulations of similar models performed on a 20-qubit quantum hardware show that a proper accounting of such errors is essential to obtain a reasonable agreement with expected results smith_simulating-manybody_2019 .

Going beyond spin models, the superconducting quantum hardware has been tested as a UQS of the Fermi-Hubbard model las_heras_fermionic_2015 , which has been experimentally achieved with up to 4 modes barends_digital_2015 .

A hybrid approach combining adiabatic to digital quantum computation was reported in Ref. Nature_Barends_2016 , where a superconducting circuit made of frequency-tunable qubits and implementing two-qubit controlled-phase gates according to scheme (i) was used to probe the dynamics of a chain of 9 interacting spins, starting in a factorized state and evolving along an adiabatic trajectory by switching on an anisotropic exchange interaction term. Recent results in hybrid quantum-classical approaches, such as the VQE method already recalled at the end of the previous Section, have been efficiently employed to show quantum chemistry simulations on superconducting NISQ processors omalley_scalable_2016 ; moll_optimizing_2016 ; kandala_hardware-efficient_2017 , in which the ground state energy of multi-atomic molecules was calculated with precision approaching the chemical accuracy limits PRApplied_Ganzhorn_2019 . Using an evolution of the VQE algorithm, nuclear physics quantum simulations have also been reported in superconducting quantum hardware, with the cloud computing of the deuteron binding energy dumitrescu_cloud_2018 . Along the same lines of trapped ions quantum simulators martinez_real-time_2016 ; kokail_self-verifying_2019 , a similar quantum classical algorithm has been used to solve for the Schwinger model dynamics on a superconducting quantum hardware klco_quantum-classical_2018 .

Interesting comparisons between trapped ions and superconducting quantum processors in applying such hybrid quantum-classical approaches have been reported, in which the same algorithm was simulated on different platforms, showing a substantial equivalence of the two leading architectures when the same number of qubits could be used, but trapped ions processors allowing for a larger system size to be simulated shehab_toward_2019 . Compared to trapped-ion technologies, superconducting circuits show much larger gates speed, even if the ratio of coherence time to gate operation remains smaller. This can be a remarkable advantage when comparing the performance of quantum and classical devices in terms of absolute execution time of a given algorithm. While progress in the last few years has led to a substantial increase of the coherence times in transmon qubits, the leading errors in multi-qubit architectures are still coherent chiesa_quantum_2019 . Indeed, the need to address several qubits with only slightly different transition frequencies, as well as the small anharmonicities required to keep long coherence times, require long and often not completely resonant driving pulses PRA_Gambetta_2016 . In addition, inter-qubit interactions are never completely switched off, even when qubits are in the idle phase. This leads to unwanted evolution of the multi-qubit wave-function (cross-talk), an effect increasing with the system size and strongly depending also on the chip connectivity. For quantum simulation purposes, the latter should be as close as possible to that of the target system, in order to reduce the circuit depth and avoid cumbersome SWAP gates. In fact, trapped ions are more promising for establishing entanglement between distant pairs of qubits, and hence also to simulate fermionic Hamiltonians involving multi-qubit interaction terms, although the first attempts to go beyond nearest-neighbors coupling have been reported on three-qubit quantum processors roy_programmable_2018 . Finally, compared to trapped ion-based technologies in which all qubits are identical, superconducting qubits are all characterized by different parameters (such as transition frequencies and mutual couplings), which are also affected by thermal cycling and hence require a detailed and frequent characterization to accurately calibrate the control pulses klimov_thermal-cycling_2018 .

As mentioned above, considerable work is currently focusing on understanding the main sources of error, and on developing error mitigation techniques to enhance the overall quantum simulation fidelities temme_error_2017 ; mcardle_error_2019 ; chiesa_quantum_2019 . These theoretical proposals have been recently and successfully applied to improve the accuracy of the observables extracted from VQE calculations on a superconducting chip kandala_error_2019 . Although demonstrated on elementary single- and two-qubit gates in a VQE experiment, this protocol can be applied to mitigate errors of any quantum algorithm since it is problem-independent and does not lead to any hardware overhead. However, its application requires the experimenter to control the amount of noise on the hardware (which in Ref. kandala_error_2019 is mapped on the evolution under a scaled drive) and it is limited to incoherent errors.

In the ongoing effort to develop a fault tolerant quantum computing architecture, these results pave the way towards reaching the quantum advantage, possibly already within the NISQ time frame. As a summary of all the previous results, IBM tested the concept of quantum volume on three of its devices, finding that it doubled each year (from 2017 to 2019), from 4 to 8 to 16, a trend similar to Moore’s Law for classical computers and promising for future perspectives. Nevertheless, similar conclusions apply here as already given for the trapped ion quantum hardware: despite considerable progress and the first attempts at performing full digital quantum simulations, several challenges need to be addressed before quantitative and not only qualitative accuracy be reached, especially on larger system sizes and running deeper quantum circuits. Specifically, a deeper understanding of error and noise sources, as well as suitable strategies to mitigate them on hardware with larger number of qubits is crucial at this NISQ stage smith_simulating-manybody_2019 , which will be the focus of intense research and development in the coming years.

IV.3 Prospective technologies for UQS

While the two leading technologies outlined above are currently the mainstream in practical quantum computing, it is still unclear if they will be able to overcome the challenges to reach a large number ($N>100$) of logical qubits and a significantly larger amount of error correcting ones. Alternative technologies might also start playing a role in the meantime. In this respect, it is interesting to follow recent progress in semiconductor-based technologies. After the widespread success of semiconductors in microelectronics applications, they have been a little behind the scene in the quest for practical quantum computing devices despite the early proposals loss_quantum_1998 . In particular, semiconductor based quantum dots have long been considered as potential spin qubits. Single-spin read-out and manipulation has been shown quite early elzerman_single-shot_2004 ; petta_coherent_2005 , but scalability has been hindered so far, mostly due to coherence times being limited by nuclear spin dephasing and spin-orbit coupling hanson_spins_2007 . However, recent advances in silicon-based quantum dots have renewed interest in the actual possibilities of these technologies: CNOT gates between two quantum dots have been shown with about $78\%$ fidelity zajac_resonantly_2018 , and two-qubit gates with fidelities in the order of $80\%$ have also been reported, with single qubit rotation precision of $\sim 99\%$ veldhorst_two-qubit_2015 ; watson_programmable_2018 ; noiri_fast_2018 ; huang_fidelity_2019 . Gating times are below 100 ns, and dephasing times about 200 ns watson_programmable_2018 . These results, although still far from the required performances achieved, e.g., in the trapped ion or superconducting circuits architectures, are quite promising and set a stepping stone on the development of a fully semiconductor-based quantum technology, evidently interesting for a number of reasons (from cheap costs to the potential for mass scale manufacturing).

In parallel to research on quantum dots, controlled impurities and defect ions in silicon have been lately considered for a potentially low-cost quantum technologies awschalom_quantum_2013 . After the demonstration of single-spin read-out and manipulation of localized donor impurities in silicon morello_single-shot_2010 ; pla_single-atom_2012 , two-qubit quantum gates have been proposed for this prospective platform kalra_robust_2014 ; tosi_silicon_2017 . Very recently, the first $\sqrt{\mathrm{SWAP}}$ operation between Phosphorus donor electrons in silicon, with gating time of $800\,$ps and fidelity around $94$% has been reported he_two-qubits-Phosphorus_2019 . These results hold promise that further developments might be seen in the near future, if challenges related to scalability will be overcome.

While electronic states in engineered potentials or in impurity states are naturally emerging as reasonable candidates for a qubit-based architecture, it is less obvious that photons, in particular photonic integrated circuits, could play a significant role as UQS. Photonic circuits have been largely explored as analog quantum simulators aspuru-guzik_photonic_2012 . The main limitation to exploit photonic states for quantum computing lies in their weak interactions, due to intrinsically small material nonlinearities. While it may be argued that suitable electromagnetic confinement in nonlinear materials might lead to single-photon sensitivity ferretti_single-photon_2012 ; flayac_all-silicon_2015 , no such effect has been measured at time of writing. Mixed radiation-matter excitations in semiconductors, also called exciton-polaritons, have been shown to be sensitive at the single quantum level munoz-matutano_emergence_2019 ; delteil_towards_2019 , which might play a role in realizing analog simulators of strongly interacting photonic lattices gerace_quantum-optical_2009 ; carusotto_fermionized_2009 , but their effective use as qubits is still immature. On the other hand, a few companies are investing in a photonic-based quantum computer, which could then be used as a UQS employing continuous variable cluster states menicucci_universal_2006 , but we are not aware of any proof-of-principle demonstration at the moment.

Magnetic molecules manipulated through electromagnetic pulses have also been proposed as a potential platform for quantum information processing santini_molecular_2011 , thanks to their long coherence times and high degree of chemical tunability. This allows to engineer suitable structures of elecronic ferrando-soria_modular_2016 or nuclear atzori_two-qubit_2018 spin qubits in which the qubit-qubit interaction is effectively switched on and off by electromagnetic pulses. Furthermore, the richness of the molecular Hilbert space can be exploited to directly encode logical qubits with embedded quantum error correction in single molecules hussain_coherent_2018 .

Finally, also arrays of optically or magnetically trapped neutral atoms, typically excited in Rydberg states, have been proposed as a platform for quantum simulations weimer_rydberg_2010 , although several challenges need to be addressed for quantum computing applications, particularly related to the short coherence times as compared to trapped ions saffman_rydberg_2016 . While analog simulators with more than 50 Rydberg atoms have been reported bernien_probing_2017 , and two-qubit entanglement has been recently shown with $\sim 97\%$ fidelity levine_rydberg-qubits_2018 , it is still premature to expect digital quantum simulators based on fully controlled Rydberg-atomic qubits.

Together with existing technologies that are moving their first steps into actual quantum computing applications, it is worth concluding this brief overview by mentioning a few potentially promising hybrid technologies, which are usually aimed at merging the best characteristics of two or more existing approaches xiang_hybrid_2013 ; kurizki_quantum_2015 . The philosophy behind these proposals is simple: it is quite likely that a hybrid technology will be in the best position to simultaneously meet all the requirements in terms of scalability (possibly in multi-dimensional arrays), chip-scale integration, and high operational reliability (i.e., long qubits coherence and short gating times). A number of proposals for prospective quantum technologies have been reported, which we will hereby summarize briefly and refer to the original references for further details.

For instance, spin ensembles coherently coupled to superconducting microwave resonators have been proposed as a backbone of a novel hybrid quantum technology ping_measurement_2012 ; carretta_quantum_2013 ; chiesa_robustness_2014 . This hybrid architecture would exploit the long coherence times of spin ensembles and the easy manipulation of photons in tunable resonators. A full digital quantum computing architecture has been devised chiesa_digital_2015 ; chiesa_long-lasting_2016 , for which we report an example in Fig. 12(a), where the transverse field Ising model of 3 spins (see Hamiltonian model in the inset) is theoretically shown to be simulated with a large overall fidelity ($\sim 95\%$ on average) when realistic dissipation parameters are assumed chiesa_digital_2015 . This kind of architecture could be even built with single magnetic molecules strongly coupled to the quantized resonator field jenkins_coupling_2013 ; jenkins_scalable_2016 .

Along similar lines, hybrid architectures based on Nitrogen Vacancy (NV) centers coupled to Carbon nanotubes have also been proposed li_hybrid_2016 . Mechanical degrees of freedom have also been considered to be part of hybrid platforms, due to their intrinsically low dephasing rates. In particular, quantum information processing has been theoretically shown in optical devices in which qubits are encoded in the lowest lying mechanical levels rips_quantum_2013 , as well as in NV centers coupled to mechanical resonators and superconducting waveguides li_hybrid_2015 . Recently, mechanical qubits encoding has been considered in a hybrid set up coupling vibrating nanoresonators to superconducting circuits tacchino_electromechanical_2018 . An example of such an electromechanical quantum computing architecture is reported in Fig. 12(b), with theoretical simulations of the digitized evolution of the transverse field Ising model of 2 spins (see Hamiltonian model in the inset) performed on such hypothetical platform, with very interesting fidelities in the order of $99.9\%$ for the overall quantum simulation if the EMR dephasing is neglected, which reduces to about $99\%$ for realistic pure dephasing rates. These numbers are extremely promising, and could motivate further experimental efforts towards realization of the required building blocks.

V Outlook and perspectives

We have given a brief summary of the current status on quantum simulators, restricting our overview to the use of quantum computers as general purpose machines that can be programmed to solve for the exact time evolution of an arbitrary Hamiltonian model. The only prescription for such a universal quantum simulator is that the physical model under analysis be mapped onto an effective local Hamiltonian obeying the algebra of Pauli matrices, which is then encoded directly in a qubit-based quantum computer through a quantum circuit model. A number of observables, such as spectra and correlation functions, can be accessed upon measurement in the computational basis. After giving a pedagogic introduction to the theoretical background allowing to translate the digitized unitary evolution in discrete time steps into the corresponding quantum algorithm made of a sequence of one- and two-qubits gates, we have reviewed recent experimental results on the two leading quantum technology platforms. Finally, we have outlined a few existing technologies that might develop into quantum computing hardware, and hence be useful for quantum simulation, as well as prospective hybrid approaches that might eventually be tested.

Since most Hamiltonian models of physical interest can be expressed in terms of locally interacting spin terms, we have focused this review on the most widespread spin-type models, such as the Heisenberg and Ising models in an external magnetic field. Getting acquainted and applying the basic techniques developed for such models allows to quickly grasp the quantum simulation of more general quantum many body systems that are typically intractable with classical simulations, mostly due to exponential scaling of the required resources with the system size, such as the Fermi-Hubbard model. The road to quantum advantage is an exciting targeted goal to be fulfilled in the coming years, following the availability of NISQ quantum hardware with few tens of non-error corrected qubits. While the advent of quantum error correction will most probably have a transformative impact allowing to realize universal quantum simulations with arbitrary digital precision, in the intervening years error mitigation techniques and further technological improvement will bring interesting results also from current experiments, which are limited by noisy gates and qubits coherence times. In this respect, increasing gate fidelities and reducing gate duration are some of the technological challenges to be faced in the near term. In fact, these developments should all lead to an increase in “quantum volume” moll_quantum_2018 ; cross_validating_2018 ; wei_verifying_2019 , i.e. a larger number of actual qubits usefully participating in a given quantum computation, which is currently limited to less than 10 in essentially all of the available platforms. When such a number will actually be on the order of 30 to 40, quantum advantage will finally be within reach of such NISQ digital quantum simulators, at least for some targeted applications or models.

In the meantime, a great deal of work is ongoing to devise new potential use cases and algorithms to be run on these machines, for which learning techniques of quantum circuit programming might turn being useful. A brief, non-exhaustive list, with a bit of personal taste, is given in the following. Restricting to problems of academic interest, the dynamical localization of quantum Hamiltonians that have a classical chaotic behavior benenti_dynamical_2003 ; montangero_dynamically_2004 , requiring simulation of the quantum Fourier transform, could be run on a universal quantum computer. More recently, universal quantum computers have been receiving attention from machine learning applications, in particular to develop quantum neural networks, with the aim of processing an exponentially large amount of data with polynomial resources biamonte_quantum_2017 . The first attempts in this direction have been reported havlicek_supervised_2019 ; schuld_quantum_2019 ; tacchino_artificial_2019 . Universal quantum simulators might also help solving problems in open quantum system dynamics, for which novel numerical approaches have already been developed finazzi_corner-space_2015 ; mascarenhas_matrix-product-operator_2015 . Simulating the digitized non-unitary evolution of an open quantum system on a quantum computer is a topic of current interest lloyd_engineering_2001 ; wang_quantum_2011 ; muller_simulating_2011 ; schindler_quantum_2013-1 ; di_candia_quantum_2015 ; cleve_efficient_2017 ; lamm_simulation_2018 . Finally, the huge body of knowledge accumulated in the past half a century to classically simulate the many body dynamics of quantum systems of increasing complexity, such as quantum Monte-Carlo, molecular dynamics, and density matrix renormalization group could be integrated into quantum algorithms to be run on digital quantum computers, with far-reaching and still unknown consequences.

VI Acknowledgments

This review is based on a series of lectures on Universal Quantum Simulators DG has given at the Ph.D. school in Physics, University of Pavia, in January 2019. A number of collaborators and friends are acknowledged for useful discussions and inputs over the years: G. Amoretti, L. C. Andreani, A. Auffèves, D. Bajoni, P. Barkoutsos, G. Benenti, S.-Y. Chang, M. Grossi, M. D. LaHaye, C. Macchiavello, S. Montangero, P. Santini, M. F. Santos, V. Savona, I. Tavernelli, F. Troiani, H. E. Tureci. This work was partly funded by the Italian Ministry of Education and Research (MIUR) through PRIN Project 2015 HYFSRT “Quantum Coherence in Nanostructures of Molecular Spin Qubits,” PRIN Project 2017 INPhoPOL “Interacting Photons in Polariton Circuits”, and by EU H2020 funded QuantERA ERA-NET Cofund in Quantum Technologies through projects CUSPIDOR and SUMO, both cofunded by Italian MIUR.