*Article* **Quantum Photonic Simulation of Spin-Magnetic Field Coupling and Atom-Optical Field Interaction**

#### **Jesús Liñares \*,†, Xesús Prieto-Blanco †, Gabriel M. Carral † and María C. Nistal †**

Quantum Materials and Photonics Research Group, Optics Area, Department of Applied Physics, Faculty of Physics/Faculty of Optics and Optometry, Campus Vida s/n, University of Santiago de Compostela, E-15782 Santiago de Compostela, Galicia, Spain; xesus.prieto.blanco@usc.es (X.P.-B.); gabrielmaria.carral@rai.usc.es (G.M.C.); mconcepcion.nistal@usc.es (M.C.N.)

**\*** Correspondence: suso.linares.beiras@usc.es

† These authors contributed equally to this work.

Received: 8 November 2020; Accepted: 7 December 2020; Published: 10 December 2020

**Abstract:** In this work, we present the physical simulation of the dynamical and topological properties of atom-field quantum interacting systems by means of integrated quantum photonic devices. In particular, we simulate mechanical systems used, for example, for quantum processing and requiring a very complex technology such as a spin-1/2 particle interacting with an external classical time-dependent magnetic field and a two-level atom under the action of an external classical time-dependent electric (optical) field (light-matter interaction). The photonic device consists of integrated optical waveguides supporting two collinear or codirectional modes, which are coupled by integrated optical gratings. We show that the single-photon quantum description of the dynamics of this photonic device is a quantum physical simulation of both aforementioned interacting systems. The two-mode photonic device with a single-photon quantum state represents the quantum system, and the optical grating corresponds to an external field. Likewise, we also present the generation of Aharonov–Anandan geometric phases within this photonic device, which also appear in the simulated systems. On the other hand, this photonic simulator can be regarded as a basic brick for constructing more complex photonic simulators. We present a few examples where optical gratings interacting with several collinear and/or codirectional modes are used in order to illustrate the new possibilities for quantum simulation.

**Keywords:** integrated photonics; quantum optics; quantum simulation

#### **1. Introduction**

One of the most promising tasks in quantum science and technology is the implementation of quantum simulations. Its physical foundation is based on the fact that the dynamics of a quantum system is governed by its Hamiltonian *H*ˆ (time evolution) or momentum operator *M*ˆ (spatial evolution), that is given a Hilbert space H and some input state |Ψ(0)∈H, the full evolution of the system is given by the action of the evolution operator on such a state. The evolution operator can be either the time evolution one *<sup>U</sup>*ˆt <sup>=</sup> exp (−*iHt* <sup>ˆ</sup> /¯*h*), which comes from the Schro¨dinger equation <sup>−</sup>*ih*¯ *<sup>∂</sup>*|Ψ/*∂<sup>t</sup>* <sup>=</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>|</sup>Ψ, that is <sup>|</sup>Ψ(*t*) <sup>=</sup> *<sup>U</sup>*ˆt|Ψ(0), or its spatial counterpart *<sup>U</sup>*<sup>ˆ</sup> <sup>s</sup> <sup>=</sup> exp (*iMs* <sup>ˆ</sup> /¯*h*), which comes from the momentum operator in the position representation, that is *ih*¯ *<sup>∂</sup>*|Ψ/*∂<sup>s</sup>* <sup>=</sup> *<sup>M</sup>*<sup>ˆ</sup> <sup>|</sup>Ψ, where *<sup>s</sup>* is the spatial variable that defines the direction along which the system evolves (it can be, for instance, the *z*-direction) then <sup>|</sup>Ψ(*s*) <sup>=</sup> *<sup>U</sup>*<sup>ˆ</sup> <sup>s</sup>|Ψ(0). One is perhaps more familiar with the temporal case <sup>|</sup>Ψ(*t*) <sup>=</sup> *<sup>U</sup>*ˆt|Ψ(0), but note that the spatial case is totally analogous [1,2]. Now, as it occurs in nature that very different and unrelated quantum systems share analogous Hamiltonians or momentum operators, the dynamics of these systems will be analogous. This implies that if we have some system of interest, we will be able to mimic its behavior (evolution) by means of other very different system. This last system, which is

required to be fully controllable by the experimenter, constitutes a quantum simulator [3,4] of the system of interest.

On the other hand, the importance of making quantum simulations obeys various reasons. First of all, whether classical or quantum, simulation is often an indispensable tool in science. Physical systems can be very complicated to study. If we find a way to mimic their behavior in a manageable way, their dynamics can be analyzed with much less difficulty, much less expensively, and much faster [3]. For instance, we can imagine that the system of interest is one whose working conditions are fragile and very sensitive to small perturbations. It could also be that those conditions are very specific and/or hard (even impossible) to achieve in the laboratory. Simulating these complicated systems will, in principle, improve our knowledge about them with much ease. Secondly, quantum simulation offers a crucial advantage over classical simulation. Quantum systems have a greater information storing capacity than classical systems; thus, a quantum simulator, which is a quantum system itself, will require much fewer physical resources than a classical simulator to store the same amount of information. Furthermore, classical simulations find difficulties simulating quantum systems when strong entanglement is present [3]. Finally, it is possible that the potential applications of a quantum system may be better implemented using a quantum simulation [5], if it happens that the quantum simulation parameters are easier to handle, compared to those of the quantum system of interest, in the sense of greater tunability [4]. All of this strongly suggests the need for quantum simulations. Formally, one has a quantum device, called the quantum simulator, which reproduces the evolution of the system of interest, the quantum system. Following [4], let us call the initial input state of the quantum system |*φ*(0), which evolves to |*φ*(*ζ*), where *ζ* stands either for time or some spatial coordinate, under the evolution operator *<sup>U</sup>*<sup>ˆ</sup> <sup>=</sup> exp (*iO*<sup>ˆ</sup> *<sup>ζ</sup>*/¯*h*). Here, *<sup>O</sup>*<sup>ˆ</sup> can be a Hamiltonian (−*H*<sup>ˆ</sup> ) or a momentum operator (*M*<sup>ˆ</sup> ). The quantum simulator, on the other hand, starts in the state <sup>|</sup>*ψ*(0) and evolves to <sup>|</sup>*ψ*(*ζ*) under the action of the operator *<sup>U</sup>*<sup>ˆ</sup> <sup>=</sup> exp (*iO*<sup>ˆ</sup> *ζ*/¯*h*). Since we have a system and a simulator, there exists a correspondence relating these elements, that is a correspondence between <sup>|</sup>*φ*(0)↔|*ψ*(0), <sup>|</sup>*φ*(*ζ*)↔|*ψ*(*ζ*), *<sup>U</sup>*<sup>ˆ</sup> <sup>↔</sup> *<sup>U</sup>*<sup>ˆ</sup> , and thus, *<sup>O</sup>*<sup>ˆ</sup> <sup>↔</sup> *<sup>O</sup>*<sup>ˆ</sup> . This is revealed via measurements on both systems. The greater the accuracy of this correspondence, the more one can trust the simulator [3].

In this work, we present an integrated quantum photonic simulator for atom-field quantum interacting systems. It is based on optical gratings and can be regarded as a basic brick for constructing more complex photonic simulators. First of all, we would like to stress the high interest in quantum photonic simulators; for instance, integrated photonic structures allow simulating a number of condensed matter effects such as Anderson localization, Mott transition, etc. [6], and more recently, anyonic interaction has been simulated by using an array of channel optical waveguides with a helically bent axis [7]. In our case, we will simulate in an integrated photonic device the dynamical and topological properties of two very well-known quantum interacting systems. This is relevant since a way to check the fidelity of a simulation is testing the same physical model of the dynamics with different systems and then comparing the results [5]. The two systems to be simulated are the well-known interaction of a spin-1/2 particle under a time-dependent external classical magnetic field (see for instance [8]) and the interaction of a two-level atom with an external classical electric field (see for instance [9]) which are used, at present, for implementing quantum information devices for quantum processing, quantum computation, and so on [8]. It is also well known that these mechanical systems require a very complex technology. We will take a standard photonic device, that is a two-mode planar guide modulated by an integrated grating (see for instance [10]) in which a single-photon quantum state is propagated. We will show that it can emulate both the dynamical and geometrical properties of the two aforementioned systems. As mentioned, this photonic simulator based on integrated optical gratings can be regarded as a basic brick for more complex photonic simulators; thus, the results obtained can be extended to multi-level systems or to several concatenated two-mode systems simulating concatenated temporal operations. Addressing practical issues, this simulator will also benefit from some advantages of classical and quantum integrated photonics [6,11,12], such as the fast and energetically efficient operation and miniaturization capabilities of integrated photonics, which favor scalability. Moreover, integrated photonic devices are becoming some of the most powerful and appealing technologies for classical and quantum information [13,14]. Finally, we must stress that although we present a photonic device for quantum simulation, it can be also used to implement quantum operations and/or effects with the photonic device itself such as logic gates, geometric phases, and so on.

The plan of the paper is as follows: In Section 2, we briefly present the Hamiltonian of the quantum interacting systems, in a suitable form for their photonic simulation, along with the physical parameters and some useful solutions. In Section 3, starting from quantum states as single-photon states, an integrated photonic device for the quantum simulation of both a spin-magnetic field coupling and light-matter interaction is presented, along with the limits of this simulation. In Section 4, quantum simulation of geometric phases is studied. Finally, a summary is presented in Section 5.

#### **2. Mechanical Interacting Systems**

In this section, we present the main results about the mechanical interacting systems whose dynamic and topological properties are going to be simulated with an integrated photonic device. The primary aim is to write the Hamiltonians of these systems in a suitable form to facilitate the study of their photonic simulations.

#### *2.1. Spin-1/2 Particle Interacting with a Magnetic Field*

The Hamiltonian of a spin particle interacting with a magnetic field is given by *<sup>H</sup>*<sup>ˆ</sup> <sup>=</sup> <sup>−</sup>*<sup>μ</sup>* · *<sup>B</sup>*, where *μ* is the magnetic moment of the particle, which is proportional to the spin operator, and *B* is the magnetic field. Let us consider a spin- <sup>1</sup> <sup>2</sup> particle, then *<sup>μ</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup>*μσ*, where *μ* = *h*¯ *γp*, with *γ<sup>p</sup>* the gyromagnetic ratio of the particle (*p* = e for an electron, *p* = n for a neutron, etc.), and *σ* = (*σx*, *σy*, *σz*) are the Pauli matrices. Moreover, we choose a time-dependent magnetic field, which, by assuming a sinusoidal dependence in time with frequency *<sup>ω</sup>*, is written as *<sup>B</sup>* <sup>=</sup> *<sup>B</sup>*0(sin *<sup>θ</sup>* cos *<sup>ω</sup>t*, <sup>−</sup> sin *<sup>θ</sup>* sin *<sup>ω</sup>t*, cos *<sup>θ</sup>*), representing a magnetic field of modulus *B*<sup>0</sup> rotated an angle *θ* with respect to the *z*-axis and spinning around this same axis with a frequency *ω*. We can then write the well-known spin-magnetic Hamiltonian as (see for instance [8]):

$$\hat{H}(t) = -\frac{\hbar}{2}\gamma\_p B\_0 \cos\theta \sigma\_z - \frac{\hbar}{2}\gamma\_p B\_0 \sin\theta (\cos\omega t \,\sigma\_x - \sin\omega t \,\sigma\_y). \tag{1}$$

By taking into account the expressions of the Pauli matrices, then the above Hamiltonian gives rise to the following time-dependent Schrodinger equation: ¨

$$i\hbar \frac{\partial}{\partial t} |\Psi(t)\rangle = -\frac{\hbar}{2} \gamma\_p B\_0 \begin{pmatrix} \cos\theta & \sin\theta \exp\left(i\omega t\right) \\ \sin\theta \exp\left(-i\omega t\right) & -\cos\theta \end{pmatrix} |\Psi(t)\rangle = $$

$$\equiv \hat{H}(t) |\Psi(t)\rangle. \tag{2}$$

This time-dependent Hamiltonian is enough for our simulation purposes. Note that for the neutron case, we have the well-known Nuclear Magnetic Resonance (NMR). The solution of the equation above can be obtained as a time-dependent linear combination of down and up spin states, |0 and |1, of the particle, that is |Ψ(*t*) = *co*(*t*)|0 + *c*1(*t*)|1. Thus, by using the vector representation of |Ψ(*t*) and the matrix Hamiltonian given by Equation (2), the Schro¨dinger equation can be rewritten as follows:

$$i\hbar \frac{d\mathbf{c}\_n(t)}{dt} = E\_n \cos\theta \mathbf{c}\_n(t) + i \sum\_{n \neq m} \mathbb{C}\_{nm}(t) \mathbf{c}\_m(t). \tag{3}$$

with *<sup>n</sup>* <sup>=</sup> 0, 1, *Eo* <sup>=</sup> <sup>−</sup>*h*¯ *<sup>γ</sup><sup>p</sup>* <sup>2</sup> *Bo* and *<sup>E</sup>*<sup>1</sup> <sup>=</sup> *<sup>h</sup>*¯ *<sup>γ</sup><sup>p</sup>* <sup>2</sup> *Bo* the eigenvalues of spin states when a static magnetic field along *z*-direction is applied, and *C*01(*t*) = *C*<sup>∗</sup> <sup>10</sup>(*t*) = <sup>−</sup>*h*¯ *<sup>γ</sup><sup>p</sup>* <sup>2</sup> *Bo* sin *θ* exp(*iωt*) the coupling coefficients. On the other hand, it is interesting, for the sake of expositional convenience, to present the main results about the dynamics of the system. First of all, in order to simplify things, we can shift to a reference frame that is rotating at a frequency *ω* by means of the following unitary transformation:

$$|\eta(t)\rangle = \mathcal{U}(t,\omega)|\Psi(t)\rangle = \exp\left(-i\omega\sigma\_z t/2\right)|\Psi(t)\rangle;\tag{4}$$

therefore, by inserting |Ψ(*t*) = exp (*iωσzt*/2)|*η*(*t*) into the Schro¨dinger equation, we obtain, after some simple calculations, the following Hamiltonian acting on |*η*(*t*):

$$
\langle \hat{H}' | \eta(t) \rangle = \frac{\hbar}{2} \begin{pmatrix} -\gamma\_p B\_0 \cos \theta + \omega & -\gamma\_p B\_0 \sin \theta \\ -\gamma\_p B\_0 \sin \theta & \gamma\_p B\_0 \cos \theta - \omega \end{pmatrix} | \eta(t) \rangle. \tag{5}
$$

This Hamiltonian *<sup>H</sup>*ˆ can be rewritten as a product of the form <sup>−</sup>*h*¯ <sup>2</sup>*γp*(*<sup>σ</sup>* · *<sup>B</sup>e f*), with *<sup>B</sup>e f* an effective static magnetic field, that is,

$$\hat{H}' = -\frac{\hbar}{2}\gamma\_p(\boldsymbol{\sigma} \cdot \mathbf{B}\_{\ell f}) = -\frac{\hbar}{2}\gamma\_p B\_0 \Delta \cos\theta' \sigma\_z - \frac{\hbar}{2}\gamma\_p B\_0 \Delta \sin\theta' \sigma\_x = -\frac{\hbar}{2}\Delta\_0 \begin{pmatrix} \cos\theta' & \sin\theta'\\ \sin\theta' & -\cos\theta' \end{pmatrix}, \tag{6}$$

with:

$$
\Delta\_o = \gamma\_p B\_0 \Delta = \gamma\_p B\_o' \tag{7}
$$

and:

$$
\Delta = \gamma\_p B\_0 \sqrt{1 - \frac{2\omega}{\gamma\_p B\_0} \cos \theta + \frac{\omega^2}{\gamma\_p^2 B\_0^2}};\tag{8}
$$

therefore, the effective magnetic field is given by *Be f* = (*B* <sup>0</sup> sin *θ* , 0, *B* <sup>0</sup> cos *θ* ), where *B* <sup>0</sup> and *θ* are related to *B*<sup>0</sup> and *θ* by *B* <sup>0</sup> = *B*0Δ, with sin *θ* = *γpB*<sup>0</sup> sin *θ*/Δ*<sup>o</sup>* = sin *θ*/Δ and cos *θ* = *γpB*0[cos *θ* − (*ω*/*γpB*0)]/Δ*<sup>o</sup>* = [cos *θ* − (*ω*/*γpB*0)]/Δ. As seen from the expression for Δ, it would seem that for this reparametrization to have physical meaning, restrictions on the parameters would appear, as Δ contains a possible non-positive term under the square root. The worst case would be likely to happen if cos *<sup>θ</sup>* <sup>=</sup> 1, but one can find that the resulting quantity 1 <sup>−</sup> (2*h*¯ *<sup>ω</sup>*/*μB*0)+(*h*¯ *<sup>ω</sup>*/*μB*0)<sup>2</sup> is never negative, for any value of the parameters.

Next, it is interesting to compute the eigenstates <sup>|</sup>*η* and eigenvalues *<sup>E</sup><sup>η</sup>* of Hamiltonian *<sup>H</sup>*ˆ given by Equation (6), that is <sup>|</sup>*η*(*t*) <sup>=</sup> exp − i *<sup>h</sup>*¯ *Eηt* |*η*(0). After a standard calculation, we have:

$$|\eta\_{+}(0)\rangle = \cos\frac{\theta'}{2}|0\rangle + \sin\frac{\theta'}{2}|1\rangle, \quad |\eta\_{-}(0)\rangle = \sin\frac{\theta'}{2}|0\rangle - \cos\frac{\theta'}{2}|1\rangle. \tag{9}$$

where <sup>|</sup>0 ≡ (1, 0)*<sup>T</sup>* and <sup>|</sup>1 ≡ (0, 1)*T*, with *<sup>T</sup>* denoting the transpose. The eigenvalues are *<sup>E</sup>*<sup>±</sup> <sup>=</sup> <sup>±</sup>*h*¯ <sup>2</sup>Δ*o*, that is *<sup>E</sup>*<sup>±</sup> <sup>=</sup> <sup>±</sup>*h*¯ 2*γpB* <sup>0</sup> <sup>=</sup> <sup>±</sup>*<sup>μ</sup>* <sup>2</sup> *B* <sup>0</sup>. Therefore, by taking into account Equation (4), the full evolved states |Ψ(*t*) can be easily obtained.

#### *2.2. Two-Level Atom Interacting with an Electric (Optical) Field*

On the other hand, let us consider a two-level atom coupled to a harmonic external classical electric field, which describes semiclassical light-matter interaction. As is usually done [9], we label the atomic levels as |*g* (ground) and |*e* (excited). They have energies *h*¯ *ω*<sup>g</sup> and *h*¯ *ω*e, respectively. Their energy difference is given by *hω*<sup>0</sup> = *h*¯ (*ω*<sup>e</sup> − *ω*g). The expression for this non-interacting part of the Hamiltonian can be formally written as follows:

$$
\hat{H}\_o = \hbar\omega\_{\tilde{\mathbf{g}}}|\mathbf{g}\rangle\langle\mathbf{g}| + \hbar\omega\_{\mathbf{e}}|\mathbf{e}\rangle\langle\mathbf{e}| = \hbar\bar{\omega}\mathbf{I} + \hbar\omega\_o\sigma\_z/2,\tag{10}
$$

where we used the matrix representation of |*gg*| and |*ee*|, with *ω*¯ = (*ω*<sup>g</sup> + *ω*e)/2, and **I** the two-dimensional identity matrix. As for the interacting part, it is given by the dipole interaction (electric dipole approximation) between an external electric (optical) field and the atom. Indeed, by assuming, for the sake of simplicity, that the field is propagating along *z*, has a sinusoidal dependence in time with frequency *ω*, and is linearly polarized along the *x*-direction, then the electric (optical) field can be written as follows *E* = *E*<sup>0</sup> cos *ω*(*z*/*c* − *t ux*. Moreover, we disregard the spatial dependence of the electric (optical) field because the wavelength *λ* = *c*/2*πω* is considered much larger than the atomic dimensions. Therefore, we apply the dipole or long-wavelength approximation, that is by assuming without lost of generality *ωz*/*c* = 2*mπ*, with *m* an integer, then *E* = *E*<sup>0</sup> cos *ωt ux*, the electric-dipole interaction is given by *<sup>H</sup>*<sup>ˆ</sup> *<sup>I</sup>* <sup>=</sup> <sup>−</sup>*<sup>d</sup>* · *<sup>E</sup>* <sup>=</sup> <sup>−</sup>*dxEx* <sup>=</sup> <sup>e</sup>*xEx*, where *<sup>d</sup>* is the atomic dipole operator *<sup>d</sup>* <sup>=</sup> *<sup>q</sup>* · *<sup>r</sup>*, with *<sup>q</sup>* <sup>=</sup> <sup>−</sup>e, and *<sup>r</sup>* is the electron's position vector (operator). This interaction term allows transitions between the two levels. The form of the dipole operator can be calculated by using twice the closure relationship <sup>ˆ</sup>*<sup>I</sup>* <sup>=</sup> <sup>|</sup>*gg*<sup>|</sup> <sup>+</sup> <sup>|</sup>*ee*|, that is <sup>ˆ</sup>*I*(−e*x*)ˆ*I*; therefore:

$$d\_{\mathbf{x}} = -\langle \mathbf{g} | \mathbf{ex} | \mathbf{g} \rangle \langle \mathbf{g} \rangle \langle \mathbf{g} | - \langle \mathbf{e} | \mathbf{ex} | \mathbf{e} \rangle \langle \mathbf{e} \rangle \langle \mathbf{e} | - \langle \mathbf{g} | \mathbf{ex} | \mathbf{e} \rangle (|\mathbf{g}\rangle \langle \mathbf{e}| + |\mathbf{e}\rangle \langle \mathbf{g}|) \equiv -\mathbf{c} + \mathbf{I} + \mathbf{c} - \mathbf{c}\_{2} - d\_{0} \sigma\_{\mathbf{x}}.\tag{11}$$

with *c*<sup>±</sup> = (*g*|e*x*|*g*±*e*|e*x*|*e*)/2, *d*<sup>0</sup> = *g*|e*x*|*e*, and where we have used the matrix representations (|*ge*| + |*eg*| ≡ *σx*, |*gg*| ≡ (**I** + *σz*)/2, and |*ee*| = (**I** − *σz*)/2. We must stress that sometimes, parity arguments [9] narrow the form of the dipole operator, leading up to the expression *dx* = −*g*|*ex*|*e*(|*ge*| + |*eg*|) = −*doσx*. Likewise, it is customary to introduce the Rabi frequency Ω = *Eodo*/¯*h*. In short, the total Hamiltonian is given by *H*ˆ *<sup>o</sup>* + *H*ˆ *<sup>I</sup>*, and therefore, the corresponding Schrodinger equation, by using this new variables, is given by: ¨

$$i\hbar\frac{\partial}{\partial t}|\Psi(t)\rangle = \left(\hbar(\bar{\omega}+\mathbb{C}\_{+})\,\mathrm{I}-\hbar(\frac{\omega\_{0}}{2}-\mathbb{C}\_{-})\sigma\_{z}+\hbar\Omega\sigma\_{x}\cos\omega t\right)|\Psi(t)\rangle = \hat{H}|\Psi(t)\rangle,\tag{12}$$

with *C*<sup>±</sup> = *c*<sup>±</sup> cos *ωt*. The general solution of this Schro¨dinger equation can again be obtained by using a time-dependent linear combination of fundamental and excited states, |*g*≡|0 and |*e*≡|1, of the particle, that is |Ψ(*t*) = *co*(*t*)|0 + *c*1(*t*)|1. However, the high value of the frequency *ω* means that the electric field is rapidly oscillating, which suggests to make the following change |Ψ(*t*) = *f*0(*t*) exp(*iωt*/2)|0 + *f*1(*t*) exp(−*iωt*/2)|1 ≡ exp (*iωσzt*/2)|*η*(*t*). Moreover, in many cases, by parity arguments, it is fulfilled that *c*<sup>±</sup> = 0. Therefore, by substituting this state into Equation (12), we obtain the following Hamiltonian for the state *f*0(*t*), *f*1(*t*) *<sup>T</sup>* ≡ |*η*(*t*),

$$
\langle \hat{H}' | \eta(t) \rangle \approx \left( \hbar \bar{\omega} \, \mathbf{I} - \frac{\hbar \delta}{2} \sigma\_z + \frac{\hbar \Omega}{2} \sigma\_x \right) | \eta(t) \rangle,\tag{13}
$$

with *δ* = (*ω*<sup>0</sup> − *ω*) the detuning parameter and where we have neglected terms rapidly oscillating of the form exp(±*iωt*); that is, a temporal Rotating Wave Approximation (RWA) has been made, and a time independent Hamiltonian has thus been obtained. We must stress that under this approximation, the terms *C*<sup>±</sup> in Equation (12) can be also neglected independently of the wave-function parity. Finally, note that Hamiltonians given by Equations (6) and (13) have the same algebraic structure.

#### **3. Quantum Photonic Simulations**

In this section, we present the quantum simulation, which can be implemented by an integrated photonic device, that is integrated optical gratings supporting two collinear guided modes, that is two mode guides assisted by a periodic perturbation. It can be considered the basic brick for constructing more complex simulators.

#### *3.1. Classical Study of the Photonic Device*

Let us consider a standard integrated photonic device consisting of, for example, an integrated waveguide 1D (one-dimensional) with refractive index *n*(*x*) (slab guide) that supports two optical modes *e*0(*x*) and *e*1(*x*). These modes are collinear and travel in the *z*-direction with propagation constants *β*<sup>0</sup> and *β*1, that is wave vectors in the *z*-direction defined as *β* = (*ω*/*c*)*N*, where *ω* is the frequency of the mode and *N* is the effective index [10]. The mentioned modes are coupled by an integrated grating present in a region of the slab guide, as shown in Figure 1a), where relevant parameters are indicated, that is substrate index *ns*, film index *nf* , cover index *nc* = 1, and film thickness *d*. The grating is represented, for example, by a periodic modulation (perturbation) of the electrical permittivity [10],

$$
\Delta\varepsilon(\mathbf{x}, z) = \Delta\varepsilon(\mathbf{x})\cos(\gamma z + \alpha\_o),
\tag{14}
$$

with Δ(*x*) the modulation strength of the optical grating, *α<sup>o</sup>* an initial phase, if required, *γ* = 2*π*/Λ the frequency of the perturbation, and Λ its period. This index profile can be obtained by different technologies of integrated optics, for instance ion-exchange in glass [15,16] could be used, or even by optical fiber technology [17]. Likewise, in crystals such as lithium niobate, these optical gratings can also be reconfigurable due to acousto-optic or electro-optic effects [18].

**Figure 1.** (**a**) Integrated optical grating with period Λ on a two-mode slab waveguide with modes *e*0(*x*) and *e*1(*x*). The inset shows the simulated mechanical device: atom-field interaction. (**b**) Integrated optical grating with period Λ on a two-mode channel waveguide with modes, for example *e*00(*x*, *y*) and *e*10(*x*, *y*).

As we assume that only two guided modes are excited in our integrated photonic structure, the perturbed electric field amplitude is given by *e*(*x*, *z*) = *a*0(*z*)*e*0(*x*) + *a*1(*z*)*e*1(*x*). It is well known that the general set of equations that describe the coupling between *n* copropagating optical modes with propagation constants *β<sup>n</sup>* in a perturbed waveguide and that allow calculating the amplitude coefficients *an*(*z*) is given by (see for instance [10]):

$$-i\frac{da\_n(z)}{dz} = \beta\_n(z)a\_n(z) + \sum\_{n \neq m} \mathbb{C}\_{nm}(z)a\_m(z),\tag{15}$$

where *β*˜*n*(*z*) = *β<sup>n</sup>* + *Cnn*(*z*) are the corrected propagation constants due to the self-coupling coefficients *Cnn* and *Cnm* are the coupling coefficients between the *n* mode and each of the other *m* modes. All these coefficients are calculated as follows [10]:

$$\mathbb{C}\_{\text{nW}}(z) = \frac{\omega}{2} \int \Delta \epsilon(\mathbf{x}, z) \, e\_n(\mathbf{x}) \, e\_m^\*(\mathbf{x}) \, \text{d}\mathbf{x},\tag{16}$$

with *ω* the temporal frequency of the modes and *en*(*x*) and *em*(*x*) the normalized optical *n* and *m* modes of a planar guide.

The study with 2D guides (integrated channel guides or optical fibers) can be also made, but no new relevant result would be obtained. Indeed, a channel guide can be defined starting from a planar guide whose width is reduced up to a size *a* of the same order as its depth, that is *a* ≈ *d*, as shown in Figure 1b). In such a case, the optical modes are characterized by two subscripts, one for each spatial direction, that is *enp*(*x*, *y*) and *emq*(*x*, *y*); therefore, the general coupling coefficients are given by:

$$\mathbb{C}\_{\text{npm}\eta}(z) = \frac{\omega}{2} \int \Delta \epsilon(\mathbf{x}, \mathbf{y}, z) \, e\_{\text{np}}(\mathbf{x}, \mathbf{y}) \, e\_{\text{mq}}^\*(\mathbf{x}, \mathbf{y}) \, \text{dxdy},\tag{17}$$

where we assumed that the grating modulation can also have a *y*-dependence. Finally, the mode coupling equations can be obtained by applying the formal changes *n* → *np*, *m* → *mq* in Equation (15). Accordingly, the results for planar guides can be easily transferred to channel guides.

#### *3.2. Quantum Study of the Photonic Device*

A canonical quantization procedure [1,2] proves that the *a*(*z*) coefficients become the photon absorption (or emission) operators *a*ˆ(*z*) (correspondence principle). Therefore, the coupled mode equations are in fact the Heisenberg equations, which, in general, give the evolution of the operators in time. In this case, they give the spatial evolution of the operators. Moreover, the relevant operator here responsible for the spatial propagation of quantum states of the device is the momentum operator, which is the generator of spatial translations [2,19], and not the Hamiltonian, which is the generator of temporal translations. In short, we can study the integrated photonic device in a fully quantum mechanical way, by solving the equations for absorption operators (spatial Heisenberg equations). That is, by performing the change *a*(*z*) → *h*¯ *a*ˆ(*z*) in Equation (15), we obtain [2,20]:

$$-i\hbar \frac{d\mathfrak{d}\_n(z)}{dz} = \hbar \mathfrak{J}\_n \mathfrak{d}\_n(z) + \hbar \sum\_{n \neq m} \mathbb{C}\_{nm}(z) \mathfrak{d}\_m(z). \tag{18}$$

We must stress that modal coupling preserves energy; therefore, Equation (18) corresponds to a unitary transformation, and accordingly, *Cnm* = *Cmn*. On the other hand, we only consider single-photon states, which is enough for our simulation purposes. We must stress that the linear momentum of a single photon without modal coupling, that is Δ(*x*, *y*) = 0, is given by *p*(0,1) = *h*¯ *β*(0,1) depending on whether the photon is excited in mode *β*<sup>0</sup> or mode *β*<sup>1</sup> [19,20]. Obviously, more general quantum states could be used such as multiphoton states, entangled states of two photons, and so on. These states would give rise to more complex quantum simulations, which fall outside of the scope of this work. In general, the single-photon state is a quantum superposition because the photon can either be excited in the mode *β*0, that is |10, or in the mode *β*1, that is |11. Hence, the general quantum state is given by:

$$|L(z)\rangle = a\_0(z)|1\_0\rangle + a\_1(z)|1\_1\rangle,\tag{19}$$

where *a*0(*z*) and *a*1(*z*) are the quantum complex amplitudes and fulfil the normalization condition |*a*0(*z*)| <sup>2</sup> <sup>+</sup> <sup>|</sup>*a*1(*z*)<sup>|</sup> <sup>2</sup> <sup>=</sup> 1. States <sup>|</sup>10 and <sup>|</sup>11 must be understood as single-photon states at a distance (plane) *z*. It can be checked that the solutions of Equation (18) for the spatial propagation of emission operators are the same as for single-photon states |*L*(*z*) [21]. The main reason is that single-photon states are proportional to emission operators, that is <sup>|</sup>10 <sup>=</sup> *<sup>a</sup>*ˆ† <sup>0</sup>|0, <sup>|</sup>11 <sup>=</sup> *<sup>a</sup>*ˆ† <sup>1</sup>|0. Indeed, let us consider free propagation, that is non-coupling case *Cnm* = 0, then the solutions of Equation (18) are *a*ˆ0,1(*z*) = exp(*iβ*0,1*z*)*a*ˆ0,1(0). Now, let us consider, for the sake of simplicity, a single-photon state at *z* = 0, for instance |10, then the optical propagation can be obtained as follows: <sup>|</sup>10 <sup>=</sup> *<sup>a</sup>*ˆ† <sup>0</sup>(0)|0 (one-photon emission), but by taking into account the *z*-propagation, we have <sup>|</sup>*L*(*z*) <sup>=</sup> exp(*iβ*0*z*)*a*ˆ† <sup>0</sup>(*z*)|0 = exp(*iβ*0*z*)|10 ≡ *a*0(*z*)|10; therefore, the coefficients *a*0,1(*z*) of a single-photon state have the same optical propagation solution as the operators *a*ˆ0,1(*z*). This can be proven for a general unitary transformation after a certain algebra. We will take advantage of this property for simulation purposes. Accordingly, for single-photon states, Equation (18) can formally be rewritten as follows:

$$i\hbar \frac{\mathbf{d}|\mathbf{L}(\mathbf{z})\rangle}{\mathbf{d}\mathbf{z}} = -\hbar \begin{pmatrix} \vec{\beta}\_0(z) & \mathbf{C}\_{12}(z) \\ \mathbf{C}\_{21}(z) & \vec{\beta}\_1(z) \end{pmatrix} |\mathbf{L}(z)\rangle = -\hat{M}(z)|\mathbf{L}(z)\rangle,\tag{20}$$

where <sup>|</sup>*L*(*z*) is given, in vector representation, by (*ao*(*z*), *<sup>a</sup>*1(*z*))*<sup>T</sup>* and *<sup>M</sup>*<sup>ˆ</sup> is the matrix representation of the so-called momentum operator. This is equivalent to the matrix equation given by Equation (2) for the Hamiltonian. We must stress that Equations (2), (3), (13), (18) and (20) are the main results for implementing the quantum simulations in this work.

#### *3.3. Photonic Simulation of Spin-Magnetic Field Interaction*

Let us consider an integrated optical grating characterized by the function given by Equation (14). It will be useful to work with the slowly varying operators *<sup>A</sup>*<sup>ˆ</sup> *<sup>n</sup>*, defined as *<sup>A</sup>*<sup>ˆ</sup> *<sup>n</sup>* <sup>=</sup> *<sup>a</sup>*ˆ*<sup>n</sup>* exp(−*iβnz*). On the other hand, the coupling coefficients for a grating with initial phase *α<sup>o</sup>* = 0 can be written as *C*<sup>00</sup> = *c*<sup>00</sup> cos *γz*, *C*<sup>11</sup> = *c*<sup>11</sup> cos *γz*, and *C*<sup>01</sup> = *C*<sup>10</sup> = *Co* cos *γz*, where *cnm* = (*ω*/2) <sup>3</sup> <sup>Δ</sup>(*x*)*en*(*x*) · *e*∗ *<sup>m</sup>*(*x*)*dx*, and so on. Therefore, from Equation (18), we obtain, for the two-mode case, the following spatial Heisenberg equations:

$$i\hbar \frac{d\hat{\mathbf{A}}\_0(\mathbf{z})}{d\mathbf{z}} = -\hat{A}\_0 c\_{00} \cos \gamma z - \hat{A}\_1 \frac{c\_0}{2} e^{-i\Lambda \beta\_{01} z} [e^{i\gamma z} + e^{-i\gamma z}],\tag{21}$$

$$i\hbar \frac{d\hat{\mathbf{A}}\_1(\mathbf{z})}{d\mathbf{z}} = -\hat{A}\_1 c\_{11} \cos \gamma z - \hat{A}\_0 \frac{c\_0}{2} e^{i\Lambda \beta\_{01} z} [e^{i\gamma z} + e^{-i\gamma z}],\tag{22}$$

where Δ*β*<sup>01</sup> = *β*<sup>0</sup> − *β*1. The above equation reveals that we have oscillating terms coming from the cosine of the self-coupling terms with arguments ±*γz* and also terms with arguments (*γ* − Δ*β*01)*z* and (*γ* + Δ*β*01)*z*. We make the assumption that *γ* is of the same order as Δ*β*01, so that (*γ* − Δ*β*01) is small. The other terms are rapidly oscillating and, thus, will average to zero on a sufficiently large *z*-scale. We must stress that what we do here is essentially a spatial RWA, which is well-known in light-matter interaction in the time domain. Therefore, the above equations become, to a good approximation,

$$i\frac{d\boldsymbol{A}\_0(z)}{dz} = -\hat{A}\_1 \frac{\mathbb{C}\_o}{2} \exp\left[-i(\Delta\beta\_{01} - \gamma)z\right],\tag{23}$$

$$i\frac{d\mathbf{A}\dot{A}\_1(z)}{dz} = -\dot{A}\_0 \frac{\mathbb{C}\_o}{2} \exp\left[i(\Delta\beta\_{01} - \gamma)z\right].\tag{24}$$

Next, we make the following relabeling Δ*β*<sup>01</sup> ≡ Δ*β* and define the new absorption operators *A*ˆ <sup>0</sup> = ˆ *<sup>b</sup>*<sup>0</sup> exp(−*i*Δ*βz*/2) and *<sup>A</sup>*<sup>ˆ</sup> <sup>1</sup> <sup>=</sup> <sup>ˆ</sup> *b*<sup>1</sup> exp(*i*Δ*βz*/2). We rewrite the Heisenberg equations in terms of the ˆ *b*(*z*) operators,

$$i\frac{d\mathbf{d}\hat{b}\_0}{\mathbf{d}z} = -\frac{\Delta\beta}{2}\hat{b}\_0 - \frac{\mathbb{C}\_o}{2}\hat{b}\_1 \exp\left(i\gamma z\right),\tag{25}$$

$$i\frac{d\mathbf{d}\hat{b}\_1}{dz} = \frac{\Delta\beta}{2}\hat{b}\_1 - \frac{\mathbb{C}\_o}{2}\hat{b}\_0 \exp\left(-i\gamma z\right). \tag{26}$$

*Appl. Sci.* **2020**, *10*, 8850

It is interesting to note that *a*ˆ0,1 and ˆ *b*0,1 are, after the spatial RWA, the same operators, except a global phase, that is *a*ˆ0,1 = *b*0,1 exp(*iβ*¯*z*), *β*¯ = (*β*<sup>0</sup> + *β*1)/2. Finally, as mentioned above, we can write the above equation in a matrix form acting on the single-photon state |*Lb*(*z*) = *b*0(*z*)|10 + *b*1(*z*)|11, where <sup>|</sup>*L*(*z*) <sup>=</sup> <sup>|</sup>*Lb*(*z*) exp(*iβ*¯*z*), that is,

$$i\hbar \frac{\partial}{\partial z}|L\_b(z)\rangle = -\frac{\hbar}{2} \begin{pmatrix} \Delta \beta & \mathbb{C}\_o \exp\left(i\gamma z\right) \\ \mathbb{C}\_o \exp\left(-i\gamma z\right) & -\Delta \beta \end{pmatrix} |L\_b(z)\rangle = 0$$

$$\mathbf{E} = -\frac{\hbar}{2}\mathbf{B} \begin{pmatrix} \cos a & \sin a \exp\left(i\gamma z\right) \\ \sin a \exp\left(-i\gamma z\right) & -\cos a \end{pmatrix} \left| L\_b(z) \right\rangle = -\hat{M}(z) \left| L\_b(z) \right\rangle,\tag{27}$$

with cos *<sup>α</sup>* = <sup>Δ</sup>*β*/B, sin *<sup>α</sup>* = *Co*/B, and B = <sup>Δ</sup>*β*<sup>2</sup> + *<sup>C</sup>*<sup>2</sup> *<sup>o</sup>* . This equation is just the spatial equivalent of Equation (2). A Hamiltonian operator is replaced by a momentum operator. Obviously, the dynamic properties are identical under the formal changes: *ω* ↔ *γ*, *γpBo* ↔ B, *θ* ↔ *α*. Note that a rotating equivalent vector (*γpB*)*eq* <sup>≡</sup> **<sup>B</sup>** <sup>=</sup> <sup>B</sup>(sin *<sup>α</sup>* cos *<sup>γ</sup>z*, <sup>−</sup> sin *<sup>α</sup>* sin *<sup>γ</sup>z*, cos *<sup>α</sup>*) is obtained. In short, we have achieved a photonic simulator of aspin-magnetic field coupling. However, we must stress some limitations for this simulator. The momentum operator *M*ˆ defined by (27) and simulating the coupling spin-magnetic field only is valid under the spatial RWA, that is the values of Δ*β* and *γ* have to be high and not too different. Therefore, we will be able to simulate the interaction with magnetic fields with a large *z*-component and oscillating with a frequency *ω* of the same order as the term *γpBo* cos *θ*.

Finally, it is interesting to obtain a constant momentum operator by applying a unitary transformation (rotating reference system) to the quantum state |*Lb*(*z*), that is,

$$|l(z)\rangle = \hat{\mathcal{U}}(z,\gamma)|L\_{\flat}(z)\rangle = \exp\left(-i\gamma vz/2\right)|L\_{\flat}(z)\rangle. \tag{28}$$

Therefore, by using |*Lb*(*z*) = exp (*iγσzz*/2)|*l*(*z*) in Equation (27), we obtain, after a certain, but straightforward calculation, the following momentum operator for the state |*l*(*z*):

$$\left| -\hat{M}'|l(z)\right\rangle = \frac{\hbar}{2} \begin{pmatrix} -\Delta\beta + \gamma & -\mathbb{C}\_o \\ -\mathbb{C}\_o & \Delta\beta - \gamma \end{pmatrix} \left| l(z) \right\rangle = -\frac{\hbar}{2} \text{D}\_o \begin{pmatrix} \cos a' & \sin a' \\ \sin a' & -\cos a' \end{pmatrix} \left| l(z) \right\rangle,\tag{29}$$

where cos *<sup>α</sup>* = (Δ*<sup>β</sup>* <sup>−</sup> *<sup>γ</sup>*)/D*o*, sin *<sup>α</sup>* <sup>=</sup> *Co*/D*o*, and D*<sup>o</sup>* <sup>=</sup> (Δ*<sup>β</sup>* <sup>−</sup> *<sup>γ</sup>*)<sup>2</sup> <sup>+</sup> *<sup>C</sup>*<sup>2</sup> *<sup>o</sup>* . As shown later, these results are important for simulating both dynamical and topological properties. By comparison between Equations (5) and (29), we obtain the following simulation parameters:

$$(\omega - \gamma\_p B\_o \cos \theta) \leftrightarrow (\Delta \beta - \gamma), \quad \gamma\_p B\_o \sin \theta \leftrightarrow \mathbb{C}\_o. \tag{30}$$

Next, it is interesting to compute the eigenstates |*l*(*z*) and eigenvalues *β<sup>l</sup>* of the momentum operator *<sup>M</sup>*<sup>ˆ</sup> given by Equation (6), that is <sup>|</sup>*l*(*z*) <sup>=</sup> exp *<sup>i</sup> <sup>h</sup>*¯ *plz* <sup>|</sup>*l*(0) <sup>=</sup> exp *iβlz* |*l*(0). After a standard calculation, we have:

$$|l\_{+}(0)\rangle = \cos\frac{a'}{2}|1\_0\rangle + \sin\frac{a'}{2}|1\_1\rangle, \quad |l\_{-}(0)\rangle = \sin\frac{a'}{2}|1\_0\rangle - \cos\frac{a'}{2}|1\_1\rangle. \tag{31}$$

with eigenvalues *<sup>β</sup>*<sup>±</sup> <sup>=</sup> <sup>±</sup>D*o*/2, that is linear momentums *<sup>p</sup>*<sup>±</sup> <sup>=</sup> <sup>±</sup>*h*¯ <sup>2</sup> *Do* = ±*po*. Therefore, by taking into account Equation (28), the full evolved states |*Lb*(*z*) can be easily obtained. Note that we are simulating the quantum state Ψ(*t*) given by Equation (4), and thus, for example, quantum processing based on NMR could be simulated by this photonic device. For the sake of expositional convenience, we will return to this question in the next subsection.

#### *3.4. Photonic Simulation of Light-Matter Interaction*

For the case of light-matter interaction, we have a more direct photonic simulation by using a two-mode planar guide perturbed by an integrated optical grating. Thus, by defining Δ*β* = *β*<sup>0</sup> − *β*<sup>1</sup> and *β*¯ = (*β*<sup>0</sup> + *β*1)/2 and for the sake of simulation purposes, choosing an initial phase *α<sup>o</sup>* = *π*, the momentum operator defined in Equation (20) can be rewritten as follows:

$$\hat{M}|L(z)\rangle = \hbar\left[ (\bar{\beta} + \mathbb{C}\_{+})\mathbf{I} + (\frac{\Delta\beta}{2} - \mathbb{C}\_{-})\sigma\_{z} - \mathbb{C}\_{0}\cos(\gamma z)\sigma\_{x} \right]|L(z)\rangle,\tag{32}$$

where *C*<sup>±</sup> = (*C*00(*z*) ± *C*11(*z*))/2. As in the mechanical case, these terms could be zero if the optical modes *en*,*m*(*x*) and perturbation Δ(*x*) have a suitable parity, that is even or odd modes along with an odd perturbation. Next, we perform the following relevant change in the vector representation of the single-photon state <sup>|</sup>*L*(*z*) ≡ *l*0(*z*) exp(*iγz*/2), *l*1(*z*) exp(−*iγz*/2) <sup>t</sup> <sup>≡</sup> exp(*iγσzz*/2)|*l*(*z*), with t indicating the transpose. By inserting this state into the above equation, using *<sup>M</sup>*<sup>ˆ</sup> <sup>=</sup> <sup>−</sup>*ih*¯ *<sup>∂</sup>*/*∂z*, and neglecting terms that are rapidly oscillating, that is exp(±*iqγz*) with *q* = 1/2, 1, 3/2, the following momentum operator is obtained for the state <sup>|</sup>*l*(*z*) ≡ *l*0(*z*), *l*1(*z*) t :

$$
\langle \hat{M}' | l(z) \rangle \approx \left( \hbar \bar{\beta} \, \mathbf{I} + \frac{\hbar \delta\_{\mathbf{s}}}{2} \sigma\_z - \frac{\hbar \mathcal{C}\_{\mathcal{o}}}{2} \sigma\_x \right) | l(z) \rangle,\tag{33}
$$

with *δ<sup>s</sup>* = Δ*β* − *γ*. This momentum operator simulates the Hamiltonian given by Equation (13), where the following simulation parameters are obtained:

$$
\omega \leftrightarrow -\bar{\beta}, \quad \delta = (\omega\_o - \omega) \leftrightarrow \delta\_s = (\Delta \beta - \gamma), \quad \Omega \leftrightarrow \mathbb{C}\_o. \tag{34}
$$

with *δ* ↔ *δ<sup>s</sup>* the detuning simulation parameter and Ω ↔ *Co* the Rabi frequency simulation parameter. Hence, this photonic device simulates light-matter interaction under a spatial RWA. Obviously, all temporal dynamics obtained by light-matter interaction can be simulated by means of this photonic device, as for example Rabi oscillations, logic gates for one-qubit transformations, and so on. In order to make clear these possibilities, let us consider the synchronous (or resonant) case, that is Δ*β* = *γ*. The momentum operator is thus simplified, and the solutions can be easily obtained. In matrix form, the solution of the equation above is given by:

$$
\begin{pmatrix} l\_0(z) \\ l\_1(z) \end{pmatrix} = \begin{pmatrix} \cos\frac{\mathbb{C}\_2}{2}z & -i\sin\frac{\mathbb{C}\_2}{2}z \\ -i\sin\frac{\mathbb{C}\_2}{2}z & \cos\frac{\mathbb{C}\_2}{2}z \end{pmatrix} \begin{pmatrix} l\_0(0) \\ l\_1(0) \end{pmatrix} \equiv \chi(\Theta/2) \begin{pmatrix} l\_0(0) \\ l\_1(0) \end{pmatrix},\tag{35}
$$

with <sup>Θ</sup>=*Coz* and where an irrelevant global phase *<sup>e</sup>iβ*¯*<sup>z</sup>* has been omitted. As a simple example, if we choose a length of the grating (interaction length) *z* = 2*π*/*Co*, then a Xquantum logic gate is implemented. We must stress that these are transformations corresponding to the so-called Θ-pulses in atom-light temporal interaction for computing purposes [22]. For example, given an input state |10, the state propagating along the optical grating will be:

$$|L(z)\rangle = l\_0(z)|1\_0\rangle + l\_1(z)|1\_1\rangle. \tag{36}$$

On the other hand, for the case <sup>|</sup>*δs*<sup>|</sup> *Co* and by omitting again the global phase *<sup>e</sup>iβ*¯*z*, we obtain, as a solution to Equation (33), a phase gate, that is,

$$
\begin{pmatrix} l\_0(z) \\ l\_1(z) \end{pmatrix} = \begin{pmatrix} \exp(i\frac{\delta\_s}{2}z) & 0 \\ 0 & \exp(-i\frac{\delta\_s}{2}z) \end{pmatrix} \begin{pmatrix} l\_0(0) \\ l\_1(0) \end{pmatrix} \equiv \mathcal{Z}(\Phi/2), \tag{37}
$$

where Φ = *δsz*. Therefore, a Z quantum logic gate is obtained from Equation (37) if *δsz*/2 = *π*/2, an S-gate if *δsz*/2 = *π*/4, a T-gate if *δsz*/2 = *π*/8, and so on. Moreover, by using an optical grating

implementing a transformation X(Θ/2) and two phase gates Z(Φ/2), a ZXZ-factorization of SU(2) is obtained; therefore, any unitary transformation can be implemented with the photonic device, and consequently, any one-qubit can be generated.

Next, we present solutions for the asynchronous case, that is Δ*β* = *γ*, which provides a most general solution and will be very useful to obtain geometric phases. By using standard methods to solve a linear equation system, the general solution of Equation (33) is given by:

$$
\begin{pmatrix} l\_0(z) \\ l\_1(z) \end{pmatrix} = \begin{pmatrix} \cos\frac{\delta\_r}{2}z + i\frac{\delta\_t}{\delta\_r}\sin\frac{\delta\_r}{2}z & -i\frac{\mathbb{C}\_x}{\delta\_r}\sin\frac{\delta\_r}{2}z \\ -i\frac{\mathbb{C}\_x}{\delta\_r}\sin\frac{\delta\_r}{2}z & \cos\frac{\delta\_t}{2}z - i\frac{\delta\_t}{\delta\_r}\sin\frac{\delta\_r}{2}z \end{pmatrix} \begin{pmatrix} l\_0(0) \\ l\_1(0) \end{pmatrix},\tag{38}
$$

with *<sup>δ</sup><sup>r</sup>* = *<sup>δ</sup>*<sup>2</sup> *<sup>s</sup>* + *C*<sup>2</sup> *<sup>o</sup>* and where, once more, the irrelevant global phase *<sup>e</sup>iβ*¯*<sup>z</sup>* is omitted. Note that for *δ<sup>s</sup>* = 0, the matrix given by Equation (35) is recovered. Likewise, under the condition |*δs*| *Co*, the matrix solution given by Equation (37) is obtained. We must recall that all these transformations are obtained in a spatial rotating reference system defined by Equation (28) and induced by the RWA approximation, as occurs in atom-light temporal interactions. Likewise, we must stress that by using the simulation parameters given by Equation (34), mechanical solutions for light-matter (atom) interaction are directly obtained. Finally, it is also easy to check that the same solutions are obtained for the spin-magnetic field interaction simulation given by Equation (29); therefore, such an interaction has the same quantum processing properties as the atom-optical field interaction.

Finally, it is worth paying attention to the problem of properly initializing a quantum state. For that, let us consider an SPDC (Spontaneous Parametric Down Conversion) source of biphotons <sup>|</sup>1*ka*1*k<sup>b</sup>* , that is twin photons excited in two spatial modes propagating along directions *<sup>k</sup><sup>a</sup>* and *<sup>k</sup>b*. Photon <sup>|</sup>1*k<sup>b</sup>* is directed towards an APD device, and the other one <sup>|</sup>1*k<sup>a</sup>* is directed to the prism-waveguide coupler, as shown in Figure 1a). The direction *k<sup>a</sup>* is chosen in such a way that, for example, if the fundamental mode of the planar guide is excited, that is *k<sup>a</sup>* = *k*0, then we obtain the single-photon state (or register) <sup>|</sup>10. Alternatively, direction *<sup>k</sup><sup>a</sup>* <sup>=</sup> *<sup>k</sup>*<sup>1</sup> can be chosen to excite the mode e1 of the planar guide, and thus, the single-photon state (or register) |11 is obtained. If we had channel waveguides, a similar procedure can be used, for example before the channel waveguide, there would be a planar waveguide with a prism-waveguide coupler and, next, an integrated lens or a similar integrated optical element [15] in such a way that the excited mode is focused onto the channel waveguide to excite the desired single-photon state. Next, by using optical gratings, different quantum transformations can be performed, and consequently, a general output state |*L* = *a*0(*z*)|1*k*<sup>0</sup> + *a*1(*z*)|1*k*<sup>1</sup> is obtained. Finally, at a distance *z*, another prism-waveguide coupler can be placed at the end of the device to detect the quantum state. Photon detections can be made by using coincidences between the output photon and the photon |1*k<sup>b</sup>* reaching the APD.

#### *3.5. Implementation of Photonic Simulators*

In this subsection, we present more complex simulators by using several integrated optical gratings along with other integrated components as Directional Couplers (DCs) made with Single-Mode (SMWs) or Two-Mode channel Waveguides (TMWs). We present a simulator of the interaction between an atom with four levels and Θ-pulses, next the interaction of a particle (or physical system) with spin 3/2 and a magnetic field, and finally, an arbitrary unitary transformation SU(4) implemented with optical gratings what allows reducing the number of paths by half or even to a quarter, which can be generalized to SU(N) transformations, which can be of interest for photonic simulators.

Let us consider, as the first example, the Optical Grating (OG) studied above, but with a number *d* = 4 of collinear guided modes whose propagation constants are *β<sup>j</sup>* with *j* = 1, 2, 3, 4. This simple device can simulate a four-level atom, which can be used to implement a CNOT logic gate operation under atom-laser interaction [22]. The optical grating fulfills *γ* = *β*<sup>2</sup> − *β*3; therefore, it will produce an efficient transition between Modes 2 and 3 if *z* = *π*/*Co* (*π*-pulse in atom-laser interaction) according to Equation (35). We can identify each single-photon state excited in each optical mode as a computational state (two-qubit single-photon [23]); thus, we have the state |*L* = *c*00|00 + *c*01|01 + *c*10|10 + *c*11|11. When this state goes through the optical grating, the output state is |*L* = *c*00|00 + *c*01|01 + *c*11|10 + *c*10|11, that is a CNOT operation is obtained. Obviously, this is not the best way to implement logic gates for two-qubits, but it makes clear how OGs can be used to simulate interaction between an atom with several levels and the so-called Θ-pulses.

The second example consists of using single-mode channel guides coupled by OGs in order to simulate the interaction between a particle (or physical system) with spin *s* and a magnetic field. This can be implemented by using a number *N*=2*s*+1 of optical modes excited in *N* channel waveguides and coupled by OGs, as shown in Figure 2 for the particular case of four guides (spin *s*=3/2). We must stress that coupling between parallel SMWs can be also described by the Heisenberg quantum equations given by Equation (18). In order to simulate spin-magnetic field interaction, either the strength of the OGs or the separation between consecutive SMWs has to be adjusted, and therefore, proper coupling strength values *Cj*,*j*+<sup>1</sup> with *j*=1, ...*N* − 1, between consecutive guides, are obtained. Likewise, the SMWs have to be designed with different propagation constants *β<sup>j</sup>* (asynchronous guides) by changing, for example, the depth of each channel waveguide. Finally, it is well known in integrated optics that for asynchronous DCs assisted by optical gratings, the coupling due to the overlapping fields is negligible. As a particular example, let us consider spin *s* = 3/2. The general equation system is given by:

$$-i\hbar \frac{d\mathfrak{d}\_{\vec{\jmath}}(z)}{dz} = \hbar \vec{\beta}\_{\vec{\jmath}}(z)\mathfrak{d}\_{\vec{\jmath}}(z) + \hbar \sum\_{j' \neq j} \mathbb{C}\_{j'\vec{\jmath}'}(z)\mathfrak{d}\_{j'}(z). \tag{39}$$

with (|*j* − *j* | < 2), that is coupling only exists between consecutive channel guides. Next, by using the same procedure followed for the case of spin-1/2, we obtain, after a long, but straightforward calculation, the following equation system:

$$-i\hbar \frac{\partial}{\partial z} |L\_b(z)\rangle = -\Delta\_o \mathbb{I}\_{4 \times 4} + \dots$$

$$\begin{pmatrix} 3\Lambda\_{\boldsymbol{\nu}} & \mathsf{C}\_{12}\exp\left(i\gamma z\right) & \mathbf{0} & \mathbf{0} \\ \mathsf{C}\_{21}\exp\left(-i\gamma z\right) & \Lambda\_{\boldsymbol{\nu}} & \mathsf{C}\_{23}\exp\left(i\gamma z\right) & \mathbf{0} \\ \mathbf{0} & \mathsf{C}\_{32}\exp\left(-i\gamma z\right) & -\Lambda\_{\boldsymbol{\nu}} & \mathsf{C}\_{34}\exp\left(i\gamma z\right) \\ \mathbf{0} & \mathbf{0} & \mathsf{C}\_{43}\exp\left(-i\gamma z\right) & 3\Lambda\_{\boldsymbol{\nu}} \end{pmatrix} \left| L\_{b}(z) \right\rangle \end{pmatrix} \begin{pmatrix} L\_{b}(z) \end{pmatrix} \tag{40}$$

where <sup>I</sup>4*x*<sup>4</sup> is the identity matrix and the following propagation constants are used: *<sup>β</sup>*1, *<sup>β</sup>*2=*β*1−Δ*o*, *β*3=*β*1−2Δ*o*, *β*4=*β*1−3Δ*o*. These values can be achieved by adjusting, for example, the width of the channel waveguides. On the other hand, the coupling coefficients for the gratings fulfil the relationships *C*12=*C*21=*C*34=*C*43= <sup>√</sup>3*Co* and *<sup>C</sup>*23=*C*32=2*Co*, which can be achieved by adjusting the separation between consecutive guides with optical modes ej(x, y), *j*=1, 2, 3, 4. In our case, Guides 2−3 are closer than 1−2 and 3−4 because the coupling between Guides 2*and*3 is larger, as shown in Figure 2. Finally, by defining cos *θ* = Δ*o*/Γ and sin *θ* = *Co*/Γ, with Γ = (Δ<sup>2</sup> *<sup>o</sup>* + *C*<sup>2</sup> *<sup>o</sup>* )1/2, and a fictitious magnetic field *<sup>B</sup><sup>f</sup>* = (sin *<sup>θ</sup>* cos *<sup>γ</sup>z*, <sup>−</sup> sin *<sup>θ</sup>* sin *<sup>γ</sup>z*, cos *<sup>θ</sup>*), we can write the above equation system as follows:

$$-i\hbar \frac{\partial}{\partial z} |L\_b(z)\rangle = \left\{-\hbar \Delta\_o \mathbb{I}\_{4\ge 4} + \frac{\hbar}{2} \Gamma \left| J(3/2) \right.\\ \left. \mathcal{B}\_f \right\} |L\_b(z)\rangle \right. \tag{41}$$

where *J*(3/2)=(*Jx*(3/2), *Jy*(3/2), *Jz*(3/2)) are the spin matrices for spin *s* = 3/2. In short, we constructed a photonic simulator for interaction between a spin-3/2 particle, or physical system, and a periodic magnetic field.

**Figure 2.** Four non-synchronous Single-Mode channel Waveguides (SMWs) assisted by an optical grating to simulate the interaction between a spin *s* = 3/2 particle and an oscillating magnetic field.

Next, we present a general unitary transformation SU(4) by using optical gratings. With this example, we want to show that OGs allow reducing the number of paths required for constructing a simulator. Indeed, an SU(4) simulator is formed by four paths implemented by SMWs with optical modes e(j) <sup>00</sup> (x, y) (*j*=1, ...,4), six DCs, and six Phase Shifters (PSs), as shown in Figure 3 (bottom). However, if OGs are used, we only need two paths, that is TMWs with modes e (j) <sup>00</sup> (x, y), e(j) <sup>10</sup> (x, y) (*j* = 1, 2) and six OGs, as shown in Figure 3 (up). In this case, we also use Selective Directional Couplers (SDCs) (a SDC always performs an X transformation in one mode, and the other mode undergoes an identity transformation). Note that the number of paths was reduced by half because OGs act on collinear modes, unlike directional couplers, which act on codirectional modes. Obviously, if we had used OGs with four modes, we could reduce the number of paths by 1/4. Overall, we will be able to reduce the number of paths by 1/*d* if we use OGs with a number *d* of collinear modes. Ultimately, we can take advantage of the OGs to increase integration and thus to implement more flexible and scalable photonic simulators such as for example boson sampling ones [24], where the number of paths would be reduced by half. In short, the number *N* of paths of any required SU(*N*) transformation [25] could be reduced up to *N*/*d*.

**Figure 3.** Standard implementation (bottom) of an arbitrary unitary transformation SU(4) by using Single-Mode channel Waveguides (SMWs) and Directional Couplers (DCs) (**bottom**). Alternative implementation (**top**) by using Two-Mode channel Waveguides (TMWs), Selective Directional Couplers (SDCs) and Optical Gratings (OGs). PS, Phase Shifter.

Finally, it is important to indicate that quantum photonic devices have their own limitations [4]. Thus, the difficulty to implement two-qubit logic gates is well known, which is, at present, an important drawback in general purpose quantum photonic computation. However, quantum photonic simulation, or simply quantum photonic computation, for specific purposes can provide efficient technological

solutions, particularly if new degrees of freedom are incorporated. In our case, we have just shown that integrated OGs allow processing with several collinear modes, which improves the optical integration for high-dimensional problems, that is it provides a moderate increase of the on-chip flexibility and scalability for photon-based quantum simulation. Moreover, OGs enable simulating quantum devices under variable perturbations and, in particular, periodic perturbations, which usually appear in quantum systems interacting with fields like spin (*N*=2*s*+1 modes)-magnetic field interaction, atom (*d* modes)-optical field interaction, and so on.

#### **4. Quantum Geometric Phases**

So far, we have simulated the dynamics of the spin-magnetic and light-matter interaction systems with a photonic device. However, these quantum devices also generate topological or geometric phases besides the dynamic phases. Geometric phases are precisely due to geometric properties as was originally proven by Berry in his seminal work about an adiabatic quantum system [26]. Later, these geometric phases were generalized by Aharonov and Anandan [27] to non-adiabatic processes, and their calculation is made by using the geometric properties of the projective Hilbert space. Finally, a useful extension to geometric phases associated with non-cyclic circuits on the projective Hilbert space was also proposed [28]. On the other hand, topological phases in optics have also been extensively studied in bulk devices with polarization modes [29] and also in integrated optics with spatial modes [30]. At present, geometrical phases have regained interest for their possible application to geometric quantum computation [31,32]. Non-adiabatic spatial propagation on the Hilbert space generate the geometric phase known as the Aharonov–Anandan (AA) phase. It is well known that the spin-magnetic field interaction, as for example NMR, produces AA phases. We must stress that the geometric phase for a two-dimensional projective Hilbert space can be calculated in a geometric way as *φ*<sup>g</sup> = (1/2)**Ω**(C), where **Ω**(C) is the solid angle subtended by the circuit C followed by the quantum state on the Bloch sphere. In this section, we prove that quantum geometrical phases can be obtained by an integrated photonic grating, and therefore, it simulates the geometric phases produced by both a spin-magnetic field system and an atom-optical field system.

#### *4.1. Geometric Phases in Spin-Magnetic Field Photonic Simulation*

Let us consider the eigenstates |*l*±(0) given by Equation (31). Spatial propagation of these states is given, according to Equations (28) and (31), by |*Lb*(*z*)≡|*L*±(*z*) = exp(*iγσzz*/2)|*l*±(*z*) = exp(*ip*<sup>±</sup> *<sup>h</sup>*¯ *z*) exp(*iγσzz*/2)|*l*±(0). For instance, let us take |*L*+(*z*) after a propagation distance *z* = *ν*Λ, that is for a photonic integrated grating with a length *ν*Λ, where Λ = 2*π*/*γ* and *ν* is the number of cycles taken, then we have:

$$\langle |L\_{+}(\nu\Lambda)\rangle = \exp\left(\frac{ip\_{+}}{\hbar}\nu\Lambda\right) \left(\cos\frac{a'}{2}\exp(i\nu\pi)|1\_{0}\rangle + \sin\frac{a'}{2}\exp(-i\nu\pi)|1\_{1}\rangle\right) \tag{42}$$

By following the same procedure for the state |*L*−(*z* = *ν*Λ) and taking into account that *ν* is an integer, we obtain:

$$<|L\_{\pm}(\nu\Lambda)\rangle = \exp\left(\frac{ip\_{\pm}}{\hbar}\nu\Lambda \pm i\nu\pi\right)|l\_{\pm}(0)\rangle = \exp(i\phi^{\pm})|l\_{\pm}(0)\rangle.\tag{43}$$

The phases *φ*<sup>±</sup> obtained above are the full phases acquired by the states |*L*± after *ν* cycles in the optical grating, that is,

$$
\phi^{\pm} = \nu(\pm p\_o \Lambda/\hbar \pm \pi) \tag{44}
$$

A photonic Bloch sphere is shown on the left in Figure 4 where each point corresponds to a single-photon state given by Equation (19). For comparison purposes, an NMR Bloch sphere is also shown on the right. A single-photon state propagating along *z* can be represented by the following general expression |*L*(*z*) = *c*0(*z*)|10 + *c*1(*z*)|11; therefore, each point (*x*, *y*, *z*) of the photonic Bloch sphere is defined as follows: *x* = 2Re*c*0(*z*)*c*1(*z*), *y* = 2 Im*c*0(*z*)*c*1(*z*), *z* = |*c*0(*z*)| <sup>2</sup> − |*c*1(*z*)<sup>|</sup> 2, where Re and Im stand for real and imaginary parts, respectively.

It is easy to check that eigenstates |*L*±(*z*) follow the curves *Cc* and *C <sup>c</sup>* corresponding to spherical caps, as shown in Figure 4. The solid angle subtended by a spherical cap with an angular extension *α* is given by **Ω**(*C*) = 2*π*(1 − cos *α* ).

As the full phase *φ* can be decomposed into a dynamical part *φ*<sup>d</sup> and a geometric one *φ*g, then the geometric phase can be obtained from the relationship:

$$
\Phi\_{\rm g}^{\pm} = \phi\_{\pm} - \phi\_{\rm d}^{\pm} \tag{45}
$$

We first compute the dynamical phase, that is *φ*± <sup>d</sup> <sup>=</sup> <sup>1</sup> *h*¯ 3 *<sup>ν</sup>*<sup>Λ</sup> <sup>0</sup> *L*±(*z*)|*M*<sup>ˆ</sup> (*z*)|*L*±(*z*)dz; therefore, by taking into account the transformation (28) and expression (29) for *M*ˆ , we obtain:

$$\begin{split} \boldsymbol{\Phi}\_{\text{d}}^{\pm} &= \frac{1}{\hbar} \int\_{0}^{\eta \Lambda} \langle L\_{\pm}(z) | \hat{\mathcal{M}}(z) | L\_{\pm}(z) \rangle \mathrm{d}z = \frac{1}{\hbar} \int\_{0}^{\eta \Lambda} \langle l\_{\pm}(0) | \hat{\mathcal{M}}(0) | l\_{\pm}(0) \rangle \mathrm{d}z = \\ &= \frac{1}{\hbar} \int\_{0}^{\eta \Lambda} \left( \langle l\_{\pm}(0) | \hat{\mathcal{M}}' | l\_{\pm}(0) \rangle - \frac{\gamma}{2} \langle l\_{\pm}(0) | \boldsymbol{\sigma}\_{\mp} | l\_{\pm}(0) \rangle \right) \mathrm{d}z = \\ &= \frac{p\_{\pm}}{\hbar} \nu \Lambda \pm \nu \pi \cos \boldsymbol{\mathfrak{a}}', \end{split} \tag{46}$$

and by taking into account the total phase given by Equation (43), the geometrical phases acquired by the eigenstates are given by:

$$\phi\_{\mathbf{g}}^{\pm} = \pm \nu \pi (1 - \cos a') = \pm \nu \pi \left( 1 - \frac{\Delta \beta - \gamma}{\sqrt{(\Delta \beta - \gamma)^2 + \mathbf{C}\_o^2}} \right). \tag{47}$$

Note that the geometric phase is half the solid angle subtended by the circuit, where the sign ± depends on the direction of rotation followed by the state on the Bloch sphere. The corresponding spherical cap circuits C*<sup>c</sup>* and C *<sup>c</sup>* are shown in Figure 4. Moreover, the geometric phase depends on the device parameters, that is Δ*β*, *γ*, *Co*. Likewise, it is important to indicate that the dynamical phase is of order *ω*<sup>1</sup> (note that *β* ∝ *ω*), but the geometrical phase is of order *ω*0; therefore, the geometric phase is less sensitive to errors in the distance propagation.

**Figure 4.** On the left, a photonic Bloch sphere shows the evolution of the single-photon states |*l*+ and |*l*− (spherical cap circuits C*<sup>c</sup>* and C *<sup>c</sup>*). Likewise, a spherical wedge circuit C*<sup>w</sup>* is shown. On the right, the simulated mechanical Bloch sphere for an interacting atom-field system is shown, with similar spherical cup circuits for states |*η*+ and |*η*− and also a spherical wedge circuit C*w*.

In short, a single-photon state acquires an AA geometric phase under propagation in an integrated photonic grating. The same expression is found for a spin-1/2 particle in a magnetic field; therefore, topological simulations can be made. Thus, by applying the simulation parameters given by Equation (30), geometric phases can be obtained.

#### *4.2. Elimination of the Dynamical Phase in Spin-Magnetic Field Photonic Simulation*

It is well known in the mechanical case that the geometric phase is hidden in the spin-magnetic field interaction because it is combined with the dynamical phase within *φ*±; therefore, the dynamical phase has to be eliminated in order to take advantage of the properties of a geometric phase. We present a photonic solution, which is similar to the one used in the mechanical case, that is if the quantum state, after evolution under a first Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> <sup>1</sup> <sup>=</sup> *<sup>H</sup>*<sup>ˆ</sup> for time *<sup>T</sup>*, finds a second Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> <sup>2</sup> <sup>=</sup> <sup>−</sup>*H*<sup>ˆ</sup> , then the dynamical phase are mutually canceled; however, the eigenstates do not change, therefore neither does the geometrical phase.

In the photonic case, we have to find a new momentum operator, that is a new integrated optical grating, such as *<sup>M</sup>*<sup>ˆ</sup> <sup>2</sup> <sup>=</sup> <sup>−</sup>*M*<sup>ˆ</sup> . We assume that such a new optical grating has the same frequency *γ*; therefore, we have the same rotating system, that is the same transformation (28). Accordingly, the condition for eliminating the dynamical phase is obtained from the operator *M*ˆ , that is,

$$
\hat{M}'\_2 = \frac{\hbar}{2} \begin{pmatrix} \Delta \p' - \gamma & \mathbf{C}'\_o \\ \mathbf{C}'\_o & -\Delta \p' + \gamma \end{pmatrix} = -\hat{M}' = -\frac{\hbar}{2} \begin{pmatrix} \Delta \p - \gamma & \mathbf{C}\_o \\ \mathbf{C}\_o & -\Delta \p + \gamma \end{pmatrix};\tag{48}
$$

therefore, *C <sup>o</sup>* = −*Co* and (Δ*β* − *γ*) = −(Δ*β* − *γ*). The first condition can be achieved by introducing an initial phase in the second grating, that is Δ(*x*, *z*) = Δ(*x*) cos(*γz* + *π*), and the second one is achieved if Δ*β* = 2*γ* − Δ*β*. These results indicate that we need an additional grating with a new difference between propagation constants (linear momentum of the photon) Δ*β* = (*β* <sup>0</sup> − *β* 1). In Figure 5 the system for eliminating the dynamical phase is shown . Therefore, the total phases are *φ*<sup>±</sup> = *ν*(∓*p*Λ/¯*h* − *π*), and the dynamical phases are:

$$\phi\_{\rm d}^{\pm} = \frac{p\_{\pm}}{\hbar} \nu \Lambda \mp \nu \pi \cos a' \tag{49}$$

Therefore, the total dynamical phase is Φ*<sup>d</sup>* = 0, and the total geometrical phase after the single-photon state propagates through the two integrated gratings is twice the value acquired in the first grating, that is,

$$
\Phi\_{\vec{g}}^{\pm} = \mp 2\nu \pi (1 - \cos a') = \mp \Phi\_{\vec{\chi}}.\tag{50}
$$

Alternatively, propagation constants can be unchanged, and the grating frequency can be modified, that is *γ* = 2Δ*β* − *γ*; however, in this case, the transformation (28) must be applied with the factor *γ* . In short, we eliminated the dynamical phase; therefore, these results could be used for implementing logic gates or transformations based on topological phases, which are much more insensitive to fabrication errors, unlike dynamical phases, which as mentioned are of order *ω*. With these results, robust P-gates (Phase gates) can be designed; thus, the following transformation is implemented between the eigenstates:

$$P = \begin{pmatrix} \exp(-i2\nu\pi\cos a') & 0\\ 0 & \exp(i2\nu\pi\cos a') \end{pmatrix} = \exp(-i\Phi\_{\mathcal{S}}\sigma\_z) \tag{51}$$

Note that for 4*ν* cos *α* = ±1, a Z-gate is obtained, and for 4*ν* cos *α* = ±1/2, an S-gate is implemented, and so on.

**Figure 5.** Elimination of the dynamical phase by using two consecutive integrated optical gratings. The second grating has an initial phase *π*. Likewise, a prism is used to make a projective measure of states |10 and |11 for obtaining the probabilities *P*<sup>0</sup> and *P*<sup>1</sup> and, therefore, the geometric phase Φ*g*.

#### *4.3. Geometric Phases with Other Quantum States*

On the other hand, the eigenstates given by Equation (31) can be written as single-photon states excited in rotated optical modes, that is *e*+(*x*, *y*) = cos *<sup>α</sup>* <sup>2</sup> *eo*(*x*, *<sup>y</sup>*) + sin *<sup>α</sup>* <sup>2</sup> *e*1(*x*, *y*) and *e*−(*x*, *y*) = <sup>−</sup> sin *<sup>α</sup>* <sup>2</sup> *eo*(*x*, *<sup>y</sup>*) + cos *<sup>α</sup>* <sup>2</sup> *e*1(*x*, *y*); therefore:

$$|l\_{+}(0)\rangle = \cos\frac{a'}{2}|1\_0\rangle + \sin\frac{a'}{2}|1\_1\rangle = |1\_+\rangle,\quad |l\_{-}(0)\rangle = -\sin\frac{a'}{2}|1\_0\rangle + \cos\frac{a'}{2}|1\_1\rangle = |1\_-\rangle\tag{52}$$

The elimination of the dynamical phase means that these eigenstates have undergone the transformation *<sup>e</sup>*±*i*Φ*<sup>g</sup>* <sup>|</sup>1±; therefore, the following formal relationships can be written:

$$e^{i\Phi\_{\mathcal{S}}}|1\_{+}\rangle = e^{i\Phi\_{\mathcal{S}}}\hat{a}\_{+}^{\dagger}|0\_{\pm}\rangle, \quad e^{-i\Phi\_{\mathcal{S}}}|1\_{-}\rangle = e^{-i\Phi\_{\mathcal{S}}}\hat{a}\_{-}^{\dagger}|0\_{\pm}\rangle \tag{53}$$

with Φ*<sup>g</sup>* = −2*νπ*(1 − cos *α* ). Accordingly, the following transformations, induced by geometric phases, for the absorption operators are obtained:

$$\hbar\_{\pm}(z=n\Lambda) = e^{\pm i\Phi\_{\mathbb{S}}}\hbar\_{\pm}(0) \tag{54}$$

This result can be used for obtaining geometric phases of other quantum light states. As an example, we present two states, that is the number photon state or Fock state |*n*+ and the coherent state |*α*+, where subindex + indicates that the state is excited in the optical mode *e*+(*x*, *y*). The Fock state under propagation becomes:

$$|n\_{+}(0)\rangle = \frac{1}{\sqrt{n\_{+}!}} (\mathfrak{d}\_{+}(0)^{\dagger})^{n\_{+}}|0\rangle \longrightarrow \frac{1}{\sqrt{n\_{+}!}} e^{i n\_{+} \Phi\_{\mathfrak{F}}} (\mathfrak{d}\_{+}^{\dagger})^{n\_{+}}|0\rangle = e^{i n\_{+} \Phi\_{\mathfrak{F}}} |n\_{+}\rangle \tag{55}$$

Therefore, the quantum state has acquired a geometric phase *n*+Φ*g*. Likewise, the coherent state can be rewritten by using the complex displacement operator, that is,

$$|a\_{+}(0)\rangle = e^{(a\_{+}\hbar\_{+}^{\dagger}(0) - c.h.)}|0\rangle \longrightarrow e^{(a\_{+}\varepsilon^{j\Phi\_{\xi}}\hbar\_{+}^{\dagger} - c.h.)}|0\rangle = |\varepsilon^{j\Phi\_{\xi}}a\_{+}\rangle\tag{56}$$

Therefore, the geometric phase is Φ*g*, that is the same as the one acquired by a single photon. The same procedure can be applied to any other quantum light state excited in the integrated photonic grating.

#### *4.4. Optical Measurement of Geometric Phases*

Finally, we present how to measure the geometric phase starting from the measurements of the single-photon detection probability, which can be extracted by a prism-waveguide coupler [10], as shown in Figure 5. We focus on the spin-magnetic field interaction simulation case. By assuming that the input state is a single photon excited in the mode *e*0(*x*, *y*) and taking into account the relationships given by Equation (52), the following final state is obtained after the two gratings:

$$|L(z)\rangle = (\cos^2\frac{a'}{2} + \sin^2\frac{a'}{2}e^{i\Phi\_\mathcal{S}})|10\rangle + \cos\frac{a'}{2}\sin\frac{a'}{2}(1 - e^{i\Phi\_\mathcal{S}})|1\_1\rangle \tag{57}$$

The probability of the detection of a photon in Mode 0, that is *P*0, or in Mode 1, *P*1, is a function of the geometric phase, that is,

$$P\_1 = \frac{\sin^2 \alpha}{2} \left( 1 - \cos \Phi\_{\mathcal{g}} \right), \quad P\_0 = 1 - P\_1 \tag{58}$$

If Φ*<sup>g</sup>* = 2*π*, then *P*<sup>1</sup> = 0 and *P*<sup>0</sup> = 1, and if Φ*<sup>g</sup>* = *π*, then *P*<sup>1</sup> = sin2 *α* and *P*<sup>0</sup> = cos2 *α*. Therefore, from the measurement of *P*<sup>1</sup> and *P*0, the phase Φ*<sup>g</sup>* is obtained. In Figure 5 it is shown the projective measure of states |10 and |11 by a prism-waveguide coupler, which projects these state in different spatial directions. Finally, note that by using a coherent state, these probabilities are proportional to the intensity of the light, what can be called a semiclassical optical characterization, or in the most technical way, the geometric phase is also acquired by the classical fields, but would rigorously correspond to the so-called Hannay phase [33].

#### *4.5. Geometric Phases in Light-Matter Photonic Simulation*

Finally, we check that light-matter simulation can be also used to obtain geometric phases. For the sake of simplicity, we show geometric phases for wedge circuits, although more general cases can be studied. Let us consider an initial state |*L*(0) = |10, that is the point (0, 0, 1) on the photonic Bloch sphere. Next, let us consider an asynchronous optical grating with *δ<sup>s</sup> Co*; therefore, according to Equation (37) (phase gate), for *δszo* = 3*π*/2, we obtain |*L*(*zo*) = (1/ <sup>√</sup>2)(|10 <sup>+</sup> *<sup>i</sup>*|11), that is the state reaches the Bloch sphere point (0, 1, 0). Next, we consider that the grating has a greater coupling coefficient *Co*, then the single-photon state is given by the expression:

$$|L(z)\rangle = \frac{1}{\sqrt{2}} [ (\cos \frac{\delta\_r z}{2} + b \sin \frac{\delta\_r z}{2} + ia \sin \frac{\delta\_r z}{2}) |1\_0\rangle + (a \sin \frac{\delta\_r z}{2} + i \cos \frac{\delta\_r z}{2} - ib \sin \frac{\delta\_r z}{2}) |1\_1\rangle ] \tag{59}$$

where *a* = *δs*/*δ<sup>r</sup>* and *b* = *Co*/*δr*. If we choose a distance *z* = *z*<sup>1</sup> such as *δrz*1/2 = *π*/2, then, after a certain calculation, the following state is obtained:

$$|L(z\_1)\rangle = \frac{e^{i\phi\_\sigma}}{\sqrt{2}}(|1\_0\rangle + i|1\_1\rangle)\tag{60}$$

where *φ<sup>o</sup>* = atn(*a*/*b*) = atn(*δs*/*Co*), that is *a* = sin *φ<sup>o</sup>* and *b* = cos *φo*. The state has reached the point (0, −1, 0) of the photonic Bloch sphere. Now, we show that *φ<sup>o</sup>* is a geometric phase. Indeed, for the sake of symmetry, the state before reaching the above state has crossed the meridian *y* = 0 when *δrz*/2 = *π*/4, that is, the state:

$$|L'\rangle = \frac{(1+b+ia)}{2}|1\_0\rangle + \frac{(a+i(1-b))}{2}|1\_1\rangle \equiv m\_0 \epsilon^{ic\_0} |1\_0\rangle + m\_1 \epsilon^{ic\_1} |1\_1\rangle. \tag{61}$$

It is easy to prove that both states have the same phase, that is <sup>0</sup> = 1, and the modulus is given by *m*<sup>0</sup> = sin(*φo*/2) and *m*<sup>1</sup> = cos(*φo*/2); therefore, the state crosses the point *P* = (sin *φo*, 0, cos *φo*) of the photonic Bloch sphere as indicated in Figure 4 for *ϕ* = *φo*. Now, we must recall that partial

cycles also generate geometric phases, which can be calculated by closing the end points of the open cycle by a geodesic line [28,30]. In our case, the corresponding geodesic lies along the meridian at the plane *x* = 0, from point (0, −1, 0) to initial point (0, 0, 1). In short, the state has followed a wedge circuit C*w*, as shown in Figure 4. Now, we calculate the subtended solid angle by this wedge circuit. It is easy to check that a wedge circuit with angle *φ<sup>o</sup>* subtends a solid angle **Ω**(C*w*) = 2*φo*; therefore, the geometric phase is Φ*<sup>g</sup>* = (1/2)**Ω**(C*w*) = *φ<sup>o</sup>* = atn(*δs*/*Co*), which is just the global phase obtained in Equation (60). It is worth underlining that in this case, the geometric phase is not hidden, and therefore, the dynamical phase does not have to be eliminated. Moreover, this geometric phase can be also obtained in a atom-optical field system out of resonance by using Θ-pulses, with the first optical field with low amplitude *Eo* (Z gate) and the next with a higher amplitude *Eo*. By using the simulation parameters given by Equation (34), the geometric phase would be Φ*<sup>g</sup>* = atn(*δ*/Ω).

#### **5. Conclusions**

We propose a quantum photonic device based on integrated optical gratings in a two-mode slab guide to simulate the interaction between external fields and atoms. By using single-photon states, we study the simulations of a spin-(1/2)-magnetic field system, as for example nuclear magnetic resonance, and a two-level atom-optical field system corresponding to light-matter simulation. Both dynamical and geometric properties are simulated, in particular the geometric phases obtained by the mentioned systems. We prove that dynamical properties can be simulated for a wide range of cases with practical interest, although in the spin-(1/2)-magnetic field system it is restricted to relatively high values of frequency and magnetic field amplitude *Bo*. Overall, atom-optical field interaction does not present these restrictions. This study of integrated optical gratings opens up possibilities to more general simulations if several modes are used. Thus, spin (*s*)-magnetic field interaction simulations could be implemented by using a number *N* = 2*s*+1 of codirectional optical modes assisted by optical gratings; multilevel atom (*n*)-optical field interaction can be simulated by using *N* = *n* collinear optical modes coupled by optical gratings, which can in turn simulate, for example, two-qubit single photon logic gates, which has a high interest in quantum information systems; likewise, optical gratings allow interaction between *d* collinear modes, and thus, simulators based on *N* codirectional modes can reduce the number of paths used up to *N*/*d*, which improves the optical integration of the photonic simulator. On the other hand, AA geometric phases have been also obtained for both systems. The spin-(1/2)-magnetic field system requires dynamic phase cancellation, which is simulated by using two optical gratings; however, in the atom-optical field system, such cancellation is not required. Obviously, we must emphasize that although the proposed integrated photonic device is intended for quantum simulation, it can also be used to implement quantum operations and/or effects with the photonic device itself, such as logic gates, geometric phases, and so on, by using single-photon states or more general quantum states, as shown.

**Author Contributions:** All authors contributed equally to this work. All authors read and agreed to the published version of the manuscript.

**Funding:** Xunta de Galicia, Consellería de Educación, Universidades e FP, Grant GRC Number ED431C2018/11; Ministerio de Economía, Industria y Competitividad, Gobierno de España, Grant Number AYA2016-78773-C2-2-P.

**Acknowledgments:** One of the authors (G.M.C.) wishes to acknowledge the financial support by Xunta de Galicia, Consellería de Educación, Universidades e FP, by a predoctoral grant co-financed with the European Social Fund.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.

*Appl. Sci.* **2020**, *10*, 8850

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Towards Quantum 3D Imaging Devices**

**Cristoforo Abbattista 1, Leonardo Amoruso 1, Samuel Burri 2, Edoardo Charbon 2, Francesco Di Lena 3, Augusto Garuccio 3,4, Davide Giannella 3,4, Zden ˇek Hradil 5, Michele Iacobellis 1, Gianlorenzo Massaro 3,4, Paul Mos 2, Libor Motka 5, Martin Paúr 5, Francesco V. Pepe 3,4,\*, Michal Peterek 5, Isabella Petrelli 1, Jaroslav Reháˇ ˇ cek 5, Francesca Santoro 1, Francesco Scattarella 3,4, Arin Ulku 2, Sergii Vasiukov 3, Michael Wayne 2, Claudio Bruschini 2,†, Milena D'Angelo 3,4,\*,†, Maria Ieronymaki 1,† and Bohumil Stoklasa 5,†**


**Abstract:** We review the advancement of the research toward the design and implementation of quantum plenoptic cameras, radically novel 3D imaging devices that exploit both momentum–position entanglement and photon–number correlations to provide the typical refocusing and ultra-fast, scanning-free, 3D imaging capability of plenoptic devices, along with dramatically enhanced performances, unattainable in standard plenoptic cameras: diffraction-limited resolution, large depth of focus, and ultra-low noise. To further increase the volumetric resolution beyond the Rayleigh diffraction limit, and achieve the quantum limit, we are also developing dedicated protocols based on quantum Fisher information. However, for the quantum advantages of the proposed devices to be effective and appealing to end-users, two main challenges need to be tackled. First, due to the large number of frames required for correlation measurements to provide an acceptable signal-to-noise ratio, quantum plenoptic imaging (QPI) would require, if implemented with commercially available high-resolution cameras, acquisition times ranging from tens of seconds to a few minutes. Second, the elaboration of this large amount of data, in order to retrieve 3D images or refocusing 2D images, requires high-performance and time-consuming computation. To address these challenges, we are developing high-resolution single-photon avalanche photodiode (SPAD) arrays and high-performance low-level programming of ultra-fast electronics, combined with compressive sensing and quantum tomography algorithms, with the aim to reduce both the acquisition and the elaboration time by two orders of magnitude. Routes toward exploitation of the QPI devices will also be discussed.

**Keywords:** quantum imaging; plenoptic imaging; quantum correlations; SPAD arrays; quantum fisher information; compressive sensing

#### **1. Introduction**

Fast, high-resolution, and low-noise 3D imaging is highly required in the most diverse fields, from space imaging to biomedical microscopy, security, industrial inspection, and cultural heritage [1–5]. In this context, conventional plenoptic imaging represents one of the most promising techniques in the field of 3D imaging, due to its superb temporal resolution:

**Citation:** Abbattista, C.; Amoruso, L.; Burri, S.; Charbon, E.; Di Lena, F.; Garuccio, A.; Giannella, D.; Hradil, Z.; Iacobellis, M.; Massaro, G.; et al. Towards Quantum 3D Imaging Devices. *Appl. Sci.* **2021**, *11*, 6414. https://doi.org/10.3390/app11146414

Academic Editors: Maria Bondani, Alessia Allevi and Stefano Olivares

Received: 5 May 2021 Accepted: 7 July 2021 Published: 12 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

3D imaging is realized in a single shot, for seven frames per second at 30 M pixel resolution, and 180 frames per second for 1M pixel resolution [5]; no multiple sensors, near-field techniques, time-consuming scanning or interferometric techniques are required. However, conventional plenoptic imaging entails a loss of resolution which is often unacceptable. Our strategy to break such limitation consists of combining a radically new and foundational approach with last-generation hardware and software solutions. The fundamental idea is to exploit the information stored in correlations of light by using novel sensors and measurement protocols to implement a very ambitious task: high speed (10–100 fps) quantum plenoptic imaging (QPI) characterized by ultra-low noise and an unprecedented combination of resolution and depth-of-field. The developed imaging technique aims at becoming the first practically usable and properly "quantum" imaging technique that surpasses the intrinsic limits of classical imaging modalities. In addition to the foundational interest, the quantum character of the technique allows for extracting information on 3D images from correlations of light at very low photon fluxes, thus reducing the scene exposure to illumination. The interest in QPI is specially motivated by the potential advantages with respect to other established 3D imaging techniques. Actually, other methods require, unlike QPI, either delicate interferometric measurements, as in digital holographic microscopy [6,7], or phase retrieval algorithms, as in Fourier ptychography [8], or fast pulsed illumination, as in time-of-flight (TOF) imaging [9–14]. Moreover, QPI provides the basis of a scanning-free microscopy modality, overcoming the drawbacks of the confocal methods [15].

In view of the deployment of quantum plenoptic cameras suitable for real-world applications, the crucial challenge is represented by the reduction of both the acquisition and data elaboration times. In fact, a typical complication arises in quantum imaging modalities based on correlation measurements: the reconstruction of the correlation function encoding the desired image requires collecting a large number of frames (30,000–50,000 in the first experimental demonstration of the refocusing capability of correlation plenoptic imaging [16]), which must be individually read and stored before elaborating the output. Therefore, to get an estimate of the total time required to form a quantum plenoptic image, the data reading and transmission times must be added to the acquisition time of the employed sensor. This problem is addressed by an interdisciplinary approach, involving the development of ultrafast single-photon sensor systems, based on SPAD arrays [17–22], the optimization of circuit electronics to collect and manage the high number of frames (e.g., by GPU) [23,24], the development of dedicated algorithms (compressive sensing, machine learning, quantum tomography) to achieve the desired SNR with a minimal number of acquisitions [25–28]. Finally, the performances of QPI will be further enhanced by a novel approach to imaging based on quantum Fisher information [29,30]. Treating the physical model of plenoptic imaging in the view of quantum information theory brings new possibilities of improving the setup towards super-resolution capability in the object 3D space. Having the optimal set of optical setup parameters enables object reconstruction close to ultimate limits set by nature. In this article, we will combine a review of the state of the art in the aforementioned fields with the discussion on how they contribute to the development of QPI and its application.

The paper is organized as follows: in Section 2, we discuss the working principle and recent advances of Correlation Plenoptic Imaging (CPI), a technique that represents the direct forerunner of QPI; in Section 3, we present the hardware innovations currently investigated to reduce the acquisition times in CPI; in Section 4, we review the algorithmic solutions to improve QPI; in Section 5, we outline the perspectives of our future work in the context of Qu3D project; in Section 6, we discuss the relevance of our research. The work presented in this paper involves experts from three scientific research institutions, Istituto Nazionale di Fisica Nucleare (INFN, Italy), Palacky University Olomouc/Department of Optics (UPOL, Czechia), and Ecole polytechnique fédérale de Lausanne/Advanced Quantum Architecture Lab (EPFL, Switzerland), and from the industrial partner Planetek Hellas E.P.E. (PKH, Greece). The activity is carried within the project "Quantum 3D Imaging at high speed and high resolution" (Qu3D), founded by the 2019 QuantERA call [31].

#### **2. Plenoptic Imaging with Correlations: From Working Principle to Recent Advances**

Quantum plenoptic cameras promise to offer the advantages of plenoptic imaging, primarily ultrafast and scanning-free 3D imaging and refocusing capability, with performances that are beyond reach for the classical counterpart. State-of-the-art plenoptic imaging devices are able to acquire multi-perspective images in a single shot [5]. Their working principle is based on the simultaneous measurement of both the spatial distribution and the propagation direction of light in a given scene. The acquired directional information translates into refocusing capability, augmentable depth of field (DOF), and parallel acquisition of multi-perspective 2D images, as required for fast 3D imaging.

In state-of-the-art plenoptic cameras [32], directional detection is achieved by inserting a microlens array between the main lens and the sensor of a standard digital camera (see Figure 1a). The sensor acquires composite information that allows identification of both the object point and the lens point where the detected light is coming from. However, the image resolution decreases with inverse proportionality to the gained directional information, for both structural (use of a microlens array) and fundamental (Gaussian limit) reasons; plenoptic imaging at the diffraction limit is thus considered to be unattainable in devices based on simple intensity measurement [5].

**Figure 1.** (**a**) the scheme of a conventional plenoptic imaging (PI) device: the image of the object is focused on a microlens array, while each microlens focuses an image of the main lens on the pixels behind. Such a configuration entails a loss of spatial resolution proportional to the gain in directional resolution; (**b**) shows the scheme of a correlation plenoptic imaging (CPI) setup, in which directional information is obtained by correlating the signals retrieved by a sensor on which the object is focused with a sensor that collects the image of the light source. The image in (**a**) is reproduced with the permission from Ref. [16], copyright American Physical Society, 2017.

Recently, the INFN group involved in Qu3D has proposed a novel technology, named Correlation Plenoptic Imaging (CPI), that enables overcoming the resolution drawback of current plenoptic devices, while keeping their advantages in terms of refocusing capability and 3D reconstruction [1,16,33]. CPI is based on either intensity correlation measurement or photon coincidence detection, according to the light source: actually, CPI can be based on the spatio-temporal correlations characterizing both chaotic sources [16,33] and entangled photon beams [34] to encode the spatial and directional information on two disjoint sensors, as shown in Figure 1b. CPI with chaotic light is based on the measurement of the correlation function

$$
\Gamma(\mathfrak{p}\_{\mathfrak{a}}, \mathfrak{p}\_{\mathfrak{b}}) = \langle l\_{\mathfrak{a}}(\mathfrak{p}\_{\mathfrak{a}}) l\_{\mathfrak{b}}(\mathfrak{p}\_{\mathfrak{b}}) \rangle - \langle l\_{\mathfrak{a}}(\mathfrak{p}\_{\mathfrak{a}}) \rangle \langle l\_{\mathfrak{b}}(\mathfrak{p}\_{\mathfrak{b}}) \rangle,\tag{1}
$$

where ... denotes the average on the source statistics, and *Ijρ<sup>j</sup>* (*j* = *a*, *b*) are the intensities propagated by the beam *j* and registered in correspondence of point *ρ<sup>j</sup>* on the sensor D*j*. Experimentally, the statistical averages are replaced by time averages, obtained by retrieving a collection of frames, simultaneously acquired by the two detectors. In CPI devices, the correlation function encodes combined information on the distribution of light on two reference planes, one of which corresponds to the "object plane" that would be focused on the sensor in a standard imaging setup, placed at a distance *so* from the focusing element. In general, given an object placed at a distance *s* from the focusing element, and characterized by the light intensity distribution A(*ρ*), its images are encoded in the function Γ(*ρa*, *ρb*), in the geometrical-optics limit, as

$$
\Gamma(\rho\_{a\prime}\rho\_{b\prime}) \sim \mathcal{A}^n \left(\frac{s}{s\_o}\frac{\rho\_a}{M} + \left(1 - \frac{s}{s\_o}\right)\frac{\rho\_b}{M\_L}\right), \tag{2}
$$

where *M* and *ML* are the magnifications of the images of the reference object plane and of the focusing element, respectively, while the power *n* is equal to 1 or 2, according to whether the object lies in only one [16,33,35] or both [36,37] optical paths.

Experimental CPI based on pseudo-thermal light is shown in Figure 2c, where both the acquired out-of-focus image and the corresponding refocused image are shown [16]. It was demonstrated in this proof of principle that CPI is characterized by diffraction-limited resolution on the object plane focused on the sensor. Details on the resolution limits are shown in Figure 2a, where one can observe an even more striking effect: thanks to its intrinsic coherent nature, CPI enables an unprecedented combination of resolution and DOF [16]. However, the low-noise sCMOS camera employed in the experiment, working at 50 fps at full resolution, requires several minutes to acquire 30,000 frames used to reconstruct the plenoptic correlation function, and a standard workstation has taken over 10 h for elaborating the acquired data and perform refocusing. The resulting image was also rather noisy, due to the well-known resolution vs. noise compromise of chaotic light ghost imaging that keeps also affecting CPI [35].

**Figure 2.** (**a**) shows the resolution limits, as a function of the longitudinal position, of the image of a double-slit mask with center-to-center distance *d* equal to twice the slit width; here, CPI outperforms both conventional imaging and standard PI with 3 × 3 directional resolution. The evident asymmetry of the CPI curve is due to the existence of two planes in which the object is focused: one at *zb* = *za* and one at *zb* = 0 (see [16]). Plots in (**b**) show a result of a simulation: the target is moved from the focused plane (top left) to an out-of-focus plane (top right). Starting from this position, we show the results of PI refocusing with 3 × 3 directional resolution (bottom left) and the CPI refocusing (bottom right); (**c**) shows the results of an experiment [16] in which the standard image of a triple slit was completely blurred (top), while the image obtained by CPI (bottom) was made fully visible by exploiting information on light direction. Plots in (**c**) are reproduced with the permission from Ref. [16], copyright American Physical Society, 2017.

We are addressing these issues by employing two kinds of sources:

• Chaotic light sources, such as pseudothermal light, natural light, LEDs and gas lamps, and even fluorescent samples, operated either in the high-intensity regime or in

the "two-photon" regime, in which an average of two photons per coherence area propagates in the setup. Chaotic light sources are well known to be characterized by EPR-like correlations in both momentum and position variables [38,39], to be exploited in an optimal way to retrieve an accurate plenoptic correlation function in the shortest possible time. In order to efficiently retrieve spatio-temporal correlations, tight filtering of the source can be necessary to match the filtered source coherence time with the response time of the SPAD arrays that can be as low as 1 ns. Alternatively, pseudorandom sources with a controllable coherence time, made by impinging laser light on a fast-switching digital micromirror device (DMD), can be employed. Interestingly, recent studies have shown that, in the case of chaotic light illumination, the plenoptic properties of the correlation function do not need to rely strictly on ghost imaging: correlations can be measured between any two planes where ordinary images (see Figure 1b) are formed [35]. This discovery has led to the intriguing result that the SNR of CPI improves when ghost imaging of the object is replaced by standard imaging [40]. In particular, excellent noise performances are expected in the case of images of birefringent objects placed between crossed polarizers. This kind of source is particularly relevant in view of applications in fields like biomedical imaging (cornea, zebrafish, invertebrates, biological phantoms such as starch dispersions), security (distance detection, DOF extension), and satellite imaging.

• Momentum–position entangled beams, generated by spontaneous parametric downconversion (SPDC), which have the potential to combine QPI with sub-shot noise imaging [41], thus enabling high-SNR imaging of low-absorbing samples, a challenging issue in both biomedical imaging and security.

The design of both quantum plenoptic devices are currently undergoing optimization by implementing a novel protocol that enables further mitigating the resolution vs. DOF compromise with respect to the one shown in Figure 2a: this protocol is based on the observation that, for any given resolution, the DOF can be maximized by correlating the standard images of two arbitrary planes, chosen in the surrounding of the object of interest, instead of imaging the focusing element [36]. Moreover, we are investigating the possibility to merge quantum plenoptic imaging with the measurement protocols developed in the context of differential ghost imaging [42].

In the following section, we will discuss all the technical solutions, both on the hardware and on the software side that we are investigating to speed up the process of generating a correlation plenoptic image. The problematic aspects can be clarified by a description of the practical process of CPI imaging that is made of the following steps:

1. Collecting *Nf* synchronized pairs of frames from a sensor D*a*, with resolution *Nax* × *Nay*, and D*b*, with resolution *Nbx* × *Nby*. Synchronization of two separate digital sensors entails technical complications that can be overcome by using two disjoint parts of the same sensor [16]. The total acquisition time is

$$T = N\_f(\tau\_{\text{exp}} + \tau\_{\text{d}}),\tag{3}$$

where *τ*exp is the frame exposure time that must be shorter than the coherence time of impinging light in order to exploit the maximal information on intensity fluctuations, while *τ*<sup>d</sup> is the dead time between subsequent frames, usually fixed by the employed sensors.


$$\Gamma(\boldsymbol{\rho}\_{a\prime}\boldsymbol{\rho}\_{b}) \simeq \frac{1}{N\_{f}} \sum\_{k=1}^{N\_{f}} I\_{a}^{(k)}(\boldsymbol{\rho}\_{a})I\_{b}^{(k)}(\boldsymbol{\rho}\_{b}) - \frac{1}{N\_{f}} \sum\_{k=1}^{N\_{f}} I\_{a}^{(k)}(\boldsymbol{\rho}\_{a}) \frac{1}{N\_{f}} \sum\_{k'=1}^{N\_{f}} I\_{b}^{(k')}(\boldsymbol{\rho}\_{b}), \tag{4}$$

where *I* (*k*) *<sup>j</sup>* (*ρj*), with *j* = *a*, *b*, is the intensity measured in frame *k* in correspondence of the pixel on the detector D*<sup>j</sup>* centered on the coordinate *ρj*. In this way, an information initially encoded in *Nf*(*NaxNay* + *NbxNby*) numbers is used to reconstruct a correlation function determined by *NaxNayNbxNby* values.

The accuracy of the correlation function estimate (4) increases with the number of frames like *Nf* [40]. However, increasing *Nf* also linearly extends the total acquisition time *T*. Therefore, the combination of (1) fast and low-noise sensors, (2) methods to extract a good quality signal from smaller number of frames, and (3) tools for efficient data transfer and elaboration, is crucial in order to speed-up the acquisition process and make the CPI technology ready for real-world applications.

#### **3. Hardware Speedup: Advanced Sensors and Ultra-Fast Computing Platforms**

To improve the performances of CPI in terms of acquisition speed and data elaboration time, we are employing dedicated advanced sensors and ultra-fast computing platforms. In this section, we describe the details of the implementation and the perspectives on these fields.

#### *3.1. SPAD Arrays as High-Resolution Time-Resolved Sensors*

A relevant part of the speedup that we are seeking is determined by replacing commercial high-resolution sensors, like scientific CMOS and EMCCD cameras, with sensors based on cutting-edge technology such as single-photon avalanche diode (SPAD) arrays. A SPAD is basically a photodiode which is reversely biased above its breakdown voltage, so that a single photon which impinges onto its photosensitive area can create an electron– hole pair, triggering in turn an avalanche of secondary carriers and developing a large current on a very short timescale (picoseconds) [17,18]. This operation regime is known as Geiger mode. The SPAD output voltage is sensed by an electronic circuit and directly converted into a digital signal, further processed to store the binary information that a photon arrived, and/or the photon time of arrival. In essence, a SPAD can be seen as a photon-to-digital conversion device with exquisite temporal precision. SPADs can also be gated, in order to be sensitive only within temporal windows as short as a few nanoseconds, as shown in Figure 3. Individual SPADs can nowadays be used as the building blocks of large arrays, with each pixel circuit containing both the SPAD and the immediate photon processing logic and interconnect. Several CMOS processes are readily available and allow for tailoring both the key SPAD performance metrics and the overall sensor or imager architecture [43,44]. Sensitivity and fill-factor have for some time lagged behind those of their scientific CMOS or EMCCD counterparts but have been substantially catching up in recent years.

Based on the requirements of QPI, we have chosen to employ the SwisSPAD2 array developed by the AQUA laboratory group of EPFL, characterized by a 512 × 512 pixel resolution (see Figure 4), which is one of the widest and most advanced SPAD arrays to date [20,22]. The sensor is internally organized as two halves of 256 × 512 pixels to reduce load and skew on signal lines and enable faster operation. It is a purely binary gated imager, i.e., each pixel records either a 0 (no photon) or a 1 (one or more photons) for each frame, with basically zero readout noise. The sensor is controlled by an FPGA generating the control signals for the gating circuitry and readout sequence and collecting the pixel detection results. In the FPGA, the resulting one-bit images can be further processed, e.g., accumulated into multi-bit images, before being sent to a computer/GPU for analysis and storage. The maximum frame rate is 97.7 kfps, and the native fill factor of 10.5% can be improved by 4–5 times, for collimated light, by means of a microlens array [21] (higher values are expected from simulation after optimization); the photon detection probability is 50% (25%) at 520 nm (700 nm) and 6.5 V excess bias. The device is also characterized by low noise (typically less than 100 cps average Dark Count Rate per pixel at room temperature, with a median value about 10 times lower) and advanced circuitry for nanosecond gating. A detailed comparison of SwissSPAD2 with other large-format CMOS SPAD imaging cameras is presented in Ref. [22].

**Figure 3.** SwissSPAD2 gate window profile. The transition times and the gate width are annotated in the figure. The gate width is user-programmable, and the minimum gate width in the internal laser trigger mode is 10.8 ns. The image is reproduced with permission from Ref. [45], copyright The Authors, published by IOP Publishing Ltd., 2020.

**Figure 4.** SwissSPAD2 photomicrograph (**left**) and pixel schematics (**right**). The pixel consists of 11 NMOS transistors, 7 with thick-oxide, and 4 with thin-oxide gate. The pixel stores a binary photon count in its memory capacitor. The in-pixel gate defines the time window, with respect to a 20 MHz external trigger signal, in which the pixel is sensitive to photons. The image is reproduced with permission from Ref. [22], copyright IEEE, 2018.

Currently, we are using the available version of SwissSPAD2, with a 512 × 256 pixel resolution, to generate sequences of frames and store them in an on-board 2 GB memory, before transferring them to a computer by a standard USB3 connection, which can be done using existing hardware. We are integrating this sensor in the prototype of chaoticlight base quantum plenoptic camera in a way that two disjoint halves of the sensor (of 256 × 256 pixels each) are used for retrieving the images of the two reference planes. The high speed of this sensor is expected to reduce the acquisition time of quantum plenoptic images by two order of magnitudes with respect to the first CPI experiment [16], in which the region of interest on the sensor was made of 700 × 700 pixels for the spatial measurement side, and 600 × 600 for the directional measurement side; 2 × 2 binning on both sides during acquisition and a subsequent 10 × 10 binning on the directional side led to the effective spatial resolution of 350 pixels, and angular resolution of 25 × 25 pixels.

We are also working towards several further optimizations of the sensor system, e.g., by developing gating for noise reduction in QPI devices. In order to employ the full sensor (512 × 512 pixels), we are implementing a synchronization mechanism for a pair of imagers by means of two FPGAs, so as to operate on a common time-base at the nanosecond

level (to this end, two control circuits shall operate from a single clock and have a direct communication link). Finally, the SwissSPAD2 arrays are being integrated with a fast communication interface in order to speed up data transfer and make it possible to deliver, in a sustained way, full binary frame sequences to a GPU. The latter will run advanced algorithms for data pre-processing, image reconstruction and optimization, as we will discuss in more detail in the next subsection.

#### *3.2. Computational Hardware Platform*

In a QPI device integrated with SwissSPAD2, the acquired data rate for a single frame acquisition can be estimated to 26 Gb/s, which is beyond the reach of standard data buses. This poses great challenges in both hardware and software design. Our approach is to employ the expertise of PKH to seek careful design of electronic interconnections (buses) between sensor control electronics and processing device, theoretical refinement, and optimization of algorithms (e.g., compressive sensing [46,47]), porting to an efficient computational environment, and design of a specific acquisition electronics for optimizing data flow from the light sensors to a dedicated processing system, able to guarantee the required computational performance (e.g., exploiting GPUs or FPGA) [23,24].

The introduction of an embedded data acquisition- and processing board, integrating a GPU, aims at data pre-processing, thus significantly reducing the amount of data to be transferred to (and saved on) an external workstation. GPUs exploit a highly parallel elaboration paradigm, enabling to design algorithms that run in parallel on hundreds or thousands of cores and to make them available on embedded devices. A great advantage of GPUs is programmability: many standard tools exist (e.g., OpenCL and CUDA) that allow fast and efficient design of complex algorithms that can be injected on the fly in the GPU memory for accomplishing tasks ranging from simple filtering to advanced machine learning. Efforts are made to design a processing matrix, so that each line and column of the sensor will be managed by a dedicated portion of the heterogeneous processing platform (CPU/GPU/FPGA): the pixel series processor. Those dedicated units will be interconnected to one another to cooperate for implementing algorithms that require distributed processing on a very small scale.

The embedded acquisition-processing board is designed to best fit to SPAD array and SW application needs. The optimal system design will be evaluated by considering theoretical algorithms, engineered SW implementations, HW set-up, and HW/SW tradeoff. We will identify a preliminary set of possible configurations and perform a trade-off by comparing overall performances, considering the requirements in data quality, processing speed, costs, complexity, etc.

Based on the challenging objectives to be achieved, a preliminary analysis based on COTS (Commercial-Off-The-Shelf) solutions was performed, in order to identify a set of accelerating devices addressing QPI requirements in terms of computing capability and portability. The option offered by the NVIDIA Jetson Xavier AGX board shown in Figure 5 is considered a promising candidate to achieve our goal. This device indeed offers an encouraging performance/integration ratio with low power usage and very interesting computing capabilities. Its main characteristics are reported in Table 1.

Considering the listed HW capabilities, despite the ARM processor having a limited computing power when compared to high-end desktop processors, it allows for leveraging multi-core capabilities for implementing that part of the code that will feed the quite powerful GPU device on board. In addition, the foreseen optimization strategy will require a dedicated implementation able to consider the maximum amount of memory of 16 GB available, shared with the GPU. In addition, the solution to be developed should take into account the bandwidth of about 136 GB/s of the on-board memory, which may represent a limiting factor when GPU and CPU exchange buffers. An implementation based on the CUDA framework—over OpenCL or other technologies—will be preferred to best use the NVIDIA device. Finally, given the capability of the NVDLA Engines to perform multiplications and accumulations in a very fast way, we consider it interesting to perform

an assessment of how to exploit these devices for implementing the QPI-specific correlation functions and/or other multiplication/sum intensive computations.

**Figure 5.** Qu3D scenario for GPU parallel processing based on NVIDIA Volta architecture as provided by an NVIDIA Jetson AGX Xavier device.


**Table 1.** NVIDIA Jetson Xavier AGX: Technical Specifications.

Along with the assessment of a dedicated HW solution, further optimizations applicable at the algorithmic level have been considered in order to enrich the engineered device with a highly customized algorithmic workflow able to exploit the peculiarities of the CPI technique and its related input data. More specifically, we are analyzing those steps of QPI processing that appear as more computationally demanding, thus representing a bottleneck for performances. To this end, tailored reshaping operations applied over the original three-dimensional multi-frame structure of the input data were explored, to facilitate the development of a parallelized elaboration paradigm for evaluating the CPI-related correlation function. In addition, the peculiar feature of the input dataset to be acquired by SPAD sensors as one-bit images will be valued through dedicated implementations able to gain from intensive multiplication/sum math among binary-valued variables.

#### **4. Quantum and Classical Image-Processing Algorithms**

Further reduction of the acquisition speed and the optimization of the elaboration time is addressed by exploiting dedicated quantum and classical image processing, as well as novel mathematical methods coming on one hand from compressive sensing, and on the other hand from quantum tomography and quantum Fisher information.

#### *4.1. Compressive Sensing*

In order to reduce the amount of required data (at present 103–104 frames) by at least one order of magnitude, we are exploring different approaches. First, we are investigating the opportunity of implementing compression techniques for improving bus bandwidth utilization, thus acting on data transfer optimization. Data compression may also rely on manipulation of raw input data to determine only the relevant information in a sort of information bottleneck paradigm, in which software nodes in a lattice provide their contribution to a probable reconstruction of the actual decompressed data. This is very similar to artificial neural network structures, where the network stores a representation of a phenomenon and returns a response based on some similarity rating versus the observed data. Moreover, the availability of advanced processing technologies allows the investigation and implementation of alternative techniques able to exploit the sparsity of the retrieved signal to reconstruct information from a heavily sub-sampled signal, by compressive sensing techniques [46,47].

As in other correlation imaging techniques, such as ghost imaging (GI), in the CPI protocol, the object image is reconstructed by performing correlation measurements between intensities at two disjoint detectors. Katz et al. [48] demonstrated that conventional GI offers the possibility to perform compressive sensing (CS) boosting the recovered image quality.

CS theory asserts the possibility of recovering the object image from far fewer samples than required by the Nyquist–Shannon sampling theorem, and it relies on two main principles: sparsity of the signal (once expressed in a proper basis) and incoherence between the sensing matrix and the sparsity basis.

In conventional GI, the transmission measured for each speckle pattern represents a projection of the object image and CS finds the sparsest image among all the possible images consistent with the projections. In practice, CS reconstruction algorithms solve a convex optimization problem, looking for the image which minimizes the *L*1−norm in the sparse basis among the ones compatible with the bucket measurements, see Refs. [49–52] for a review.

We are developing a novel protocol reducing the number of measurements required for image recovery by an order of magnitude. Once properly refocused, a single acquisition can be fed into the compressive sensing algorithm several times thus exploiting the plenoptic properties of the acquired data and increasing the signal-to-noise ratio of the final refocused image. We tested the CS-CPI algorithm with numerical simulations, as summarized in Figure 6. In (a), a double-slit mask is reconstructed by correlations measurements considering *N* = 6000 frames; in (b), the standard reconstruction is repeated considering only the 10% of available frames, while, in (c), the CS reconstruction using the same reduced set of measured data. In addition to a data-fidelity term corresponding to a linear regression, we penalized the *L*1−norm of the reconstructed image to account for its sparsity in the 2*D* − *DCT* domain. The resulting optimization problem is known as the LASSO (Least Absolute Shrinkage and Selection Operator) [53]. We employed the coordinate descent algorithm to efficiently solve it, and we set the regularization parameter, controlling the degree of sparsity of the estimated coefficients, by cross-validation. In this proof-of-concept experiment, we simply use Pearson's correlation coefficient to measure the similarity between the reconstructions obtained using the restricted dataset and the image obtained considering all the *N* = 6000 frames.

**Figure 6.** (**a**) double-slit image reconstruction obtained by correlations measurements, considering *N* = 6000 frames and two detectors characterized by a 128 × 128 and a 10 × 10 pixel resolution; (**b**) the standard reconstruction is repeated considering only the 10% of the available frames, chosen randomly; (**c**) compressive sensing reconstruction using the same dataset as in (**b**). While in the first case, the Pearson's correlation coefficient is *rred* = 0.55; in the latter case, the coefficient is increased to *rCS* = 0.81.

#### *4.2. Plenoptic Tomography*

Novel reconstruction algorithm with the advantage of real 3D image lack of artifacts is based on the idea of recasting plenoptic imaging as an absorption tomography. The classical refocusing algorithm is based on superimposing images of the 3D scene from different viewpoints, which necessarily results in the contribution of out-of-focus parts of the scene, the effect responsible for the existence of the blurred part of 3D image reconstruction. The tomography approach is based on a different principle and provides sufficient axial and transversal resolution without artifacts. For this purpose, the object space is divided into voxels and the goal is to reconstruct absorption coefficients of each voxel. The measured correlation coefficient of a pair of points from the correlated planes is transformed to be linearly proportional to the attenuation along the ray path connecting those two points. Classical inverse Radon transform can be used in this scenario to obtain a 3D image, but using the Maximum Likelihood absorption tomography algorithm [54] further enhances the quality of the tomography reconstruction and performs well even for a small range of projections angels. In fact, advanced tomographic reconstruction algorithms based on the Maximum Likelihood principle are more resistant to noise and require fewer acquisitions, for a given precision, in comparison to standard tomographic protocols. We are investigating special tools to deal with informationally incomplete detection schemes for very high resolution and optimal methods for data analysis based on convex programming tools.

#### *4.3. Quantum Tomography and Quantum Fisher Information*

A further quantum approach to image analysis and detection schemes will be employed to achieve super-resolution (or, eventually, for maintaining the desired resolution and speeding up acquisition and elaboration of the quantum plenoptic images) and to compare correlation plenoptic detection scheme to the ultimate quantum limits. The basic concept underpinning the Fisher Information super-resolution imaging is the formal mathematical analogy between the classical wave optics and quantum theory, which makes it possible to apply the advanced tools of quantum detection and estimation theory to classical imaging and metrology [29,55]. The advantage of this approach lies in the ability to quantify the performance of the imaging setup based on rigorous statistical quantities. Inspired by the quantum theory of detection and estimation, quantum Fisher information, a quantity connected with the ultimate limits allowed by nature, is computed for simple 3D imaging scenarios like localization and resolution of two points in the object 3D space [30,56–58]. For example, one might be interested in measuring the separation of two point-like sources and seeking the optimal detection scheme, extracting the maximum amount of information about this parameter. We shall thus employ quantum Fisher information to design optimal measurement protocols within the quantum plenoptic devices, able to extract specific relevant information, for enhanced resolution, with a minimal number of acquisitions. The question of the minimal number of detections over the correlation planes which achieves acceptable reconstruction quality is of relevance because of the direct connection to the amount of processed data and can lead to reduced time for data processing.

#### **5. Perspectives**

In the context of the Qu3D project, all the developments and technologies presented in the previous sections will be integrated into the implementation of two quantum plenoptic imaging (QPI) devices, namely,


Our objective is to achieve, in both devices, high resolution, whether diffractionlimited or sub-Rayleigh, combined with a DOF larger by even one order of magnitude compared to standard imaging. The science and technology developed in the project will contribute to establishing a solid baseline of knowledge and skills for the development of a new generation of imaging devices, from quantum digital cameras enhanced by refocusing capability to quantum 3D microscopes [37] and space imaging devices.

#### **6. Conclusions**

We have presented the challenging research directions we are following to achieve practical quantum 3D imaging: minimizing the acquisition speed without renouncing to high SNR, high resolution, and large DOF [1–5]. Our work represents a significant advance with respect to the state-of-the-art of both classical and quantum imaging, as it enhances the performances of plenoptic imaging and dramatically speeds up quantum imaging, thus facilitating the real-world deployment of quantum plenoptic cameras.

This ambitious goal will be facilitated by working toward the extension of the reach of quantum imaging to other fields of science, and opening the way to new opportunities and applications, including the prospects of offering new medical diagnostic tools, such as 3D microscopes for biomedical imaging, as well as novel devices (quantum digital cameras enhanced by 3D imaging, refocusing, and distance detection capabilities) for security, space imaging, and industrial inspection. The collaboration crosses the traditional boundaries between the involved disciplines: quantum imaging, ultra-fast cameras, low-level programming of GPU, compressive sensing, quantum information theory, and signal processing.

**Author Contributions:** Conceptualization, M.D., Z.H., J.R., S.B., P.M., E.C. and C.B.; methodology, ˇ F.V.P., G.M., M.P. (Martin Paur), B.S., S.B., P.M., E.C. and C.B.; software, F.S. (Francesco Scattarella), G.M., L.M., M.P. (Michal Peterek), M.I. (Michele Iacobellis), I.P., F.S. (Francesca Santoro), P.M., A.U. and M.W.; firmware, P.M., A.U. and M.W.; validation, F.D.L., G.M., M.I. (Michele Iacobellis), P.M., A.U. and M.W.; formal analysis, F.D.L., M.I. (Michele Iacobellis), I.P. and F.S. (Francesca Santoro); setup implementation, P.M., A.U. and M.W.; setup design, F.D.L., D.G. and S.V.; writing—original draft preparation, F.V.P., I.P., F.S. (Francesca Santoro), S.B. and C.B.; writing—review and editing, F.D.L., F.V.P., I.P., F.S. (Francesca Santoro) and C.B.; supervision, M.D., A.G., F.V.P., C.A., L.A., E.C. and C.B.; project administration, M.D., C.A. and L.A.; funding acquisition, M.D., C.A., M.I. (Maria Ieronymaki) and C.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** Project Qu3D is supported by the Italian Istituto Nazionale di Fisica Nucleare, the Swiss National Science Foundation (grant 20QT21\_187716 "Quantum 3D Imaging at high speed and high resolution"), the Greek General Secretariat for Research and Technology, the Czech Ministry of Education, Youth and Sports, under the QuantERA program, which has received funding from the European Union's Horizon 2020 research and innovation program.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Relevant data are available from the authors upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:



#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com