Digital quantum simulation of gravitational optomechanics with IBM quantum computers

We showcase the digital quantum simulation of the action of a Hamiltonian that governs the interaction between a quantum mechanical oscillator and an optical field, generating quantum entanglement between them via gravitational effects. This is achieved by making use of a boson-qubit mapping protocol and a digital gate decomposition that allow us to run the simulations in the quantum computers available in the IBM Quantum platform. We present the obtained results for the fidelity of the experiment in two different quantum computers, after applying error mitigation and post-selection techniques. The achieved results correspond to fidelities over 90%, which indicates that entanglement was indeed generated through the simulation of a quantum gravitational field.


Introduction
In the last years a novel approach to the open debate of quantum gravity has emerged, where the focus has shifted to just prove the quantum nature of gravity, without disclosing the underlying full quantum theory [1][2][3]5].The key aspect is that, since two quantum systems cannot be entangled via a classical channel if they were not entangled beforehand, if we are able to show that the gravitational field is able to generate entanglement, we can conclude that it must arise from genuine quantum properties, and therefore the field must be a quantum entity.
While the prospects for an experimental implementation of these proposals are encouraging -as compared to a direct proof of the full quantum theory of gravityactual experimental results are still missing.In this context, a quantum simulation in a quantum computer [6,7] could be helpful by informing the ongoing experimental efforts, leveraging the experimental amenability and tunability of modern quantum platforms [8][9][10].Some of us have recently introduced a recipe for the digital quantum simulation of bosonic Hamiltonians [11], which combines Trotter [7] and gate-decomposition techniques with a boson-qubit mapping [12,13] to encode the Hamiltonian into a sequence of single-qubit and two-qubit gates.We have applied this scheme to the high-fidelity digital quantum simulation of paradigmatic quantum-optics interactions, such as beam-splitters or single-mode and two-mode squeezers [14].
The procedure carried out in this project is similar to the recent work by one of the authors in [15], with the main difference laying in the fact that, in that work, arXiv:2401.08370v3[quant-ph] 5 May 2024 the author reports the digital quantum simulation of the Hamiltonian involved in the generation of quantum entanglement between a pair of quantum harmonic oscillators by gravitational means [3], while here we apply this scheme in order to study the quantum gravitational Hamiltonian that represents the interaction of a matter system and an optical field.The study of this Hamiltonian is related to the field of gravitational optomechanics, which is a term that was introduced in [5] as a way to refer to any optomechanical interaction that has the gravitational field as its source.
This interaction gives rise to a bosonic two-mode Hamiltonian which we translate into a multiqubit one and then decompose into one and two-qubit quantum gates in order to launch it to IBM quantum devices.Even though the interaction studied in this work is exactly solvable, it provides a framework towards simulating more complex scenarios that could not be studied with analytic treatments.The techniques presented here are the basic building blocks of interactions including more modes and more excitations per mode, which would give rise to those scenarios.
We consider a simplified scenario allowing a maximum of one excitation per mode -a sensible approximation given the weakness of the interaction-and leverage error mitigation and postselection techniques to achieve fidelities above 90%.Due to some technical limitations, which are explained in detail later, we were not able to reproduce the current experimental range of parameters, which purely depends on the particular experimental setup that we have studied, so we consider slightly larger values, which amplify the amount of quantum gravitational entanglement generated.Thus, these preliminary four-qubit experiments can be considered as a first step towards larger and more elaborate quantum simulations, with more modes and more photons per mode, beyond the capabilities of classical computers.
The structure of the paper is the following.In the next section, we introduce the Hamiltonian of gravitational optomechanics and develop its digitization.After that, we show how to realize a full quantum state tomography and compute the fidelity of this state with respect to theoretical predictions.Then, we present experimental results for the fidelity of the digital quantum simulation in IBM quantum computers.We conclude with a summary and discussion of our results.

Digitization of gravitational optomechanics
We are interested in performing the digital quantum simulation of the interaction presented in [5] as the one that rules the interaction betweeen a matter-wave and a photonic field mode.This interaction is governed by the following Hamiltonian: where g 0 is an optomechanical coupling whose value is related to the parameters of the experiment and b † ( b) and â † (â) are the creation (annihilation) operators of the mechanical oscillator and photon, respectively.What the authors in [5] realised is that (1) is formally the same as the Hamiltonian interaction that arises in the field of optomechanics when studying a quantum system composed of a cavity field interacting with a movable mirror [16,17], which means we can apply existing protocols from quantum optomechanics to our interaction.One of these procedures is usually performed when the cavity is strongly driven to some large coherent amplitude α = ⟨â⟩ [18], and it consists in linearising the interaction by writing â → α + δâ, where δâ represents the small quantum fluctuations around this mean value α.This approximation is only valid when two conditions are met [4], |α| ≫ 1 and µσ ≪ 1, where µ ∼ g 0 t and σ ∼ Var( b + b † ).Both of them are fulfilled in our case, so, by renaming δâ as â, we get: with g = g 0 α.This is formally equivalent to the well-known cavity optomechanical interaction of quantum optics in the linearized regime [16][17][18], but now the coupling is induced by the interaction with the gravitational field encoded in the coupling g.
Thus, in order to perform and analyze the simulation of our Hamiltonian, our goal is to develop a quantum circuit that replicates its action when applied on a certain initial state with no entanglement -such as the ground state-, prove that it generates entanglement between the matter and photonic systems and give an estimate of the accuracy of the simulation by calculating the fidelity of the state with respect to an ideal theoretical prediction.
We will now work upon the Hamiltonian presented in (2).The time evolution of the initial state will be (ℏ = c = 1): The main goal of this work will be the digitalization of the following unitary operator: where we have defined: For this purpose, we will need to use the boson mapping techniques presented in [12,13].It is shown there that it is possible to map N bosonic modes with N P excitations each to N (N P + 1) qubits.In our case, we are dealing with two modes, corresponding to the matter and photon states in (2).For the rest of this work, we will consider N P = 1 for both bosonic modes.We then will need four qubits in all our quantum circuits, which we will label from 0 to 3. Following [12,13], we will have that: where the subscripts m and p refer to the matter and photonic states, respectively.The mapping of the operators for Now, we need to express the ( b + b † ) and the (â + â † ) terms in ( 2) with respect to the Pauli matrices σ x , σ y , σ z .For the first one, we have: and therefore: Similarly: and finally we get the expression for the full Hamiltonian: It can be shown that all of the terms in (11) commute with each other, so we can factorize the unitary without the need of Trotter approximations.For instance, let us consider the first and last terms.Clearly, the operators acting on different qubits commute, so we have: In a similar way: Recalling the commutation relationships of Pauli matrices, we get: and thus the two terms in (12) vanish and so does the full commutator.An equivalent procedure can be carried out for the rest of the terms.Therefore we can write: where the U (i) ε are obtained as the exponential of each of the terms in (11).For instance: with ε = ε/4.In principle, these are already unitary gates that we could implement in a quantum circuit to simulate our Hamiltonian.However, these unitary gates are not available in the IBM quantum devices, which means that we have to carry out a gate decomposition to express each of them in terms of gates that are actually available in the IBM software [11,14,15].Let U be some unitary operation such that, for a Hamiltonian H, it satisfies: where H 0 is some convenient single-qubit operation.Then, we can always write: Our goal is now to carry out this procedure for each piece of our Hamiltonian.For example, for U ε , by defining H 0 = −σ 0 z , we have that, using ( 17): Thus, we find that U (1) ε can be rewritten as: with The decompositions for the rest of the U (i) ε can be found by adding the necessary e iπ/4σ (i) z to the operator in (21), in order to turn the corresponding σ x into σ y in each case.As stated earlier, the purpose of this gate decomposition technique was to be able to express our Hamiltonian operator in terms of gates that can be implemented in the IBM quantum computers.Taking a look at the decomposition of each of our terms (see Table 1), it can be noticed that there are only three different types of operators that were needed throughout the whole process: The ones with σ i x , σ j z or σ i x σ j z in the exponent.All three of these operators can be converted into gates available in the IBM architecture [14,15].For the first two, we can use the single-qubit rotations R x (θ) and U 1 (λ): where the U 1 gate takes the particular form of the S gate for λ = π/2: For the last type of gates, the conversion is given by: where the CNOT gate is available in IBM for some pairs of qubits, depending on the particular device connectivity.
We can now put all this together to express each of the terms of the unitary with a combination of these three types of gates.For U (1) ε , we will have that, up to an irrelevant global phase and after some cancellations and simplifications: For the sake of completeness, we leave the rest of the digital decompositions for the It is also worth mentioning that when we send a quantum circuit to an IBM quantum computer, the circuit where the measurements are actually performed is usually not exactly the same as the one we designed, since the circuit is first transpiled, that is, optimized to the particular quantum architecture used.
The gate set of the computers that our computations are being carried out in can be checked accessing their configuration.All of the accessible quantum computers possess the same gate set, formed by the X, CNOT, SX (also referred to as √ X, and equivalent to a R x (π/2) gate) and R z gates, on top of the identity and the reset gates.Therefore, every single gate we place in our circuit will be transformed to a combination of some of these gates.
The transpilation of a circuit might increase or decrease the number of gates that are needed to fulfill the action of the circuit with respect to the original one.We are in general interested in having the circuit simplified after the transpilation, in order to reduce the number of gates and therefore the global error.However, connectivity issues might translate into additional gates.
In  Note that the initial two X-gates in Figure 1 are introduced in order to prepare the system into the initial state which we worked with, which is the ground state: By observing these three images, we can see that the transpilation process helps slightly simplify the circuit by decreasing the total number of single-qubit gates required to simulate the action of U ε without increasing the number of two-qubit gates.For both quantum computers, the circuit consists of 32 single-qubit gates and

Fidelity
In order to determine the quality of our digitalization, we calculate the fidelity as [14,15,19]: where |ψ⟩ is the state at which the system will be after the action of the full quantum circuit (that is, after the action of U ε ), while ρ will be the density matrix that represents the state actually obtained in the experiment.We will now analyze the procedure that is required to obtain an expression for the fidelity as a function of ε.
Considering again U ε , if we restrict ourselves to perturbative values of ε, we can perform a Taylor expansion of U ε up to second order in ε and provide an estimate to the state of the system after a certain time, by applying this expression for U ε to our initial state |ψ 0 ⟩.We arrive to: which is, as we can see, an entangled state involving the ground and first excited states of the matter system and the optical field.In particular its concurrence, defined for pure states as [20,21]: where σ m is the reduced density matrix of the matter oscillator subsystem, can be quiqkly found to be, up to second order: Let us now consider the experimental density operator.We perform a full state tomography of ρ.For a two qubit state, we have [19]: where the sum is over vectors ⃗ v = (v 1 , ..., v n ) with entries v i chosen from the set 0, 1, 2, 3. Plugging this into (28), we have: where i, j = 0, 1, 2, 3 and σ 0 is defined as the identity matrix.Using eq. ( 29) we obtain that the nonzero terms in eq. ( 33) are: which leads us to: The traces in eq. ( 35) are expected values of observables of the form σ i,m ⊗ σ j,p over ρ, and we can express them as [14]: where these quantities are the probabilities of each mode being in state |0⟩ or |1⟩ of the corresponding basis in each case.For instance, p XY (|0⟩ m |0⟩ p ) represents the probability of the matter and photon systems of being in the |0⟩ state of bases X and Y respectively.Nonetheless, the measuring operator available in the IBM system is equivalent to measuring in the Z basis and there is no way to perform a measure directly in the other two.This can be solved by implementing quantum gates in the relevant qubits before the measurement, flipping them to X or Y basis.In particular, we apply a Hadamard gate before measurement in the case of X-measurements, and S † and Hadamard for Y -measurements.
With this, we are already able to perform all the necessary measurents for the computation of the fidelity of our experiment.Moreover, in order to improve the fidelity, we apply error mitigation in the readout measurements and postselect the results to the subspace of the 4-qubit Hilbert space with physical meaning in our simulation, namely the one spanned by the four states in (6).

Results
We have now developed all the necessary tools to implement the quantum simulation of the Hamiltonian in ( 2), but we still must consider an appropriate range of values of ε.
From (5), we recall that the value of ε is determined by that of the gravitational optomechanical coupling g in the Hamiltonian and the experimental time for which we are assuming we would let the system evolve.In [5], it is mentioned that, with the currently available experimental setups, a value of 2π • 10 −11 for g could be reached and adequate values of t would be of the order of seconds, giving rise to a value of ε around 10 −10 .However, IBM's Qiskit approximates to 0 any number smaller than 10 −9 .Notice, however, that the amount of entanglement generated by the quantum gravitational field is proportional to ε, according to eq. ( 31).Thus, it would be of interest to simulate the system for values of ε that magnify the generated entanglement while lying beyond current experimental reach -but remaining in the perturbative regime.In Figure 4, we show the results for the fidelity of the state after the action of the circuit for values of ε that range from 10 −7 to 10 −2 .The corresponding measurements were carried out in two different IBM quantum computers, namely Belem and Nairobi, with 5 and 7 qubits and 16 and 32 quantum volume [22] respectively.The average error rates of the IBM Belem (Nairobi) quantum computer are of around 2.76 (3.28) × 10 −4 for the single-qubit gates,  As we can see, we have achieved values for the fidelity between 0.9 and 0.95 in both cases, after applying error mitigation and postselection.The results are improved to a small extent when the measurements are done in the IBM Nairobi computer, raising the average fidelity from around 92% to 94%, due to its larger quantum volume.

Conclusions
In this work, we carry out the digital quantum simulation of a Hamiltonian governing the interaction between a mechanical oscillator and an optical field which interact through a quantum gravitational field generating entanglement between them.Therefore, we simulate an experiment for the generation of gravitational entanglement to prove the quantum nature of gravity.This is achieved by implementing the corresponding boson-qubit mapping and performing a suitable gate decomposition in order to write our time evolution operator U ε in terms of gates that can be used in the IBM Quantum platform computers.Exploiting post-selection and error mitigation techniques, we achieve high-fidelity -above 90%-in two different devices, and observe an improvement of the fidelity in the device with larger quantum volume.Given the rate at which IBM is increasing the quantum volume, it could soon be possible to increase the number of qubits while keeping high fidelity.This would allow us to relax the restriction to one maximum excitation per mode.While this should be a very good approximation given the weakness of the coupling, increasing the allowed excitations and thus the number of qubits will enable to reach the post-classical regime, going beyond the capabilities of classical computers.As mentioned in the introduction, this is the reason why we found that it would be convenient to carry out this simulation protocol in a quantum computer: While the interaction that was presented here can be studied analytically, it would be computationally harder to do so with more complex versions of the process, either with more modes or excitations per mode.For instance, with the currently available 127-qubit computers, we could consider more than 50 excitations per mode, which would also allows us to consider larger coupling values and therefore simulate the gravitational entanglement generated by large gravitational fields.While the quantum volume of these devices have already raised to 128, it is likely that it would not be enough to compensate for the dramatic improvement in the number of quantum gates that would be required to simulate such a large number of excitations.However, novel and more sophisticate error mitigation techniques are currently available [23], enabling error mitigation in the full circuit and not only in the readout, which would presumably improve the fidelity of those large post-classical experiments.Meanwhile, our preliminary four-qubit results represent a first step along that direction, since larger simulations would essentially require the same techniques described here.

Appendices
Appendix A: Gate decomposition The decompositions of the U (i) ε in terms of gates that are available in the IBM software are: (37)

Appendix B: Fidelity error
In order to estimate the error of the fidelity for each value of ε, we recall the expression in (35), in which we express the fidelity in terms of probabilities that are determined experimentally by performing measurements on the circuit in different bases.These measurements have an associated error, which we can use to estimate the error in the fidelity through standard error propagation: As a good approximation, we can consider that the readout error is the same for all qubits -the averager readout error provided by IBM-and doesn't depend on whether the measurement is meant to result in 0 or 1.Let us refer to this magnitude as λ.Since these erorr rates are small, the rate of error in the process of measuring n qubits can be taken, in the perturbative regime, as δ = n • λ.Therefore, every time that we perform N measures on the circuit to calculate each of these probabilities, the number of "counts" that carry along a mistake in the readout will be ∆ = δ • N .
Let us now consider the tr ZZ ρ term as an example.It was calculated as: where the N represent the counts associated to each of the probabilities that are calculated to estimate the value of the trace.Through linear error propagation, we can see that the error associated to both the numerator and the denominator will be ∆ tot ≡ ∆a = ∆b = ∆ + ∆ ′ + ∆ ′′ + ∆ ′′′ .The total error associated to this term will be: (42) Tables Table 1: Gate decomposition of each term in our Hamiltonian evolution operator.

Figures 1 ,
2 and 3, we show the full quantum circuit that is implemented to simulate the action of our time evolution operator and the transpiled version of said circuit in the quantum computers IBM Belem and IBM Nairobi.The transpiled circuits look different due to the different number of qubits and connectivity.It is also worth remarking that the value of ε does not affect the structure of the original quantum circuit or the transpiled versions, but only results in changes in the parameters in some of the quantum gates.

Figure 1 :
Figure 1: Quantum circuit for the simulation of U ε .
8.75 (14.92) × 10 −3 for the CNOT gates and 2.11 (3.06) × 10 −2 for the readout -however, recall that we have made use of an error mitigation technique for our measurements.It is important to remark that these parameters change on a daily basis.The experiments in Belem were carried out on June 10, 2023, while the ones in Nairobi were done on June 18, 2023.The error bars in the graphs are related to the readout error and the way they are computed is detailed in Appendix B.

Figure 4 :
Figure 4: Fidelity results of the digital quantum simulation procedure in the IBM Belem and IBM Nairobi computers.