Digital Quantum Simulation of Linear and Nonlinear Optical Elements

: We provide a recipe for the digitalization of linear and nonlinear quantum optics in networks of superconducting qubits. By combining digital techniques with boson-qubit mappings, we address relevant problems that are typically considered in analog simulators, such as the dynamical Casimir effect or molecular force ﬁelds, including nonlinearities. In this way, the beneﬁts of digitalization are extended in principle to a new realm of physical problems. We present preliminary examples launched in IBM Q 5 Tenerife.

Afterdecades of both theoretical and experimental efforts, a new generation of technologies is on the brink of delivering the heralded quantum revolution. For instance, large networks of superconducting qubits are close to proving quantum supremacy [1][2][3], opening a new era in quantum computing and quantum simulation.
While these promising applications are gaining a great deal of attention, qubits are not the only key players of the quantum world. Indeed, the current scenario has only been possible due to the impressive developments in quantum optical and quantum information setups in the last few decades, which managed to reach the single-atom and single-photon level. Therefore, electromagnetic fields or, more generally, bosonic field modes are also central to modern quantum setups. An example of the richness offered by the physics of the electromagnetic field is boson sampling [4], namely a computation of the number statistics of the output photons of a linear optics network, which can in principle be implemented in a quantum optical setup, but is widely believed to be intractable by classical means. Under certain conditions for the number of photons and modes, a boson sampling experiment would also prove quantum supremacy. However, while some small-scale experiments have been realized [5][6][7][8] and a number of promising proposals are available [9,10], a post-classical boson sampler has not yet been implemented in the laboratory. The ingredients needed in a boson sampling architecture are simple linear optics elements, such as beam-splitters and phase-shifters. Going a step further, we can also consider other paradigmatic ingredients such as two-mode squeezers or Kerr nonlinearities. Indeed, the combination of controllable two-mode squeezers and beam-splitter interactions between two bosonic modes is the basis of quantum information processing and quantum computation with continuous variables [11,12], such as the realization of Fredkin and exponential SWAPgates, which have recently been implemented experimentally with fidelities ranging from 60% to 90% [13]. Kerr nonlinearities can also find interesting applications, such as quantum metrology [14], whose simulation in a quantum computer might thus be of interest.
In this work, we aim to bring linear and nonlinear optics to the realm of cutting-edge qubit-based quantum computers and quantum simulators. To this end, we consider digital quantum simulation techniques [15], which combined with a boson-qubit mapping [16][17][18], allow us to encode generic multimode bosonic interactions into a sequence of single-qubit and two-qubit gates. Applications include the digital quantum simulation of continuous-variable quantum information processing and

Digitalization of Bosonic Hamiltonians
The aim of this work is to provide a recipe for the quantum simulation of bosonic Hamiltonians of interest with networks of qubits. The first step is to encode the bosonic Hamiltonian into the qubit network by means of a suitable boson-qubit mapping.

Mapping Bosons to Qubits
As shown in [16,17], it is possible to map N bosonic modes containing a maximum number of N p excitations each to N(N p + 1) qubits. For each bosonic mode, we are able to associate an (N P + 1)-qubit quantum state to each Fock state. If we consider the j th mode, we have: . . . . . . |N P j ↔ |1 0 1 1 1 2 · · · 0 N P j where |n j denotes a quantum state with n bosons in the j th mode. Notice that the state |n j is simulated by a state where out of the N P + 1 qubits that are associated with the j th mode, the only one that is in the state |0 is the n th qubit. The bosonic creation operator maps to: where a pair (n, j) refers to the n th qubit in the chain of qubits representing the j th bosonic mode. The Pauli creation and annihilation operators are given by: in terms of the Pauli matrices σ x and σ y . From Equation (3), it is straightforward to obtain as well the annihilation operator and then all the combinations that appear in bosonic Hamiltonians of interest; we will see relevant examples in the next sections.
Once we have encoded the bosonic Hamiltonian into a suitable qubit-network Hamiltonian and before we enter into the details of examples and applications, let us recall the basic notion of digital quantum simulations, namely the Suzuki-Trotter approximation.

Trotter-SuzukiDecomposition
The idea of the Trotter-Suzuki decomposition (see for instance [15]) is to decompose an involved non-local Hamiltonian dynamics into a sequence of experimentally amenable local gates. To this end:

•
The Hamiltonian is decomposed into a number m of terms: • Time is discretized, namely divided into s steps of duration t s : • The exponential of a sum of operators is approximated by the exponential of a product of operators.
The approximation is exact only in the case that all the Hamiltonian terms commute. Otherwise, it neglects all the commutators in the Baker-Campbell-Hausdorff formula.
Putting it all together, the dynamics is approximated by: that is s repetitions of the sequence of Hamiltonian terms given by H j . The error in this approximation can be bounded and controlled [15], being suppressed by a sufficiently large number of repetitions. However, the number of repetitions also increases the number of gates, which gives rise to larger experimental errors. Therefore, in the cases where the Suzuki-Trotter is exact-namely the different terms of the Hamiltonian commute-there is no need for repetitions. The Suzuki-Trotter approach allows us to simulate any Hamiltonian dynamics as a sequence of unitary operations. The next step is now to decompose any unitary as a sequence of simple quantum gates, typically chosen from a desired set of universal single-qubit and two-qubit quantum gates.

Gate Decomposition
In general, if we find an unitary operation U such that: then we can write the dynamics governed by the Hamiltonian H as: The idea is then to relate H with a simple single-qubit H 0 via an unitary U. Then, this U can be decomposed into a sequence of single-qubit and two-qubit gates of interest, by using well-known techniques [39]. We will analyze some examples of gate decompositions below.

Examples
We will apply the techniques of the previous section to several families of interesting bosonic Hamiltonians.

Boson Sampling and Boson Sampling Hamiltonian
Boson sampling consists of the sampling of the output photon-number statistics of a linear-optics network operating over a number M of photonic modes, which are initialized in a certain Fock state containing n photons. It has been shown [4] that in the regime where M n 2 and n 20 − 50 [40], this problem is computationally hard and should be intractable by classical devices. While this feature turns a boson sampler into a compelling architecture for quantum supremacy, the experiments [5][6][7][8] have not yet reached the aforementioned post-classical regime.
Recently, a Hamiltonian formulation of boson sampling was introduced [41]: where R is the matrix representation of a random unitary transformation and a and b are the input and output modes, respectively. Evolution with this Hamiltonian transforms an initial Fock state: into the N boson sampling superposition: with a photon distribution given by the permanents |γ n | 2 = | vac|b †n 1 1 · · · b †n M M |φ(π/2) | 2 , n i ∈ {0, 1}. Therefore, our task would be to digitalize the Hamiltonian in Equation (9). In this case, we do not need to use the Suzuki-Trotter decomposition since we can directly leverage the Reck decomposition, which is standard in the boson sampling literature and entails that we can decompose the unitary evolution governed by H BS into a mesh of M(M − 1)/2 beam-splitters and appropriate phase-shifters. Each beam-splitter unitary would be given by: Since they are the building blocks of the boson sampling simulation, it is worth analyzing in detail the digitalization of these unitary beam-splitting interactions. Moreover, beam-splitter interactions are also crucial for applications in continuous-variables quantum information processing and quantum computing [11][12][13].

Beam-Splitters
The initial state of boson sampling contains a maximum number of one photon per mode; therefore, the first beam-splitter of the mesh would only require 2 modes × 2 qubits per mode = 4 qubits. Then, it is interesting to consider in detail this case.
Then, we can perform a separate gate decomposition for each U (i) +− . For instance, for the first term, we can use: x σ (2) x σ (3) x σ (4) x .
Then, by using Equations (7), (8), and (21), we find: where: See this gate decomposition in Figure 1. Note that similar decompositions can be obtained for the rest of the U (i) +− 's by adding at the end of the string the number of e iπ/4σ (i) z necessary to rotate some of the σ x to σ y ; although there might be more efficient decompositions for each particular case. With this procedure, we will have a total number of 24 single-qubit rotations and 24 two-qubit ZX-gates.
x , where iand j stand for the qubits among which the two-qubit gate operates, and X = e i π 4 σ (1) x . The rest of the terms U The relevant initial states would be for instance |0 + ⊗ |1 − and |1 + ⊗ |0 − , which via Equations (1) and (15), are mapped to |0110 and |1001 , respectively, which can be obtained by spin-flipping a pair of qubits out of the four-qubit ground state.
Putting all the above together, a single two-mode beam-splitter with one photon per mode can be simulated in a four-qubit quantum simulator. Indeed, we launched the experiment in Figure 1 for the initial state |1 + ⊗ |0 − in IBM Q 5 Tenerife by means of the qiskit tools. An extra step is needed for that, which is the conversion of the ZXgates into the CNOTgates available in the IBM architecture. This can be done by using: (and similarly with 1 − 3 and 1 − 4). Inserting Equation (25) into Equations (23) and (24) and making straightforward simplifications, we get the circuit in Figure 2.
In particular, we choose ε +− 8 = π, so the dynamics should completely transform the initial state into |0 + ⊗ |1 − , namely |0110 . We check the state of the first qubit, which is the one involved in more quantum gates (after running the circuit). We get the right state in 1717 out of 2048 shots, namely 84%. Running the whole beam-splitter-that is, the eight terms of Equation (16)-with the same parameters, the state of the first qubit should be one, and we get the right state 1158 shots out of 2048 (57%). Note that a change in the initial state is not expected to modify the fidelity significantly, since the number of required gates is not significantly affected.

Sequence of Beam Splitters
If we allow more modes and more photons per mode in the simulation, the number of required qubits increases dramatically. However, since a beam-splitter only acts upon two modes and the transition from |n j to |n ± 1 j only flips two qubits, the action of a beam-splitter onto a definite Fock state |n j ⊗ |m k with n photons in mode j and m photons in mode k can still be implemented in a few-qubit quantum simulator. However, this is not enough for a boson-sampling architecture. After a number of beam-splitting steps, the state of the modes would not be a definite Fock state, but a non-trivial combination of them. Therefore, a large number of qubits would be involved at each step, a signature of the computational complexity of boson sampling.
In the post-classical regime, N P 20 − 50 and the number of modes N 2 P . Therefore, the full simulation would require 10 4 − 10 5 qubits, a number that seems completely out of reach.

Two-Mode Squeezing
In a Gaussian boson sampler [42], the initial state is Gaussian, as opposed to the initial Fock state of the standard boson sampling. A particular architecture is a setup in which half of the output of M two-mode squeezers is input into a linear network of M optical modes, while the other half is sent directly to single photon detectors. Then, n single photons are detected in the latter half. As shown in [42], this device is able to solve a randomized version of boson sampling known as scattershot boson sampling, which possesses a similar computational complexity as the original problem [4] and therefore is widely believed to be out of reach classically.
To simulate this, we only need to add to our recipe the ability of simulating an initial set of two-mode squeezed states. For each of two modes, we only need to initialize them in the vacuum state |0101 and then apply the unitary: Restricting ourselves to moderate values of the squeezing parameter β, the probability of generating more than one photon pair would be low, and then, we would be able to assume N P = 1. In particular, the average number of photons per mode would be given by sinh 2 β [43], so for β = 0.5, we have an average number of photons of 0.3, and the approximation is still safe. Note that the approximation of a maximum of one photon per mode is also assumed in scattershot boson sampling [42]. Therefore, again only four qubits per pair of modes would suffice. For each two-mode squeezer, similar techniques as in the case of the beam-splitter can be applied, with similar results in terms of computational complexity. Indeed, for each two-mode squeezer, we get a similar expression to Equation (16), but now instead of having products of even numbers of σ x and σ y , we have products of odd numbers, such as σ (1) x σ (2) x σ (3) x σ (4) y , which give rise to similar terms as in Equations (23) and (24). Moreover, they also commute with each other. In particular, as in the beam-splitter terms, they can be obtained by adding as many e iπ/4σ (i) z as the required σ (i) y to U in Equation (24). See Figure 3 for the gate decomposition corresponding to σ (1) x σ (2) x σ (3) x σ Figure 3. Sequence of gates that implement the two-mode squeezer term σ (1) x σ (2) x σ (3) x σ (4) y .
x , where i and j stand for the qubits among which the two-qubit gate operates, and X = e i π 4 σ (1) Note that larger values of squeezing can already be achieved in other systems, such as cold atom setups [44][45][46][47]. In our case, we would need to increase the number of excitations allowed and therefore the number of qubits. Note also that we are not considering three-mode or multimode squeezing Hamiltonians, which can be achieved for instance in cold atoms [46] and superconducting circuits [48] and which would give rise to completely different multimode states [48,49].

Bogoliubov Transformations
Combining the techniques for digitalizing beam-splitters and two-mode squeezers, we can consider the simulation of arbitrary Bogoliubov transformations, which transform any annihilation operator a i to: From a Hamiltonian viewpoint, the transformation in Equation (27) amounts to a set of two-mode beam-splitters implemented by the α ji coefficients and two-mode squeezers implemented by the β ji coefficients. Therefore, in principle, any Bogoliubov transformation could be simulated using the techniques developed in this paper. Bogoliubov transformations are ubiquitous in quantum field theory, especially in relativistic applications such as the dynamical Casimir effect, the Unruh effect, and Hawking radiation. An immediate application would be a digital quantum simulation of the dynamical Casimir effect, namely the generation of photons out of the vacuum of a quantum field by means of the relativistic motion of a mirror, which generates two-mode squeezing. For the simplest case of two modes and moderate values of the squeezing-which is proportional to the velocity of the moving mirror-four qubits would be enough, as in the case of the previous subsection. Considering higher values of the squeezing and thus higher values of N P would allow considering larger velocities, presumably overcoming the limitations in analog quantum simulators [23]. Note that we only need to go to N p = 2, since the maximum value of squeezing achieved in experiments is 0.46 [23], which can be safely simulated with N p = 1, as discussed above.

Quantum Information Processing and Quantum Computing Gates
There is growing effort, both theoretically and experimentally [11][12][13], in continuous-variable quantum information and quantum computation with multiphoton states in superconducting microwave cavities. As explained for instance in [11], a key step is the ability of generating a controllable set of bilinear interactions, namely an interaction Hamiltonian of the form: which is a combination of beam-splitting and two-mode squeezing terms. This interaction Hamiltonian can be simulated with our techniques. In [13], a Hamiltonian like Equation (28) with only beam-splitting terms was used to generate experimentally several gates, with a fidelity ranging from 60 to 90%. Therefore, it could be of interest to compare these results with a digital quantum simulation, with could reach similar fidelities and benefit from quantum error correction techniques.

Molecular Force Fields
Molecular transitions can be described by using a set of bosonic modes immersed in a force field, which can be modeled by a non-harmonic oscillator [32,33]. For each mode, we will have: where n j = a † j a j is the number operator of the j th mode and the parameter χ j accounts for the anharmonicity of the Hamiltonian. The mapped number operator is [16,17]: and then: Typically, it is enough to consider three or four excitations per mode. Considering for instance N P = 4, we would need five qubits per mode. Then, expanding explicitly Equations (30) and (31), we obtain:

Conclusions
There is a growing interest in building up large networks of thousands of qubits for groundbreaking applications in quantum simulation and quantum computing. In parallel, continuous-variable quantum information processing and quantum computation has motivated a renewed effort in the construction of large linear optical grids. In this work, we brought the latter into the realm of the former, by providing a recipe for the digitalization of bosonic linear and nonlinear optics interactions. Applications range from the dynamical Casimir effect to molecular force fields. A boson sampler in the quantum supremacy regime would require a number of qubits that seems out of reach; however, we have shown that all the individual ingredients of a boson sampler can be digitalized, enabling the digital quantum simulation of a few-mode setup. All these phenomena, which are typically considered in analog simulators, would then benefit in principle from the advantages of digitalization, such as error correction. Moreover, we can digitalize quantum gates that are crucial for continuous-variable quantum information processing and quantum computing. As an example, we have launched several experiments in IBM Q 5 Tenerife, obtaining a fidelity of 84% for the digital quantum simulation of one building block of a beam-splitter and 58% for the full beam-splitter.