A quantum computing system and a method for using such a quantum computing system

Abstract

A computer-implemented method is disclosed herein which comprises: providing a Hamiltonian representation of a physical system; using a first quantum circuit to block encode the Hamiltonian into a matrix; and using a second quantum circuit to perform a Quantum Singular Value Transformation to invert the matrix. The inverted matrix is used to provide the Green's function corresponding to the physical system for determining one or more physical properties of the physical system.

Description

A QUANTUM COMPUTING SYSTEM AND A METHOD FOR USING SUCH A QUANTUM COMPUTING SYSTEM Technical field This patent specification relates to quantum computing systems and methods for using such quantum computing systems to help determine physical properties of physical systems. Background Quantum computers are inherently suitable for determining properties of complex physical systems, such as modelling the energy states of molecules, because they exploit quantum phenomena. Unlike classical digital computers in which the basic unit of computation, the bit, has in one of two discrete binary states (1 or 0), quantum computers utilise qubits which may exist in a superposition of many different computational states. This allows a quantum computer to investigate many different states in parallel. However, quantum computers are vulnerable to noise which may lead to decoherence between the different quantum states. This noise problem becomes more significant as the number of qubits increases for more complex computational processing. There is interest in developing procedures and implementations of quantum computers for determining the physical (e.g. spectroscopic) properties of physical systems such as atoms and molecules. Summary Provided are computer-implemented methods and computer systems for use when analysing physical properties of quantum systems, such as optical and electronic properties and behaviours including quantum state changes. The described methods and systems are useful, for example, when simulating a molecular system - which may be an individual molecule or a solid material. The invention is defined in the appended claims. A computer-implemented method is disclosed herein which comprises: providing a Hamiltonian representation of a physical system; using a first quantum circuit to block encode the Hamiltonian into a matrix; and using a second quantum circuit to perform a Quantum Singular Value Transformation to invert the matrix. The inverted matrix is then used to provide a Green’s functioncorresponding to the physical system, for determining one or more physical properties of the physical system. For example, the inverted matrix may provide matrix elements for calculating a single-particle, many-body Green’s function. The determined properties could be ground state energies or excitation and ionization energies, for example. The input representation of the physical system is preferably provided as a computer-readable expression including one or more operators representing the dynamic response behaviour of the physical quantum system when manipulated by the operators, such as defining a relationship between an external perturbation or internal interaction and the system’s response. For example, a Hamiltonian representation of the energy of a physical quantum system may be provided as an input to a quantum computation system for implementing the invention to calculate the energy of an excited state. For example, the Hamiltonian or other input representation may be provided as an input from a separate computer system or from a storage device. A ground state or other reference state of the quantum system may also be an input. The method of the invention can be applied to different physical systems using input representations of those physical systems that include different excitation (and annihilation or ‘de-excitation’) operators relating to the dynamic response behaviour of the respective physical system. The system and method disclosed herein may also utilise conventional computing systems, for example as an adjunct to a quantum computing system or as a simulation environment for a quantum computing system. These methods may be implemented in a hybrid classical-quantum computer system, including a classical (digital) computing apparatus that generates one or more measurable quantum circuits to block encode the Hamiltonian into a matrix and to perform the Quantum Singular Value Transformation (QSVT) to invert the matrix (that is used to calculate the Single-Particle Green’s Function) and a quantum computing apparatus that executes the quantum circuits on its qubits or qudits and outputs measurement results indicating changes to the quantum states of the physical quantum system. Allocation of appropriate tasks to each of a classical binary digital computer and a coupled quantum computer allows exploitation of the respective capabilities of each. This can provide significant advantages when the quantum circuits are executed on resource-constrained noisy intermediate-scale quantum (NISQ) computers, if using the method as described below. The methods and systems described below enable analysis of the properties and behaviours of real- world physical quantum systems, using computer-readable expressions to describe those propertiesand behaviours, and provides a computationally acceptable mitigation of the errors that result from an approximation. This method enables a reduction in the number of single-qubit and two-qubit gates compared to other techniques in the calculation of single-particle many-body Green’s Functions which takes into account the resource constrained NISQ computers. A quantum system responds to disturbances via the responses of its excited states (which are frequency-dependent, and therefore time-dependent), and Green’s Functions can be used to analyse those responses. Green’s Functions can be used to calculate physical and dynamic properties and behaviours such as: ^ the effects of particle and spin interactions within the system, ^ distribution of excited states, for investigating transitions between states as a function of frequency of the time-dependent perturbation, ^ response functions for analysing the effects of external perturbations and correlations between the applied perturbation and the system response. These measured behaviours can provide information about spectroscopic properties (photon emission, photon absorption) of molecules, energy levels, lifetimes of quantum states and conductivity. In Dynamical Mean Field Theory (DMFT), Green’s Functions describing electronic properties of materials can be used to provide information about particle interactions and propagation within a material, for example describing the behaviour of electrons at a single lattice site of an isolated impurity and using this to approximate electronic properties within a solid material. Therefore, Green’s Functions can be useful to determine the behaviour and evolution of the quantum system as a simulation that avoids the high cost of synthesizing a new material. This is highly advantageous for cost-efficient materials development and in the use of simulation to enable optimised design of a physical system, and for many industrial applications that require an understanding of the quantum behaviour of a physical system. However, calculation of the Green’s function is also known to be computationally expensive, especially as the size of the quantum system of interest scales up (which is useful for solving large quantum many-body systems that are intractable for classical computers). This is a problem as typical NISQ devices are limited in terms of the hardware resources. Previous quantum computing methods that calculate the Green’s function, such as the variational quantum eigensolver (VQE) or quantum subspace expansion (QSE), result in either a large overhead in terms of the number ofmeasurements or require large amounts of resources including lots of ancilla qubits and a circuit depth that is excessive for noisy intermediate scale quantum (NISQ) computers. Therefore, there is a need for a method that is suitable for computing the Green’s functions, which mitigates the problems and limitations outlined above. Described herein are methods and systems that achieve a qubit-efficient block encoding of a Hamiltonian into a matrix with a reduced circuit depth and perform the QSVT algorithm to invert the matrix with a computationally acceptable performance, requiring fewer qubits. A Green’s Function can be evaluated at a plurality of different frequencies which are selected based on the physical properties of the system being evaluated – i.e. frequencies related to the energy of a perturbation such as an incoming photon. For example, analyzing a Green’s Function at different frequency points can provide insights into different physical properties and interactions, such as resonances and energy transfer mechanisms. As each frequency point in the Greens Function calculation requires its own quantum circuit to implement a matrix inversion via the QSVT algorithm, the evaluation for various frequency points can be run across multiple quantum processing units (QPUs) in a hybrid quantum-classical workflow. Brief Description of the Drawings Various implementations of the approach described herein will now be described in detail by way of example with reference to the following drawings: Figure 1 depicts an example of a circuit performing block encoding of a matrix for use in a computing system and / or method as described herein. Figure 2 depicts an example of a circuit template for a Pauli operator for use in a computing system and / or method as described herein. Figure 3 (Table 1)) provides a summary for the scaling of each part of the quantum circuit, including circuit scaling summaries for different unitaries used to implement a block encoding of a linear combination of Pauli operators via the LCU method as an example of the approach described herein. Figure 4A depicts an example circuit construction of a multicontrol single-qubit (^)^Z gate, for ^ = 2, for use in a computing system and / or method as described herein. Figure 4B depicts an example circuit construction of a multicontrol single-qubit (^)^Z gate, for ^ ∈ {1, 3}, for use in a computing system and / or method as described herein. Figure 5 depicts a quantum circuit to implement the QSVT for an odd degree polynomial.Figure 6 is a plot showing the absolute error of a quantum signal processing (QSP) circuit approximating f(a) = 1 / a for two different values of k, namely k=10 (red) and k=50 (blue). (The red line is generally the lower line in the left-hand plot of Figure 6, but generally the upper line in the right-hand plot of Figure 6). Figures 7 and 8 show a spectral function of a single-particle Anderson model as utilised herein for k=10 (Figure 7) and for k=50 (Figure 8), in both cases in comparison with an exact calculation. Figure 9 is a schematic diagram of a metal to insulator Mott transition of the two-site Anderson impurity model. Figures 10-12 and 13-15 provide spectral plots based on examples of the approach described herein applied to the two-site single particle Anderson model at different U, V values. The results of Figures 10-12 are for k=10 and the results of Figures 13-15 are for k=50. Figures 16 and 17 depict a metal to insulator Mott transition of the two site Anderson impurity model. The left plots show the density of states on the Bethe lattice and the right plots are absolute error in the density of states with respect to the exact solution. The results of Figure 16 are for k=10 and the results of Figure 17 are for k=50. Figure 18 shows a sequence of steps of a method according to the invention. Figure 19 is a schematic representation of an example of an example computing system comprising a hybrid combination of a classical digital computing system and a quantum computing system. Figure 20 depicts a quantum circuit to implement the quantum chemistry Coulomb operator with mid circuit measurement. Figure 21 shows a quantum circuit for performing a Hadamard test after a matrix inversion. Detailed description Outline Section I below provides an introduction to the problem space. Section II describes background material that is helpful for understanding the approach described herein. The single particle Green’s function, the quantum singular value transformation (QSVT) algorithm and the single-impurity Anderson model are reviewed. Section III describes a method and system implementing the approach described herein. The results of a numerical study are then presented in Section IV. Section V compares the block encoding strategy described herein with prior publications. Section VI presents some additional comments and description and section VII presents an Appendix.Section 1 – Introduction In physics, many systems are described using linear inhomogeneous differential equations. In the 1820s, George Green developed tools to deal with such problems [1, 2]. His work was rediscovered by William Thomson (later Lord Kelvin) [3]. The mathematical methods Green proposed are widely used in physics now and played a major role in the development of quantum field theory. Green’s function has since become one of the most important tools in many-body theories [7]. The function allows one to calculate many physical properties of a system, for example: excitation and ionization energies, ground-state energies, transition matrix elements, absorption coefficients, and dynamical polarizabilities, as well as elastic and inelastic electron cross sections [8]. In addition, self- consistent perturbation theories can be formulated in terms of the Green’s function [8]. More details on Green’s function theory may be found in standard textbooks on many-body theory [9]. Described below are methods and systems adapted to obtain the single particle many-body Green’s function using an approach specifically designed to operate efficiently on a quantum computer. There have been several different proposals for utilising a quantum device to calculate the matrix elements of the Green’s function. Some of these have focussed on variational approaches [10–16] and other near-term methods

[0017] . Alternatively, other proposals have utilized quantum phase estimation [18, 19]. The present approach has some similarities to Tong et al. in

[0020] , in which the quantum singular value transform (QSVT) algorithm is used

[0021] . However, Tong et al. apply a specific fast inversion strategy which is used to reduce the query complexity to the block encoding by preconditioning the linear problem, whereas the present approach does not utilise the same inversion strategy. Some methods assume that the matrix to block encode H can be split as H = A + B, under the assumption that the spectral norm of A is much greater than B: ∥A∥ >> ∥B∥. By preconditioning the problem using A−1, the new query complexity depends on ∥B∥ rather than ∥H∥, thereby reducing the overall computation cost

[0020] . Described below is a qubit-efficient block encoding technique that avoid the requirement for precfonditioning of

[0020] . The approach described herein has been tested on the two-site single impurity Anderson model (SIAM)

[0022] , which is a four qubit problem. Despite its simplicity, it captures some very interesting physics, for example it can be used to approximate the Mott insulator phasetransition

[0023] . The single-particle Green’s function has been calculated for this system via an implementation of the quantum singular value transform (QSVT) which mitigates approximation errors. The results have been used to plot this phase transition and demonstrate how approximating the inverse function can lead to singular values not being inverted properly, which in turn causes calculation errors due to the matrix not being inverted properly. Described below are specific solutions to mitigate against this problem. The Quantum Singular Value Transformation (QSVT) is a technique that provides a unified framework for describing many of the quantum algorithms discovered to date. We implement a noise-free simulation of the technique to investigate how it can be used to perform matrix inversion, which is an important step in calculating the single-particle Green’s function in the Lehmann representation. Due to the inverse function not being defined at zero, we explore the effect of approximating f(x) = 1 / x with a polynomial. This is carried out by calculating the single-particle Green’s function of the two- site single-impurity Anderson model. We also consider a new circuit construction for a linear combination of unitaries block encoding technique that reduces the number of single and two-qubit gates required. Section 2 – Method for block encoding and matrix inversion In this Section, we summarize mathematical details and notation, review the one-particle Green’s function and the single-impurity Anderson model that defines a type of physical system, and introduce and describe an implementation of the quantum singular value transformation algorithm. A. One-particle Green’s Functions The time-ordered single-particle Green’s function (GF) at zero temperature in the frequency domain is defined in the Lehmann representation as [24, 25]:where:For simplicity we have assumed |Ѱ^^ to be non-degenerate, but this can be extended to degenerate ground-states at nonzero temperature. ^(^) ^^ (^) and ^(^) ^^ (^) are called the advanced and retarded Green’s function respectively, or the electron and hole excitation parts of the GF [19, 20]. Here ^^and ^^ are fermionic creation and annihilation operators of an ele^ ctron in the i-th spin orbital, is the ground-state wavefunction, ^^is the ground-state energy, H is a second quantized fermionic Hamiltonian and z = ω + iδ is a complex frequency often interpreted as an energy shift

[0020] . The imaginary part of z, given by δ, is small and required for convergence of the Fourier transform [8]. Equation 1 can be mapped to an equation involving qubit operators, by applying a fermionic-to-qubit transformation to the fermionic operators. For a given z, ^(^) ∈^ ^is an N x N matrix that isefficient to store classically, where N is the number of spin orbitals (or qubits) describing the system. Equation 1 shows how the i-th row and j-th column of G(z) is calculated. Even though this matrix is efficient to store and manipulate classically, it should be noted that each entry requires solving an exponentially large problem. This is due to the size of the Hamiltonian scaling as ^^^or exponentially with the number of spin orbitals. Classically computing each entry in the Green’s function quickly becomes intractable, as doing so requires inverting an exponentially large matrix. Such a Hilbert space is naturally expressed on a quantum computer with N qubits, but this would require the ability to perform a matrix inversion on such a device. The inventors have identified a way to perform matrix inversion using the quantum singular-value transform algorithm while mitigating errors. As will be discussed in the next section, this method provides a way to apply an (approximate) inverse of a block encoded operator onto a quantum state, enabling the calculation of Green’s Functions on a quantum computer for evaluating properties and dynamic behaviours of quantum systems. A point to note when calculating the Green’s function is that the ground-state |Ѱ^^ must be known. In this work we assume it is known a priori and can be efficiently prepared on a quantum device. However, in general the ground-state problem of a l-local Hamiltonian is QMA-complete for l ≥ 2 (for l = 1 the problem is in P)

[0026] . Currently there are no known algorithms that can find a solution in polynomial time, although solutions may in some cases be attainable in an acceptable time-frame. B. Quantum Singular-Value Transform Algorithm A comprehensive review on quantum signal processing

[0027] and the quantum singular-value transform

[0021] can be found in

[0028] . In this section we summarise the steps involved in performing matrix inversion via QSVT. The algorithm can be broken down into four major steps: 1. Construct a quantum circuit that block encodes a matrix.2. Generate the quantum signal-processing angles required to implement the desired function that will be applied to the singular values of the block encoded matrix. Here this will be an approximation of the inverse function: f (x) ≈ 1 / x. 3. Construct the quantum circuit to implement the QSVT algorithm using the outputs of steps 1 and 2. 4. Implement a Hadamard test to evaluate the real and complex parts of each entry in the Green’s function. Quantum signal processing (QSP) is a known method for exact implementation of matrix polynomials on quantum computers. The quantum signal processing angles are angles that characterize transformation of quantum states. These angles are used to manipulate the quantum states to perform operations such as transforming the information encoded in qubits. In the context of the present invention, quantum signal processing angles are the angles that define the unitary operations applied to the quantum states to perform the singular value decomposition. The following subsections review each of these steps, apart from the Hadamard test since a full analysis of step 4 is given by Tong et al. in

[0020] . 1. Linear Combination of Unitaries (block encoding) There are many different methods available to block encode a matrix [27, 29]. We focus on the linear combination of unitaries (LCU) approach, a technique to block encode any linear combinations of unitary op- erators [30, 31]. Given such a matrix A:where, without loss of generality, we can assume ^^> 0 and ^^∈ ℝ ∀i by absorbing any complex phases and signs into the unitaries ^^[27, 30, 32]. Given a list of ^^and each ^^, which are assumed to be easy to implement as controlled operations on a quantum device, the block encoding can be constructed using the oracles [27, 32]:andIn the following, the subscript s denotes the system register and p the prep (ancilla) register. The number of prep qubits will be np = [log2(k)].The PREP or “Prepare” oracle is a unitary that prepares the state |^^ =||^||^|^^ from the all zero state on the prep register - i.e. |0^^ ^ |P^. This is why in equation 4 only the first column of the PREP unitary is defined. As discussed in

[0032] , the other columns can take any value providing that PREP remains unitary. This means there is a lot of freedom in how to construct this operator. If one simply finds the quantum circuit that realises |0^^^^ |^^^, then the circuit’s action on the other basis states is automatically accounted for and the whole of PREP will be defined

[0032] . The only quantum circuit requirement is being able to generate any real quantum state from the all-zero state on the prep register. There are many different proposals on how to prepare arbitrary quantum states [33–36]. Following the approaches given in both [35, 36], a real quantum state can be generated using multiplexed Ry rotations with the number of single-qubit and CNOT gates both scaling as ^(2^^). As the number of prep qubits scales logarithmically with the number of terms in ^, ^( ^^^^|^|), the number of single-qubit and CNOT gates will scale linearly as ^(|^|)for the PREP part of the block-encoding circuit. The desired LCU block encoding is achieved by performing ^^^^^^^^^^^^^^^^ and post selecting on the all- zero state on the prep qubit register. We can check this via the following proof

[0027] :This is implemented according to the circuit in Figure 1. Figure 1 depicts a circuit performing block encoding of ^ = ∑ , where ∑^|^^| = 1. Here the PREP unitary (equation 4) prepares the following quantum state: |0^^The select operator is then performed (equation 5) followed byPREP†. The overall circuit produces the following final state:+ ^1 − ||^|^^^^^||^ |⊥^. Post selecting the all zero state on the “preparation” register(|0^^^^^^)results in M being applied to the system state[32, 37]. Note the number of qubits needed by the system registeris defined by the physical problem, and the number of preparation qubits required is [^^^^(|^|)]The probability of success for this block encoding is As discussed in

[0032] , A is not necessarily unitary and so A† A may not equal I. The probability of success therefore depends on the system state |ψ^ and the 1-norm of the block-encoded matrix. Amplitude amplification [38–40] and oblivious ampli- tude amplification [41, 42] can then be used to increase the probability of success

[0021] .In the literature, it is common to see (α, κ, ^)-block encodings. Here α is a normalisation factor of the block- encoded matrix, κ is the number of extra ancillary qubits to implement the block encoding and ^ is the error of the block encoding. In this work, the LCU is given as a linear combination of Pauli operators. The “SELECT” oracle (equation 5) applies a controlled version of each of these Pauli operators on the system register, controlled by the prep register. This involves performing multi-control Pauli operators with phases{i, −i, 1, −1}. Following the work in

[0032] , this can be achieved using the template given in Figure 2. Figure 2 depicts a circuit template to perform as an n-control Pauli operator having single-qubit Pauli Z and I matrices

[0032] . Figure 4B (see below) provides the circuit construction for the multi-control (^)^^ gate. The single-qubit gates:{^, ^^, ^}may be used to convert this circuit into a general n-control (^)^^ operator via a change of basis. The relevant phases are then obtained via the following identities:These can be implemented according to the circuit templates summarised in Figures 4A and 4B. Figures 4A and 4B depicts circuit constructions of a multicontrol single-qubit (^)^Z gate, for k ∈ {0, 1, 2, 3}. In particular, Figure 4A shows the circuit for k=2, while Figure 4B shows the circuit for k ∈ {1, 3}. A circuit for k = 0 is not explicitly pictured, as it is the trivial case of a multicontrol single-qubit Z gate. This is the same as Figure 4A with both X gates removed. By performing a change of basis on certain qubits, using the single-qubit gates {S, S†, H}, the circuit proposed in Figure 2 can be used to generate any multicontrol Pauli operator with a ±1, ±i phase. The ordering of unitaries in equation 5 above is arbitrary, but an optimal ordering may lead to significant circuit simplifications. The work in

[0032] ,

[0043] and

[0044] can readily be applied to this problem. To determine the overall circuit cost to implement ^^^^^^^(equation 5 and Figure 1), we determine the cost of implementing a multicontrol Z gate and a multicontrol Rz gate. Following the proposal by

[0045] , any n- control single-qubit gate can be decomposed with ^(n2) single qubit and CNOT gates with linear depth

[0045] . For an n-control single qubit Z gate with n ≤ 6, the approach outlined in

[0046] and

[0047] (theorem 8) uses fewertwo-qubit gates, where the number of single-qubit and CNOT gates used scales as ^(2n). In general, following

[0045] , the cost of performing a multi-control Pauli operator via the template in Figure 2 scales as: 1. O(2ns) single-qubit gates, required to implement a change of basis. 2. O(2[ns− 1]) CNOT gates, performing the ladder of CNOT gates on the system register. 3. O(n2c) CNOT and single-qubit gates for the multicontrol ikZ gate. Here ncis the number of control qubits and nsis the number of ‘system’ qubits the Pauli operator acts on. The single-qubit and CNOT gate cost per nc-controlled Pauli operator scales linearly in system qubits and quadratically in control qubits as O(n2c + ns). The overall cost of implementing ^^^^^^^via the circuits presented will depend on the number of nc-controlled Pauli operators that are used. Looking at Equation 5, it can be seen that k operators are used to bring the final cost to O(k[n2c + ns]) single and two-qubit gates. As k = |A| and nc = ⌈log2(|A|)⌉, the final scaling of single-qubit and CNOT gates can be written as ^(|A|(⌈log2(|A|)⌉2+ns). Table 1 (see Figure 3) provides a summary for the scaling of each part of the circuit. In particular, Table 1 shows circuit scaling summaries for different unitaries used to implement a block encoding of a linear combination of Pauli operators via the LCU method. Here |A| denotes the number of Pauli operators in A (see Equation 3 above), ^^denotes the number of control qubits and ^^denotes the number of system qubits. The SIAM Hamiltonian considered herein is constructed as a linear combination of Pauli operators (see equation 13 below). Each individual Pauli operator Pi is a unitary Hermitian operator that is easy to implement on a quantum computer as a controlled operation. However, rather than block encoding H, we block encode the following complex shifted Hamiltonians:The reason why the Hermitian conjugate is taken in equations 8a and 8b stems from how a matrix inverse can be obtained from a singular-value decomposition. If we write vectordecomposition of an arbitrary matrix as A then its inverse (if it exists) is Takingthe Hermitian conjugate of A yields ^^= VΣ^^and to find the inverse of A all that remains is to invert the singular values of ^^. The QSVT algorithm approximately inverts the values of a block encoded matrix. Performing this algorithm on the block encodings of ^(^)and ^(^)will therefore allow us to calculate the inverseparts of equations 2a and 2b. Each inverse may be calculated separately to reduce the circuit depth; however, it is possible to calculate the terms simultaneously using the method of adding different block encodings given in [21, 48]. Although this is a viable approach for determining an inverse which may be adopted in some implementations of the approach described herein, other implementations as described herein do not adopt this approach so as to be more amenable to early fault-tolerant quantum computers. As z is a constant complex shift (has a real and imaginary component), ^(^)and ^(^)are no longer Hermitian operators, whereas H is. The Quantum Eigenvalue Transformation (QET) implements a polynomial transformation of a block-encoded Hermitian matrix when the polynomial of interest is represented by QSP

[0027] . For a polynomial transformation of a general matrix, the quantum singular- value transformation is used

[0021] (and is also adopted herein). It is noted that for a Hermitian matrix with a block-encoding input model, the quantum circuits of QET and QSVT can be the same

[0049] . It is possible to convert the non-Hermitian problem into a Hermitian one via matrix dilation, which requires an extra qubit - for further details see Section 4.2 in

[0050] . This matrix dilation is not necessary for the present approach as described herein, it may be used if the quantum linear-system algorithm proposed by Harrow, Hassidim, and Lloyd (HHL) were to be used [11, 51]. As noted in

[0028] on block encodings, generally additional work is involved to find block encoding techniques optimal to a given physical system. The LCU method is general and so can be utilised in such situation, but does not utilise any underlying structure of a problem and requires ⌈log2(|A|)⌉ extra ancillary qubits. It is possible to make a (α, 1, 0) block encoding of any n-qubit matrix A - i.e. only involving a single ancilla qubit for a block encoding (see example 6.2 in

[0052] and appendix D in

[0053] ). These schemes require a singular-value decomposition of A and thus will not scale in a general setting. However, certain structures in particular physical problems may allow for more efficient encoding strategies. 2. Matrix inversion via quantum signal processing To perform the QSVT, the quantum signal processing angles [27, 37, 54] are generated to implement a function (or usually some approximation of a desired function). On a single qubit, QSP is usually defined as

[0054] :where a ∈ [−1, 1], W(a) = Rx [2 ^^^^^(a)] and P is a polynomial with degree at most the length of the sequence of QSP phases (≤ d). The constraints on what sort of polynomials can be implemented using this technique are covered in

[0054] . In equation 9, once the polynomial to be implemented is fixed and the QSP angles are defined, all the {^^|k = 0, 1, ..., d} ∈ ^^remain fixed. The only free variable remaining is a. It is therefore always possible to plot P(a) by calculating: ^0| Uϕ (a) |0^ = P(a), where one scans over −1 ≤ a ≤ 1. The determination of the angles in ^^can be regarded as a ‘black-box’ calculation, see [28, 55–57] for further details. Once the phases ^^have been calculated for a particular polynomial, thephasescan be reused and do not need to be calculated again. The pyqsp

[0058] and QSPPACK

[0059] open-source libraries allow users to generate different sequences of QSP angles for many different functions. An algorithm proposed by Haah in

[0055] gives a rigorous analysis of how to find the angle sequence corresponding to a supplied polynomial that has a runtime scaling as O(d3polylog(d / ^)), for a degree-d polynomial. This algorithm returns a set of QSP angles for a uniform ^-approximating polynomial over the interval [−1, 1]. The approach described herein uses animplementation of the inverse function involving 1 / x, which leads to consideration of the discontinuity at x = 0. Instead of approximating 1 / x over the full range, the present ^ ^ approach approximates it over the range [−1, − ^] ∪ [ ^ , 1]. The existence of such an odd polynomial is guaranteed in Corollary 69 of

[0060] . Importantly, the approximation of 1 / x used in QSVT requires all the singular values of the block-encoded matrix {σ} to be σ ≥ 1 / κ ∀σ, otherwise they fall into the region where the polynomial approximation of 1 / x is ill-defined. Extending the single-qubit QSP (equation 9) to higher dimensions is discussed in

[0021] (see theorem 2), where ideas from qubitization

[0027] and two-dimensional invariant subspaces coming from Camille Jordan’s Lemma

[0061] are used. Their results show how to apply certain polynomials to a block encoded matrix:Here A is written using its singular-value decomposition. The location of A in B is determined by certain projectors ∏^

[0028] , herein: |0^^^0^|. The QSVT circuit, for odd d, can be built as [21, 28, 60]:wh Note the QSP phases ^^∈ ℝ^^^(equation 9), have been modified to ^^ ^∈ for QSVT. This accountsfor W(a) not being a reflection operator, which is better suited to the qubitization formalism

[0027] . In summary, equation 11 shows how a polynomial transform is applied to the singular values {^^} of A (equation 10). The above analysis assumes A to be a square matrix, but this is not a requirement

[0028] . Figure 5 depicts a quantum circuit to implement the QSVT for an odd degree polynomial (equation 11)

[0021] . Here the polynomial approximation of ^(^)≈ 1 / ^ has odd degree. The zero controlled NOT gates have zero controls on all the ‘prep’ qubits. As the polynomial approximation of the inverse function is real, the presented approach may use the {|+〉,|−〉} signal basis on the QSP qubit

[0028] . This is why the Hadamard gates (H) are present in the circuit. Each block encoding unitary B is constructed according to the approach outlined in Figure 1. Accordingly, Figure 5 summarises the QSVT circuit, where ∏^(φ) is given by a multi zero-controlled X gate targeted on the QSP qubit and controlled by the ‘prep’ qubits (see Figure 1 in

[0057] ), followed by an Rz rotation on the QSP qubit and followed by another multi zero-controlled X gate: ^− C. Single-impurity Anderson model A common model used to describe strongly correlated electron systems in thermodynamic equilibrium is the Hubbard Hamiltonian. However, classical simulation of this model is severely limited by how many fermionic orbitals can be described, due to the exponential increase of the Hilbert space. Dynamical mean field theory (DMFT) has been developed to address this issue, where the physics of a many-body problem is captured via a single impurity that is coupled self-consistently to a fermionic host (bath)

[0062] . In the limit of alattice with infinite dimensions, for the Hubbard model with infinite coordination number (nearest neighbours), DMFT exactly maps the solution of the Hubbard model to that of the Anderson impurity model. This is because interacting electrons in the Hubbard model in the thermodynamic limit (infinite lattice sites) are modelled by a single-impurity site coupled to an electronic bath (infinite bath sites) that tunnel into the impurity site [17, 62, 63]. DMFT is derived in the limit of infinite lattice coordination; however, for finite dimensions it may still provide good approximations and allow interesting phenomena to be explored [64, 65]. For the approach described herein, we consider a two-site one-dimensional single-impurity Anderson model defined by the Hamiltonian

[0066] :Details on this Hamiltonian are provided in [10, 23, 66]. Note that H is written under the Jordan- Wigner transformation, which allows the fermionic operators to be mapped to spin operators acting on qubits

[0067] . Qubit index 1 (3) represents the impurity spin-up (spin-down) site and index 2 (4) represents the spin up (spin down) bath site. Here, U is the onsite Coulomb repulsion, µ is the chemical potential that controls the electron filling in the grand canonical ensemble, ɛ^describes the on-site energy of the non-interacting bath site 2, and V is the interaction of this bath site with the impurity. (The grand canonical ensemble is a generalization of the canonical ensemble that represents the possible states of a mechanical system in thermal equilibrium with a heat bath at a fixed temperature, where the restriction to a definite number of particles is removed. An example of this is in chemistry, where the number of each molecular species is not conserved but the number of atoms is. For example: 4A + 2B → A4B2, where there are six molecules (par^cles) on the le^ and only one on the right, but always six atoms. Returning to equation 13, H is equivalent for spin-up and spin-down electrons and so the self-energy of the impurity only needs to be calculated for one spin site

[0010] . To solve equation 13 via DMFT, i.e. to find the parameters of the effective model, we consider the Green’s function of the lattice problem ^^^^(^) and the impurity ^^^^(^).For infinite bath sites ^^^^(^) = ^^^^(^). In practice, only a finite number of bath sites can be used and so the difference between ^^^^(^) and ^^^^(^) is minimised.For the present approach, we consider 2-site DMFT under the particle-hole (ph) symmetric case, ^ where µ = ^ and ɛ^= 0

[0010] . The only impurity parameter is therefore V. For a fixed U and given threshold ζ, the following steps are taken

[0066] : 1. For a fixed U, guess an initial on-site energy V, thus determining H (equation 13). 2. Calculate the Green’s function of the Hamiltonian ^^^(^) (equation 1). (a) Herein, each element of ^^^(^) is determined by the quantum singular-value transform. 3. From ^^^(^) define ^^^^(^) (a) This is achieved by selecting the elements of ^^^(^) that correspond to the impurity site. ^^ 4. Calculate the quasi-particle weight ^^^= (1 .(a) The self energy can be obtained as: ∑ ^ ^^^(^) = ^^^^(^)^^− ^^^^(^)^^, where the noninteracting Green’s function is defined as ^ ^ ^^^(^) = (z − ^^+ µ −(b) Due to particle hole symmetry, Im[∑^^^(^^) ] is a single number due to spin-up and - down self-energies being the same for the impurity site. 5. Set6. If |^^^^− ^| ≤ ζ then the bath parameter (and so DMFT) has converged. Otherwise, set V = ^^^^and repeat from step 2. For the 2-site model considered here, there is an analytic form for V [10, 23]:In the present approach, rather than optimising for V at different fixed U, equation 14 is used to determine the optimal V before calculating the Green’s function. The goal is to investigate calculation of the Green’s function via QSVT (rather than performing DMFT self-consistent optimisation loops). Described below is a new block encoding method which implements the present invention for encoding an operator with mid-circuit measurement. The quantum chemistry Coulomb operator has been used as an example. This implementation provides a reduced circuit depth compared to previously known methods, putting rigorous bounds on the complexity of the circuit, and is particularly suitable for Quantinuum’s high-qubit fidelity quantum computer hardware (implemented using trapped ions, with connectivity between non-adjacent qubits).The rank 4 Coulomb tensor containing the electron interaction coefficients for a molecular system can be expressed in second quantisation as follows:This is very expensive (in a with the linear combination of unitaries block encoding method, due to the ^^number of terms (where ^ is the number of spin orbitals). To mitigate this problem, a double factorised form of operator can be used:In an advantageous implementation of the invention, the inventors take this factorised form of the Coulomb operator (or another operator) and implement the internal multiplication of the ^∑^^^ ^ ^^^ ^ ^^^^ in a circuit with mid circuit measurement, as shown in Figure 20,In Figure 20, the U are the unitary Givens rotation matrices, P are the PREPARE oracles and Sel are the SELECT oracles. This type of operation is only possible on devices which enable mid circuit measurement such as Quantinuum’s quantum computing hardware. Assuming the SELECT Oracle is the dominant cost, the complexity of the circuit reduces to: for a small error ^ from the level of truncation in sum over t. Therefore, this reduces the asymptotic cost of implementing the Coulomb operator from ^^to ^ which greatly reduces the circuit depth.After performing the matrix inversion using QSVT, the inverted matrix may be evaluated with a Hadamard test. Figure 21 shows a quantum circuit in which a Hadamard test is sandwiched between controlled Pauli operators (for the fermionic operators). This can be performed after a Jordan-Wigner transform. Section 3 – Implementation and Performance Evaluation We numerically investigated the performance of calculating the Green’s function for the two-site Anderson model via the QSVT algorithm. To build the qubit Hamiltonian, the InQuanto package from Quantinuum was utilized [68, 69]. The circuits for performing QSVT were then constructed using PyTket

[0070] (see also https: / / cqcl.github.io / tket / pytket / api / ). The QSVT circuit was specifically built to perform matrix inversion. First, the QSP phase angles were generated using the open-source python library QSPPACK

[0059] . Next, for each z, we built two quantum circuits that performed T ≈ (z − [H − E0])−1and W ≈ (z + [H − E0])−1via the quantum singular value transform algorithm – see Figure 5. This involved using the block-encoding circuits for (∥B(e)∥1, 3, 0) and (∥C(h)∥1, 3, 0), except for the SIAM Hamiltonian defined for U = 8 and V = 0, where block-encoding circuits for (∥B(e)∥1, 1, 0) and (∥C(h)∥1, 1, 0) were used (due to certain Pauli operators having a coefficient of zero). Each block encoding was constructed according to the template in Figure 1. Note that QSP only approximates the true inverse function via a polynomial, hence the above use of ‘approximately equal’ (≈). In all instances, the complex part of z = ω + iδ was fixed to be δ = 0.1; this was chosen to ensure all the singular values of each block encoded matrix were above 0.02. (These parameters may be varied according to the circumstances of any given implementation). After each quantum circuit was built, a noise-free classical simulation was performed giving the unitary of the whole QSVT circuit for each z value. We post-selected into the correct block of the unitary (see equation 11) to obtain the transformed matrix. For each pair of quantum circuits, the post-selected matrices are T and W. The procedure then classically determines Gij(z) by evaluating〈^^^^^^ ^Ψ^ ^^ ^ ^^^^Ψ^〉(see equations 2a and 2b above) for all i, j, where i and j run over all qubit indices using the standard linear algebra python libraries [71, 72]. The ground state |Ψ0^ used in each calculation was obtained by diagonalizing H on a classical computer for particular (U, V) parameterizations. For each QSVT simulation, the exact classical solution was also calculated, where the Green’s function was calculated via matrix inversion performed on classical hardware.Section 4 – Numerical study The QSVT applies a function, defined by the classically pre-computed^ɸ^angles, to the singular values of a (block encoded) matrix. The user does not have access to the singular values; rather the algorithm applies a function to the singular values: ^ = ^∑^^^ ^^(∑)^^. However, not knowing what the singular values may havematrix inversion via the QSVT. As discussed in Section II B 2 above, the inverse function is not defined at x = 0 and so is ^ ^ approximated over the domain ^−1, − ^^ ∪ ^ 1^. If any singular value of the matrix to invert fallsoutside of this range it will not be This leads to a dilemma, where one needs to know the singular values to determine an appropriate k; however, knowing the singular values is the same as solving the inversion problem. This issue can be resolved in at least two ways. A first strategy makes the approximation of 1 / x arbitrarily small, by using a very large value of k. This is somewhat similar to how conventional computers perform mathematical operations to machine precision. However, the degree of polynomial approximating the inverse function scales as

[0028] : Accordingly, as a better polynomial approximation is used (higher k value), the greater the degree of the resulting polynomial will be. As the circuit depth of QSVT scales as O(d), repeats of the block encoded circuit

[0028] using an arbitrarily large k could (unnecessarily) increase the circuit depth of a given problem. The second approach to determine a valid k, is to estimate the magnitude of the lowest singular value of a matrix ^ ∈ ℂ^ ^ ^. In

[0077] , it is shown that:which provides a lower bound on the magnitude of the value ^^^^of a non- singular n´n complex matrix M. Here, ^ = | det(^)| .||^||^ is the Frobeniusnorm. Evaluating equation 16, allows k to be determined as: ^ ≥ 1 / ^^^^ . However, this approach requires the determinant of the matrix to be found, which can be costly. A further approximation could be used to estimate |det(M)|, such as using the methods in [78, 79].The two-site SIAM considered in this work is defined on four qubits and so classically performing a singular value decomposition (SVD) of a (16 х 16) matrix was possible. Therefore, in order to find appropriate k, rather than using equation 16, we used the true ^^^^. We found all the singular values were above 0.02 and so used a k = 50 approximation. This represents a scenario where all the singular values would be inverted properly via the QSVT. We also simulated a k = 10 polynomial approximation, where some of the singular values lay below 0.1 and so would not be inverted properly. The k = 10 and k = 50 polynomial approximations of 1 / x were represented by 303 and 1519 degree polynomials respectively. Figure 6 illustrates the errors of these polynomial approximations compared to the true inverse function. This figure shows two plots of the absolute error of the quantum signal processing (QSP) circuit approximating f(a) = 1 / a. Note the change in scale of both the x-axis and the y-axis between the left and right plots. The left-hand diagram has the x-axis representing the value of a across the range ±1, whereas for the right-hand diagram the x-axis represents the value of a across the range ±0.1. The scaling of the y-axis changes between the left and right diagrams to allow the latter to show the larger errors closer to a= 0. ^ Each polynomial approximation is defined over the range ^−1, −∪ ^ ^ , 1^. The blue and red data represent different approximations, κ = 10 and κ = 50, of the inverse function(the red line is generally lower than the blue line in the left-hand plot and generally higher than the blue line in the right-hand plot). These approximations are based on 303 and 1519 degree polynomials for the red and blue lines respectively. The left figure is plotted over the domain^ dashed data goes over the range ^−1, − and the solid data over the range ^ ^ , 1^. The constant multipliers given in the Figure legend account for the normalization factors required by each polynomial to ensure they lie between ±1. Note no data is calculated at a = 0. For each value of k, the single particle Green’s function was calculated at different z for different U, V parameters via QSVT. The outputs were then used to plot the spectral function which is defined as:The results for U = 4 and V = 0.745 are given in Figures 7 and 8. The results for the other U, V regimes are given in Appendix A. As this problem is defined over 4 qubits, exact diagonalization solutions were possible to compute classically. In Figures 7, 8 and those in Appendix A, the lowest singular value of the matrix to undergo QSVT matrix inversion is provided. As expected, whenever the singular value lies below 1 / k, the error in the spectral function becomes large. This can be seen for k = 10, wherein Figure 7 errors due to σ being below 0.1 can sometimes lead to a difference from the real answer by ∆^(ѡ)≈ 1. The top portion of Figure 7 depicts a spectral function (corresponding to equation 17) of the single- particle Anderson model for U = 4 and V = 0.745. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 10. The middle portion of Figure 7 shows a plot of the absolute error of spectral function ∆^(⍵) =The bottom portion of Figure 7 shows the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular value below the horizontal black line is outside the region in which the approximation of 1 / x is well defined. The top portion of Figure 8 depicts a spectral function (corresponding to equation 17) of the single particle Anderson model for U = 4 and V = 0.745. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 50 (rather than with k=10 as per Figure 7). The middle portion of Figure 8 shows a plot of the absolute error of spectral function ∆^(⍵) =− ^^^^^(⍵)^. The bottom portion of Figure 8 shows the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular value below the horizontal black line is outside the region of where the approximation of 1 / x is well defined. In Figure 8 all singular values lie above 0.02 (1 / k) and the error in A(ω) remains at ∆^(ѡ)≈ 10^^^; in contrast for Figure 7, some of the singular values lie below 0.1 (1 / k) leading to much higher errors (see the plot in the middle portion of Figure 7 compared to the plot in the middle portion of Figure 8). With reference to the Green’s function definition in equation 1, it can be seen that when the real part of z is equal to ±^^∓ ^^(where ±^^is an eigenvalue of H), then the real part of the denominator in either the advanced or retarded Green’s function vanishes leading to a so-called ‘pole’. For a non- interacting system, where the eigenfunctions are represented by single-configuration states, theamplitudes and Lehmann energies ω are equal to the eigenfunctions and eigenvalues of the corresponding one-electron Hamiltonian [8]. The spectral function consists of a set of peaks at those eigenvalues and each peak is associated with a particle [8, 80]. When interactions are considered, the eigenfunctions are no longer single-configuration states, instead they are in general normalised linear combinations of them. There will now be more non-vanishing contributions to the spectral function, by merging these contributions they will form a structure which can be thought of as de- riving from peaks when the interactions are turned off [8]. This provides a particle-like picture, whereby each peak is now associated with a “quasi-particle” [8, 80]. The spreading of the peak contains information about many-body correlation effects in the interacting system 18]. The spectral function is usually peaked at each energy ^^= ^^+ ^^∆∑^(^^), with a lifetime given byΣ is the self-energy operator. Further details on this are discussed in

[0180] . From the calculated Green’s functions, the Mott phase transition can be plotted by presenting the density of states on the Bethe lattice

[0181] for different U, V regimes at different z. The k = 50 result is summarised in Figure 9, which shows the metal to insulator Mott transition of the two- site Anderson impurity model (error details and the k = 10 result are presented in the Appendix B). The plots presented in Figure 9 show the density of states on the Bethe lattice. The first (top) plot of Figure 9 gives the non-interacting system, followed by four further plots for respective (U, V) combinations of (2, 0.943), (4, 0.745), (5.99, 0.058), and (8, 0). Note that the free density of ^ states on the Bethe lattice with infinite coordination is ^^(^)=^^^^√4^^− ^^(where for present purposes we take t = 1) and the interacting density of states is ^(⍵)= The non-interacting Green’s function is defined as ^(where for present purposes we take ^^= 0 and µ = U / 2) [110, 66]. For the final (bottom) plot (8, 0) of Figure 9, particle hole symmetry is broken and the spin-up and -down parts of∑^^^[^] were treated separately and combined. Figure 9 shows that the QSVT result matches well with the true (exact diagonalization) results. The error of the QSVT approach for the k = 50 polynomial approximation, for each data point, is around the 10^^^level. In contrast, much larger errors are obtained for the k = 10 approximation, ranging from ≈10^^^to ≈ 10^^. This is due to errors occurring when the singular values are below 0.1, which are then not inverted properly leading to errors in the calculated single particle Green’s function and thus the spectral function.Section 5 - Improved LCU circuits In Section 2 B 1, the LCU technique from [30, 37] was reviewed as a potentially efficient way to construct the “SELECT” operator when defined as a linear combination of Pauli operators (see also [32, 45, 46], together with Table I (see Figure 3) which summaries the different costs). The overall cost to implement a LCU block encoding via this approach has a circuit cost scaling as:single-qubit gates. Here ns is the number of system qubits. For situations in which the Hamiltonian is defined as a linear combination of unitaries, the circuit implementation of LCU presented herein will result in quantum circuits with exponentially fewer gates compared with some other techniques, such as FABLE

[0086] , given the exponentially fewer ancillary qubits used to perform the block encoding. For example,

[0086] applies the FABLEapproach to different Hubbard Hamiltonians; in this scenario, the present approach will user exponentially fewer single- and two-qubit gates. Section 6 – Additional Comments One example of the approach described herein is illustrated in the flowchart of Figure 18 and provides a computer-implemented method comprising providing a Hamiltonian representation of a physical system 310. The Hamiltonian used mostly in the current context was the two site Anderson model (in part for the ability to run on a simulator). However, a wide range of other choices are available for the Hamiltonian, such as molecular systems (e.g. drug molecules), battery materials, etc. The method further comprises using a first quantum circuit to block encode the Hamiltonian into a matrix 320. In other words, the Hamiltonian provided at operation 310 on the classical computing system 210 is converted (translated) into a corresponding quantum circuit 260 on the quantum computing system 250 (or simulation thereof). One issue to address here is that although both the Hamiltonian and the quantum circuit 260 can be represented by matrices, the Hamiltonian matrix is generally not unitary, whereas the matrix for the quantum circuit 260 has to be unitary. To overcome this issue, the Hamiltonian is loaded into the quantum circuit 260 by block encoding. More particularly, the Hamiltonian is encoded into the quantum circuit as a block which forms only a portion of the quantum circuit, and the remaining portion of the quantum circuit is configured toensure that the overall quantum circuit is unitary. This encoding is performed using a linear combination of unitaries (LCU), The method further comprises using 330 a second quantum circuit to perform a Quantum Singular Value Transformation to invert the matrix. The inverted matrix provides the Green’s function corresponding to the physical system for determining one or more physical properties of the physical system. Figure 19 is an example of a computing system (platform) for performing the processing (operations) shown in Figure 18. The computing system comprises a classical computing system 210 and a quantum computing system 250. The classical computing system 210 may comprise one or more computers (devices), likewise the quantum computing system 250 may comprise one or more computers (devices). The quantum computing system 250 is shown including at least one quantum circuit 260 which includes at least one qubit and at least one gate 255. The quantum circuit 260 can be considered as somewhat analogous to a compiled program (low-level code) which has been adapted to run on the specific hardware implementation of the quantum computer, such as the number and connectivity of the qubits and gates 255. In Figure 19, the “Quantum” of computing system 250 is shown in brackets because in some cases, the quantum computing system 250 may be replaced (implemented) by a classical computing system which simulates (emulates) the quantum computing system. In some cases, the quantum circuit 260 may be run on such a simulator to test and control the operation of the quantum circuit 260 to support software and / or hardware development. For example, the simulator may implement a noise-free environment, allowing for quantum circuits to be tested without having to specifically address issues of decoherence. The classical computing system 210 is shown in Figure 19 as including a compiler 220 and a control facility 225. The compiler 220 may be used to process higher level code to produce desired quantum circuits 260. The control facility 225 may manage and implement various interoperations between the classical computing system 210, for example, transferring a compiled quantum circuit (program) to the quantum computing system 250 for execution. The approach described herein allows the quantum singular transform to be used to calculate the Green’s function in the Lehmann representation, given a sufficient polynomial approximation of theinverse function. In the noise-free simulations presented herein, large errors are observed when performing matrix inversion via QSVT if the singular values of the matrix to invert fall outside the well- defined domain of the polynomial approximation of the inverse function. These results indicate that care must be taken to ensure this doesn’t happen. Equation 16 above offers a route for determining a proper approximation level (as discussed above, this may involve further approximations to evaluate). The approach described herein has been used in the context of the metal insulator phase transition of the two-site Anderson model and Green’s functions obtained from noise-free QSVT simulations (see Figure 9). For k = 50 we find the absolute error from the true answer remains at the 10−12level. For the k = 10 result, where some of the singular values of the problem lie outside the domain where the polynomial approximation is correctly defined, much larger errors are obtained. Described herein is a new circuit strategy to implement the SELECT oracle in the LCU technique (following the approach of Ralli et al. in

[0032] and da Silva et al. in

[0045] ). The overall circuit cost to perform a block encoding of any matrix supplied as a linear combination of Pauli operators scales as O(|A|( ⌈ log2 (|A|) ⌉2+ ns)) CNOT and single-qubit gates. Here ns is the number of system qubits and |A| is the number of Pauli operators in the matrix defined in the LCU. The circuit cost of the SELECT operator may potentially be improved by the use of ‘phase gadgets’

[0044] . Furthermore, the ZX- calculus may potentially be used to maximize gate cancellations between multi-control Pauli gates in USELECTimplemented according to the template outlined in Figure 2 above. The approach described herein therefore includes receiving a given H, where we want to block encode different shifts of H according to equation 2a and 2b - in effect, writing a quantum circuit that gets the task onto a quantum computer or other such device. The block encoding include a ‘SELECT’ part based on a combination of teachings from

[0032] and

[0045] . Once the block encoding circuit for the shifted H has been obtained, the inverse of H is determined. This may be done by running the quantum singular value (QSVT) algorithm. If we refer to the quantum circuit for this step as C, then an evaluation of the Green’s function (see equation 2) may involve the following: a. First apply a creation / annihilation operator in a quantum circuit on qubit index i or j. b. Then apply quantum circuit C. c. Then apply the annihilation / creation operator on qubit index j or i. d. This will give the i,j th entry of the Green’s function, which can be evaluated through a Hadamard test.In order to be successful, we must select a value of k greater than the minimum singular value (sigma_min) of the block encoded matrix: (i.e. sigma_min >= 1 / k). Since finding sigma_min is the same as solving our problem, an approximation may be used as per equation 16 (or as per various alternatives to the approach of equation 16).Section 7 – Appendix Appendix A: Spectral function plots Figures 10, 11 and 12 give the spectral plots for the two-site single particle Anderson model at different U, V values with k=10 results. In this context, U is the onsite Coulomb repulsion and V is the interaction of this bath site with the impurity. The value of k determines the size of the cut-out region for the polynomial approximation for the matrix inverse 1 / x at x ~ 0. The larger the value of k, the smaller the cut-out. Having a smaller cut-out provides for determining a more accurate inverse, but involves additional processing. In particular, Figure 10 (top) shows a plot of a spectral function (as per equation 17) of the single particle Anderson model for U = 2 and V = 0.943. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 10. Figure 10 (middle) shows a plot of the absolute error of spectral function ΔA(ω) = |AQSVT (ω) − Atrue(ω)|. Figure 10 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular value below the horizontal black line is outside the region where the approximation of 1 / x is well defined. Figure 11: (top) shows a plot of spectral function (equation 17) of the single-particle Anderson model for U = 5.99 and V = 0.058. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 10. Figure 11 (middle) shows a plot of the absolute error of the spectral function ΔA(ω) = |AQSVT(ω) − Atrue(ω)|. Figure 11 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular below the horizontal black line is outside the region where the approximation of 1 / x is well defined. Figure 12 (top) shows a plot of spectral function (equation 17) of the single-particle Anderson model for U = 8 and V = 0. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 10. Figure 12 (middle) shows a plot of the absolute error of spectral function ΔA(ω) = |AQSVT(ω) − Atrue(ω)|. Figure 12 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Anysingular value below the black horizontal line is outside the region where the approximation of 1 / x is well defined. Figures 13, 14 and 15 give the spectral plots for the two-site single particle Anderson model at different U, V values. The spectral plots are directly analogous to those of Figures 10-12, except they use a higher value of k, namely k=50. Figure 13 (top) shows a plot of spectral function (equation 17) of the single-particle Anderson model for U = 2 and V = 0.943. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 50. Figure 13 (middle) shows a plot of the absolute error of spectral function ΔA(ω) = |AQSVT(ω) − Atrue(ω)|. Figure 13 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular below the black horizontal line is outside the region where the approximation of 1 / x is well defined. Figure 14 (top) shows a plot of spectral function (equation 17) of the single-particle Anderson model for U = 5.99 and V = 0.058. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 50. Figure 14 (middle) shows a plot of the absolute error of spectral function ΔA(ω) = |AQSVT(ω) − Atrue(ω)|. Figure 14 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular value below the horizontal black line is outside the region where the approximation of 1 / x is well defined. Figure 15 (top) shows a plot of spectral function (equation 17) of the single-particle Anderson model for U = 8 and V = 0. The red line and orange points (diagonal crosses) show the spectral function being calculated via exact diagonalization. The blue points (non-diagonal crosses) show the spectral function calculated via the quantum singular value transform with k = 50. Figure 15 (middle) shows a plot of the absolute error of spectral function ΔA(ω) = |AQSVT (ω) − Atrue(ω)|. Figure 15 (bottom) shows a plot of the minimum singular value of the matrix to undergo inversion via the QSVT. Any singular value below the horizontal black line is outside the region where the approximation of 1 / x is well defined.Appendix B - Mott Phase transition Figures 16 and 17 are used to summarise the Metal to insulator Mott phase transition results for different polynomial approximations of the inverse function. The results are given in Figures 16 and 16. For k = 10 it (Figure 16) it should be noted that the error in the spectral function may change by many orders of magnitude compared, whereas for k = 50, the error in the results remains much smaller, around the pico level of accuracy at all points. Figure 16 depicts a metal to insulator Mott transition of the two site Anderson impurity model. The left plots show the density of states on the Bethe lattice and the right plots are absolute error in the density of states with respect to the exact solution (obtained via classical matrix inversion). The first plot gives the non-interacting system, followed by (U, V) combinations of (2, 0.943), (4, 0.745), (5.99, 0.058), (8, 0). Note the free density of states on the Bethe lattice with infinite coordination is ρ0(x) = ^^^^^√4^^− ^^(where we take t=1). The interacting density of states is ρ(ω) = ρ0(ω + μ − Σ[ω]). The ^ non-interacting Green’s function is defined as|^|-1(ω)= (ω + iδ − ^^+ μ − ^ ^ ^^ ) (where we take ^^=0). For the final plot (8, 0), particle hole symmetry was broken and the spin up and down parts of Σ[ω]) were treated separately and combined. The yellow data points (diagonal crosses) show the results obtained from a noise free QSVT simulation with k = 10. Figure 17 depicts a metal to insulator Mott transition of the two site Anderson impurity model. The left plots show the density of states on the Bethe lattice and the right plots are absolute error in the density of states with respect to the exact solution (obtained via classical matrix inversion). The first plot gives the non-interacting system followed by (U, V) combinations of (2, 0.943), (4, 0.745), (5.99, 0.058), (8, 0). Note the free density of states on the Bethe lattice with infinite coordination is ρ0(x) = ^^^^^√4^^− ^^(where we take t=1). The interacting density of states is ρ(ω) = ρ0(ω + μ − Σ[ω]). The non-interacting Green’s function is defined as ^^^|^|^^^(ω)= (ω + iδ − ^^+ μ −-1^ ^ ^^ ) (where we take ^^=0). For the final plot (8, 0), particle hole symmetry was broken and the spin up and down parts of Σ[ω]) were treated separately and combined. The yellow data points (diagonal crosses) show the results obtained from a noise free QSVT simulation, where k = 50.References [1] G. Green, An essay on the application of mathematical analysis to the theories of electricity and magnetism, Vol.3 (author, 1889). [2] D. G. Duffy, Green’s functions with applications (Chapman and Hall / CRC, 2015). [3] L. Challis and F. Sheard, Physics Today 56, 41 (2003). [4] “The nobel prize in physics 1965,” https: / / www.nobelprize.org / prizes / physics / 1965 / summary / , accessed: 2022-08-15. [5] J. Schwinger, arXiv preprint hep-ph / 9310283 (1993). [6] M. Kline, Mathematical Thought from Ancient to Modern Times: Volume 2, Vol.2 (Oxford university press, 1990) p.684. [7] A. Szabo and N. S. Ostlund, Modern quantum chemistry: introduction to advanced electronic structure theory (Courier Corporation, 2012). [8] G. Onida, L. Reining, and A. Rubio, Reviews of modern physics 74, 601 (2002). [9] N. H. March, Electron Correlation in the Solid State (Imperial College Press and Distributed by World Scientific Publishing Co, 1999).

[0010] I. Rungger, N. Fitzpatrick, H. Chen, C. Alderete, H. Apel, A. Cowtan, A. Patterson, D. M. Ramo, Y. Zhu, N. H. Nguyen, et al., arXiv preprint arXiv:1910.04735 (2019).

[0011] X. Cai, W.-H. Fang, H. Fan, and Z. Li, Physical Review Research 2, 033324 (2020).

[0012] S. Endo, I. Kurata, and Y. O. Nakagawa, Physical Review Research 2, 033281 (2020).

[0013] H. Chen, M. Nusspickel, J. Tilly, G. H. Booth, et al., Physical Review A 104, 032405 (2021).

[0014] F. Jamet, A. Agarwal, C. Lupo, D. E. Browne, C. Weber, and I. Rungger, arXiv preprint arXiv:2105.13298 (2021).

[0015] R. Sakurai, W. Mizukami, and H. Shinaoka, Physical Review Research 4, 023219 (2022).

[0016] J. Zhu, Y. O. Nakagawa, Y.-S. Zhang, C.-F. Li, and G.-C. Guo, New Journal of Physics 24, 043030 (2022).

[0017] T. Steckmann, T. Keen, A. F. Kemper, E. F. Dumitrescu, and Y. Wang, arXiv preprint arXiv:2112.05688 (2021).

[0018] B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer, Physical Review X 6, 031045 (2016).

[0019] T. Kosugi and Y.-i. Matsushita, Physical Review A 101, 012330 (2020).

[0020] Y. Tong, D. An, N. Wiebe, and L. Lin, Physical Review A 104, 032422 (2021).

[0021] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (2019) pp.193–204.

[0022] P. W. Anderson, Physical Review 124, 41 (1961).

[0023] M. Potthoff, Physical Review B 64, 165114 (2001).

[0024] H. Lehmann, Il Nuovo Cimento (1943-1954) 11, 342 (1954).

[0025] M. Hjorth-Jensen, M. P. Lombardo, and U. Van Kolck, Springer Lecture Notes in Physics 936 (2017).

[0026] J. Kempe, A. Kitaev, and O. Regev, Siam journal on computing 35, 1070 (2006).

[0027] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).

[0028] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, PRX Quantum 2, 040203 (2021).

[0029] D. Camps, L. Lin, R. Van Beeumen, and C. Yang, arXiv preprint arXiv:2203.10236 (2022).

[0030] A. M. Childs and N. Wiebe, arXiv preprint arXiv:1202.5822 (2012).

[0031] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical Review Letters 114, 090502 (2015).

[0032] A. Ralli, P. J. Love, A. Tranter, and P. V. Coveney, Physical Review Research 3, 033195 (2021).

[0033] G.-L. Long and Y. Sun, Physical Review A 64, 014303 (2001), arXiv:0104030 [quant-ph].

[0034] M. Mottonen, J. J. Vartiainen, V. Bergholm, and M. M. Salomaa, Quantum Information and Computation 5, 467 (2004), arXiv:0407010 [quant-ph].

[0035] V. Shende, S. Bullock, and I. Markov, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25, 1000 (2006).

[0036] I. F. Araujo, D. K. Park, F. Petruccione, and A. J. da Silva, Scientific Reports 11, 6329 (2021), arXiv:2008.01511.

[0037] G. H. Low and I. L. Chuang, Physical Review Letters 118, 010501 (2017).

[0038] L. K. Grover, Physical Review Letters 80, 4329 (1998).13

[0039] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Contemporary Mathematics 305, 53 (2002).

[0040] T. J. Yoder, G. H. Low, and I. L. Chuang, Physical Review Letters 113, 210501 (2014).

[0041] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in Proceedings of the forty-sixth annual ACM symposium on Theory of computing (2014) pp.283–292.

[0042] B. Yan, S. Wei, H. Jiang, H. Wang, Q. Duan, Z. Ma, and G.-L. Long, Scientific Reports 12, 1 (2022).

[0043] M. B. Hastings, D. Wecker, B. Bauer, and M. Troyer, arXiv preprint arXiv:1403.1539 (2014).

[0044] A. Cowtan, S. Dilkes, R. Duncan, W. Simmons, and S. Sivarajah, arXiv preprint arXiv:1906.01734 (2019).

[0045] A. J. da Silva and D. K. Park, arXiv preprint arXiv:2203.11882 (2022).

[0046] S. S. Bullock and I. L. Markov, arXiv preprint quantph / 0303039 (2003).

[0047] V. V. Shende, S. S. Bullock, and I. L. Markov, in Proceedings of the 2005 Asia and South Pacific Design Automation Conference (2005) pp.272–275.

[0048] V. von Burg, G. H. Low, T. Häner, D. S. Steiger, M. Reiher, M. Roetteler, and M. Troyer, Physical Review Research 3, 033055 (2021).

[0049] Y. Dong, L. Lin, and Y. Tong, PRX Quantum 3, 040305 (2022).

[0050] S. Chakraborty, A. Gilyén, and S. Jeffery, arXiv preprint arXiv:1804.01973 (2018).

[0051] A. W. Harrow, A. Hassidim, and S. Lloyd, Physical Review Letters 103, 150502 (2009).

[0052] L. Lin, arXiv preprint arXiv:2201.08309 (2022).

[0053] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, PRX Quantum 2, 010342 (2021).

[0054] G. H. Low, T. J. Yoder, and I. L. Chuang, Physical Review X 6, 041067 (2016).

[0055] J. Haah, Quantum 3, 190 (2019).

[0056] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, arXiv preprint arXiv:2003.02831 (2020).

[0057] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Physical Review A 103, 042419 (2021).

[0058] J. Martyn, A. Tan, C. Huang, and I. chuang, “pyqsp,” https: / / github.com / ichuang / pyqsp (2022).

[0059] Y. Dong, X. Meng, J. Wang, and l. Lin, “QSPPACK,” https: / / github.com / qsppack / QSPPACK (2022).

[0060] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, arXiv preprint arXiv:1806.01838 (2018).

[0061] C. Jordan, Bulletin de la Société mathématique de France 3, 103 (1875).

[0062] G. Kotliar and D. Vollhardt, Physics today 57, 53 (2004).

[0063] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Reviews of Modern Physics 68, 13 (1996).

[0064] M. Caffarel and W. Krauth, Phys. Rev. Lett.72, 1545 (1994).

[0065] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys.68, 13 (1996).

[0066] J. M. Kreula, L. García-Álvarez, L. Lamata, S. R. Clark, E. Solano, and D. Jaksch, EPJ Quantum Technology 3, 1 (2016).

[0067] P. Jordan and E. Wigner, Zeitschrift fur Physik 47, 631 (1928).

[0068] A. Tranter, C. D. Paola, D. M. Ramo, D. G. David Zsolt Manrique, G. Greene-Diniz, G. Christopoulou, I. Polyak, J. P. Irfan Khan, J. Kirsopp, K. Yamamoto, M. Tudorovskaya, M. Krompiec, and N. Fitzpatrick, “Introduction to the inquanto computational chemistry platform for quantum computers,” https: / / medium.com / cambridge-quantum computing / 4fced08d66cc (2022).

[0069] Quantinuum, “Inquanto,” https: / / www.quantinuum .com / products / inquanto (2022).

[0070] S. Sivarajah, S. Dilkes, A. Cowtan, W. Simmons, A. Edgington, and R. Duncan, Quantum Science and Technology 6, 014003 (2020).

[0071] C. R. Harris, K. J. Millman, S. J. van derWalt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).

[0072] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).

[0073] Y. Hong and C.-T. Pan, Linear Algebra and its Applications 172, 27 (1992).

[0074] G. Piazza and T. Politi, Journal of Computational and Applied Mathematics 143, 141 (2002).

[0075] T.-Z. Huang, Computers & Mathematics with Applications 55, 1075 (2008).

[0076] L. Zou and Y. Jiang, Linear algebra and its applications 433, 1203 (2010).

[0077] L. Zou, J. Math. Inequal 6, 625 (2012).

[0078] Z. Bai, G. Fahey, and G. Golub, Journal of Computational and Applied Mathematics 74, 71 (1996).

[0079] I. C. Ipsen and D. J. Lee, arXiv preprint arXiv:1105.0437 (2011).

[0080] F. Aryasetiawan and O. Gunnarsson, Reports on Progress in Physics 61, 237 (1998).

[0081] M. Eckstein, M. Kollar, K. Byczuk, and D. Vollhardt, Physical Review B 71, 235119 (2005).

[0082] X. Zhang, M. Rozenberg, and G. Kotliar, Physical Review Letters 70, 1666 (1993).

[0083] P. Sriluckshmy, M. Nusspickel, E. Fertitta, and G. H. Booth, Physical Review B 103, 085131 (2021).

[0084] G. H. Low, V. Kliuchnikov, and L. Schaeffer, arXiv preprint arXiv:1812.00954 (2018).

[0085] K. Wan, Quantum 5, 380 (2021).

[0086] D. Camps and R. Van Beeumen, arXiv preprint arXiv:2205.00081 (2022).

Claims

Claims 1. A computer-implemented method comprises: providing a Hamiltonian representation of a physical system; using a first quantum circuit to block encode the Hamiltonian into a matrix M; and using a second quantum circuit to perform a Quantum Singular Value Transformation (QSVT) to invert the matrix, wherein the inverted matrix is used to provide a Green’s function corresponding to the physical system for determining one or more physical properties of the physical system.

2. A computer-implemented method according to claim 1, wherein the first quantum circuit implements a double-factorisation of an operator of the Hamiltonian and a mid-circuit measurement to determine the values of the matrix M.

3. The computer-implemented method of claim 1 or claim 2, wherein the Green’s function is expressed in matrix form (i, j) and comprises an advanced component according to equation 2a and a retarded component according to equation 2b.:

4. The computer-implemented method of claim 3, wherein z = ω + iδ is a complex frequency in which w >> δ.

5. The computer-implemented method of claim 2, 3 or 4, wherein the block encoding is performed using a linear combination of unitaries (LCU), optionally wherein the LCU is given as a linear combination of Pauli operators.

6. The computer-implemented method of any of claims 2 to 5, wherein the block encoding is performed with respect to the following complex-shifted Hamiltonians:

7. The computer-implemented method of claim 6, where the QSVT is configured to approximately invert the singular values of the block encoded matrix using the block encodings of^(^)and ^(^)to allow separate calculation of the inverse parts of equations 2a and 2b to reduce the cir depth.

8. The computer-implemented method of any preceding claim, wherein inverting the matrix includes generating quantum signal-processing angles to implement a desired inverse function that will be applied to the singular values of the block encoded matrix.

9. The computer-implemented method of claim 8, wherein the desired inverse function is a polynomial approximation of f(x) = 1 / x.

10. The computer-implemented method of claim 9, wherein the polynomial expression is defined in the range x = [−1, −1 / κ] ∪ [1 / κ, 1].

11. The computer-implemented method of any of claims 8 to 10, further comprising selecting k to be large enough such that there is no singular value in the range (-1 / k) – (+1 / k).

12. The computer implemented method of any of claims 8 to 11, further comprising constructing the second quantum circuit to implement the QSVT algorithm using the matrix M and the quantum signal-processing angles.

13. The computer-implemented method of any preceding claim, wherein the Hamiltonian is not unitary but the matrix M is unitary.

14. The computer-implemented method of any preceding claim, wherein the inverted matrix provides matrix elements for calculating the single-particle, many-body Green’s function corresponding to the physical system.

15. The method of any preceding claim, further comprising performing a Hadamard test to evaluate the real and complex parts of each entry in the Green’s function.

16. The method of any preceding claim, further comprising using the Green’s function to calculate one or more physical properties of the physical system, optionally wherein the physical properties include at least one of excitation and ionization energies, ground-state energies,transition matrix elements, absorption coefficients, dynamical polarizabilities, and elastic and inelastic electron cross-sections.

17. The method of any preceding claim, wherein the Green’s function is determined for a two- site, single-impurity Anderson model, optionally for approximating the Mott insulator phase transition.

18. The method of any preceding claim, wherein the method is implemented at least partly on a quantum computing system.

19. A computer system configured to: provide a Hamiltonian representation of a physical system; use a first quantum circuit to block encode the Hamiltonian into a matrix M; use a second quantum circuit to perform a Quantum Singular Value Transformation (QSVT) to invert the matrix, wherein the inverted matrix is configured to provide matrix elements for calculating a Green’s function corresponding to the physical system, for determining one or more physical properties of the physical system.

20. A computer system according to claim 19, wherein the first quantum circuit implements a double-factorisation of an operator of the Hamiltonian and a mid-circuit measurement to determine the values of the matrix M.

Images