Foundations of Photonic Quantum Computation. Part 1. Introduction to Quantum Operators

1  Introduction

1.1 Photonic Quantum Computing

Simulating quantum physics on a computer is widely believed to be an algorithmically difficult task. This was already the case as early as 1980, which led Manin [1,2] to briefly mention the idea of using quantum systems to build quantum automata. The same year, Benioff [3] introduced quantum Turing machines. One of the problems that arise when one wants to simulate quantum physics accurately is the amount of memory needed to fully represent the quantum state of the system [4]. This difficulty in simulating quantum systems led Richard Feynman to imagine a “computer” made of and for quantum systems in his foundational article “Simulating Physics with Computers” [5]. The concept was then extended to what is known as gate model quantum computing (QC) by David Deutch [6]. The first advantage of this quantum model of computation was found by leveraging the query model on Deutsch’s problem. QC allows us to check if a given function f:{0,1}→{0,1} is constant (f(0)=f(1)) or if it is balanced (f(0)≠f(1)) in one query whereas its classical counterpart needs two measurements. An extended version was then proposed in [7], achieving an even greater improvement over the classical counterpart.

Building on these early successes, researchers came up with other quantum algorithms which had a lower algorithmic complexity than their classical counterparts. The canonical and most famous examples are Grover’s [8] search algorithm and Shor’s factoring algorithm [9].

Photonic QC is one of many physical realizations of a quantum computer. The first proposal came from [10], which introduced “dual-rail encoding,” making it possible to encode qubits into Fock states [11]. This led Knill, Laflamme, and Milburn to develop an entire protocol for physically realizing a gate-model quantum computer using optical components, known as the KLM protocol [12]. This protocol was further developed in the following years in [13,14]. This model was refined in [15] into a more complex model of photonic QC, called one-way QC or measurement-based QC, which lies outside the scope of this work.

Photons have proven to be good candidates for the main applications of quantum information. First, they have been used to simulate Hamiltonian dynamics with good precision and high success probability on a 4-qubit photonic silicon chip [16]. Photonic quantum information has applications to cybersecurity with many examples of quantum key distribution being achieved using photons [17–19]. In the field of optimization, although no advantage has been found yet, several promising papers aim to reduce the qubit count required to use photonic quantum computers to solve hard optimization problems [20–22].

Another reason to consider photons a powerful quantum resource is that Scott Aaronson and Anton Arkhipov proposed a protocol for demonstrating quantum supremacy, called boson sampling, which can be implemented using photons and linear optical components [23]. In 2020, ref. [24] reported a claim of quantum computational advantage using this protocol.

1.2 Understanding Modes and Defining Qubits from Modes

To create qubits in a photonic quantum computer, the first method proposed was path encoding, introduced by [10]. The locations of the photons—referred to as spatial modes—define the state of the qubit. Spatial modes can, for instance, correspond to optical fibers. For convenience, throughout this chapter, we use the term mode to refer specifically to spatial mode.

To describe a system composed of photons and modes. We use the notation |n0,n1,…,nm−1⟩01…(m−1) where 0,1,…,(m−1) denote the spatial modes occupied by n0,n1,…,nm−1 photons. This type of states is referred to as Fock state.

To define a Qubit, let us now consider a specific case where a single photon is distributed over two modes. The definition of the two basic states of a qubit in the dual rail encoding relies on the location of the said photon.

|0⟩=|1,0⟩01

|1⟩=|0,1⟩01

When the photon is in mode 0, this corresponds to the Fock state |1,0⟩01 sometimes written simply as |1,0⟩ (although this notation should be avoided to prevent confusion between qubit states and Fock states).

Remark

When possible and in cases where ambiguity may arise, we use the notation |n0n1…nm−1⟩f to indicate explicitly that the set of modes corresponds to a Fock state.

It is also common to omit the commas between each component of the state, so that |n0,n1,…,nm−1⟩f is often written as |n0n1…nm−1⟩f.

The Fock state |1,0⟩01is encoding the qubit state |0⟩ and a representation is provided in Fig. 1.

Figure 1: Qubit state |0⟩=|10⟩01.

Similarly, the Fock state |01⟩01 is encoding the qubit state |1⟩ and a representation is provided in Fig. 2.

Figure 2: Qubit state |1⟩=|01⟩01.

This encoding is called the dual mode encoding. But there exist many other possibilities to encode qubits, qudits or even continuous variable quantum information! [25].

1.3 Creation and Annihilation Operators for the Modes

The creation and annihilation operators associated with a mode a are denoted a^a† and a^a, respectively. Their actions are as follows:

•   The creation operator a^0† allows the transition from |0⟩0 to |1⟩0 in spatial mode 0. For example, a^0† transforms the vacuum state |00⟩01 into |10⟩01.

•   The annihilation operator a^0 allows the transition from |1⟩0 to |0⟩0. For example, a^0 transforms |10⟩01 into the vacuum state |00⟩01.

These operators correspond to the emission or annihilation of a photon in a given mode. They play a crucial role when describing the operation of complex optical components.

2  The Generic Beam Splitter Component

2.1 Physical Description of the Action of the Beam Splitter

A beam splitter is an optical component used in photonic circuits. Let Ea and Eb denote the amplitudes of the two input electric fields, and Ec and Ed the amplitudes of the two output fields, as illustrated in Fig. 3. Let rac and rbd be the reflection coefficients, tad and tbc the transmission coefficients. We assume that the magnitudes are restricted to |rac|,|rbd|,|tad|,|tbc|∈[0, 1].

Figure 3: Representation of a beam splitter with two input modes and two output modes.

The output field amplitudes can be expressed in terms of the input amplitudes through the so-called continuity relations:

Ec=racEa+tbcEb(1)

Ed=tadEa+rbdEb(2)

where

rac: is the contribution of Ea found at the output Ec

rbd: is the contribution of Eb found at the output Ed

tad: is the contribution of Ea found at the output Ed

tbc: is the contribution of Eb found at the output Ec

The conservation of energy (proportional to the squared modulus of Ex) leads to:

|Ea|2+|Eb|2=|Ec|2+|Ed|2(3)

Using (1) and (2), the action of the beam splitter can be described in matrix form based on the field amplitudes of each mode, as shown in Fig. 4.

Figure 4: Action of a beam splitter.

Let us note:

U=[ractbctadrbd]

If the entire input field is injected into mode a, the input amplitudes satisfy

|Ea|2=1,|Eb|2=0,

which is consistent with the normalization condition ∣Ea∣2+∣Eb∣2=1, then we have:

(ractbctadrbd).(Ea0)=(EcEd)

It follows that:

rac.Ea=Ec(4)

tad.Ea=Ed(5)

Staying in the case of one single input (here in a), we have |Eb|2=0, applying to (3) gives:

|Ea|2=|Ec|2+|Ed|2(6)

because (4) and (5):

|Ea|2=|racEa|2+|tadEa|2

So

|Ea|2=(|rac|2+|tad|2).|Ea|2

Because we consider that we have an input in a, we have Ea≠0 and we can divide both terms of the equation by |Ea|2

1=|rac|2+|tad|2(7)

By the applying the same process this time with Ea=0 and Eb≠0:

(ractbctadrbd).(0Eb)=(EcEd)

One would obtain:

1=|rbd|2+|tbc|2(8)

Given (1) and (2), it is possible to calculate:

|Ec|2=|racEa+tbcEb|2=|racEa|2+|tbcEb|2+racEa(tbcEb)∗+tbcEb(racEa)∗

Similarly

|Ed|2=|tadEa|2+|rbdEb|2+tadEa⋅rbd∗Eb∗+rbdEb⋅tad∗Ea∗

Adding both equations together:

|Ec|2+|Ed|2=(|rac|2+|tad|2)|Ea|2+(|rbd|2+|tbc|2)|Eb|2+racEa⋅(tbcEb)∗+tbcEb⋅(racEa)∗+tadEa⋅rbd∗Eb∗+rbdEb⋅tad∗Ea∗=(|rac|2+|tad|2)⏟=1(see ((7)))|Ea|2+(|rbd|2+|tbc|2)⏟|Eb|=12+(tadrbd∗+ractbc∗)EaEb∗+(rbdtad∗+tbcrac∗)Ea∗Eb

So

|Ec|2+|Ed|2=|Ea|2+|Eb|2+(tadrbd∗+ractbc∗)EaEb∗+(rbdtad∗+tbcrac∗)Ea∗Eb

Reminder

According to (3), we should have:

|Ea|2+|Eb|2=|Ec|2+|Ed|2

By substituting |Ec|2+|Ed|2 using the above relation:

|Ea|2+|Eb|2=|Ea|2+|Eb|2+(tadrbd∗+ractbc∗)EaEb∗+(rbdtad∗+tbcrac∗)Ea∗Eb

We get:

(tadrbd∗+ractbc∗)EaEb∗+(rbdtad∗+tbcrac∗)Ea∗Eb=0

Since this expression must hold for any choice of Ea and Eb, it implies the two following relations:

(tadrbd∗+ractbc∗)=0(9)

(rbdtad∗+tbcrac∗)=0(10)

Theorem 1: U is unitary.

Proof: Consider the matrix:

U†U=[rac∗tad∗tbc∗rbd∗]⋅[ractbctadrbd]

U†U=[racrac∗+tadtad∗rac∗tbc+tad∗rbdtbc∗rac+rbd∗tadtbc∗tbc+rbdrbd∗]

U†U=[|rac|2+|tad|2rac∗tbc+tad∗rbdtbc∗rac+rbd∗tad|tbc|2+|rbd|2]

Using (7)–(10), we directly get:

U†⋅U=[1001]

Therefore, the beam splitter’s operation is unitary. □

2.2 Parameterization of the Beam Splitter Matrix U

The gate U representing the action of a beam splitter can then be denoted as BS:

BS=(ractbctadrbd)

There are several conventions for the beam splitter, some of which belong to U(2) but not to SU(2). The goal of this section is to understand how those different conventions appear and why different forms exist.

Theorem 2: ϕad−ϕbd−ϕac+ϕbc=π with ϕij is the phase of the exponential form of tij.

Proof: As a starting point, let us take the relation demonstrated in the preceding section:

tbc∗.rac+rbd∗.tad=0

As all the coefficients are complex numbers, each one can be represented in its exponential form so we can be written:

|tbc|e−i.ϕbc.|rac|.ei.ϕac+|rbd|.e−i.ϕbd.|tad|.ei.ϕad=0

|tbc|.|rac|.ei.(ϕac−ϕbc)+|rbd|.|tbd|.ei.(ϕad−ϕbd)=0

So

|tbc|.|rac|.ei.(ϕac−ϕbc)=−|rbd|.|tad|.ei.(ϕad−ϕbd)

|tbc|.|rac|=−|rbd|.|tad|.ei.(ϕad−ϕbd−ϕac+ϕbc)

|rac||tad|=−|rbd||tbc|.ei.(ϕad−ϕbd−ϕac+ϕbc)

It is now possible to analyze:

|rac||tad|=−|rbd||tbc|.ei.(ϕad−ϕbd−ϕac+ϕbc)

All magnitudes are positive real numbers, so ei.(ϕad−ϕbd−ϕac+ϕbc) can only take the value −1.

From this we find the first important relation:

ϕad−ϕbd−ϕac+ϕbc=π(11)

□

Theorem 3: r2+t2=1 with t:=|tad|=|tbc| and r:=|rac|=|rbd|.

Proof: To find the second important relation, let us start with:

|rac||tad|=|rbd||tbc|

Squaring both sides of the equation gives:

|rac|2|tad|2=|rbd|2|tbc|2

Using (7) and (8):

1−|tad|2|tad|2=1−|tbc|2|tbc|2

i.e.,

1|tad|2−1=1|tbc|2−1

It follows that

|tad|=|tbc| and we note t:=|tad|=|tbc|

The equation |rac||tad|=|rbd||tbc| can be rewritten |rac|t=|rbd|t which implies that:

|rac|=|rbd| and we note r:=|rac|=|rbd|

Recalling (8)

|rac|2+|tad|2=1

And using the new notation leads to the second important relation:

r2+t2=1(12)

□

Remark on the Physical Meaning of the Coefficients r and t

The coefficients r2=R and t2=T correspond to:

The reflectance R=PreflectedPIncoming

The transmittance T=PTransmittedPIncoming

We can also define:

The absorbance A=PAbsorbedPIncoming where all the P values represent powers.

Energy conservation requires R+T+A=1 and throughout this document, we consider lossless beam splitters (i.e., A=0)

By taking the complex form of its coefficients, the BS operator can be parameterized as:

BS=(|rac|.ei.ϕac|tbc|.ei.ϕbc|tad|.ei.ϕad|rbd|.ei.ϕbd)

And by definition of r and t:

BS=(r.ei.ϕact.ei.ϕbct.ei.ϕadr.ei.ϕbd)

Thus, the BS operator is completely defined by the six parameters:

r,t,ϕac,ϕbc,ϕad,ϕbd

And by applying the two important relations (11) and (12):

r2+t2=1

ϕad−ϕbd−ϕac+ϕbc=π

We would like to reduce the number of parameters of BS.

To do this, we start by using (11):

{r2+t2=1r=|rac|≥0ett>0⇒∃θ∈[0,+π2],r=sin⁡(θ)and t=cos⁡(θ)withθ=arctan⁡(rt)

One can now write:

BS=(r.ei.ϕact.ei.ϕbct.ei.ϕadr.ei.ϕbd)=(sin⁡(θ).ei.ϕaccos⁡(θ).ei.ϕbccos⁡(θ).ei.ϕadsin⁡(θ).ei.ϕbd)

In their paper, ref. [26] propose the change of variables:

ϕt=12(ϕad−ϕbc)

ϕr=12(ϕac−ϕbd−π)

ϕ0=12(ϕad+ϕbc)

To obtain

ϕac=ϕ0+ϕr

ϕad=ϕ0+ϕt

ϕbc=ϕ0−ϕt

ϕbd=ϕ0−ϕr−π

Verification

The verifications for ϕad and ϕbc are straightforward; let us detail ϕac and ϕbd

For ϕac:

ϕ0+ϕr=12(ϕad+ϕbc)+12(ϕac−ϕbd−π)

=12(ϕad+ϕbc+ϕac−ϕbd−π)

ϕ0+ϕr=12(ϕad−ϕac+ϕbc−ϕbd⏟π(11)−π+ϕac+ϕac)

ϕ0+ϕr=ϕac

And for ϕbd:

ϕ0−ϕr−π=12(ϕad+ϕbc)−12(ϕac−ϕbd−π)−π

ϕ0−ϕr−π=12(ϕad−ϕac+ϕbc−ϕbd⏟π(11)+ϕbd+ϕbd+π)−π=(2ϕbd+2π)2−π

ϕ0−ϕr−π=ϕbd

We can therefore rewrite by replacing each of the expressions:

BS=(sin⁡(θ).ei.ϕaccos⁡(θ).ei.ϕbccos⁡(θ).ei.ϕadsin⁡(θ).ei.ϕbd)=(sin⁡(θ).ei.(ϕ0+ϕr)cos⁡(θ).ei.(ϕ0−ϕt)cos⁡(θ).ei.(ϕ0+ϕt)sin⁡(θ).ei.(ϕ0−ϕr−π))

Thus, a common expression for the noiseless beam splitter becomes:

BS1≡BS=ei.ϕ0(sin⁡(θ).ei.ϕrcos⁡(θ).e−i.ϕtcos⁡(θ).ei.ϕt−sin⁡(θ).e−i.ϕr)

We obtained an expression of BS using only 4 parameters:

θ,ϕ0,ϕr and ϕt

2.3 Important Remark on Different Conventions

The definitions of ϕr,ϕt and ϕ0 are not universal; for example, by changing the definition of ϕr:

ϕr=12(ϕac−ϕbd+π)

We would obtain:

BS2(θ,ϕr,ϕt,ϕ0)=ei.ϕ0(−sin⁡(θ2).ei.ϕrcos⁡(θ2).e−i.ϕtcos⁡(θ2).ei.ϕtsin⁡(θ2).e−i.ϕr)

Or alternatively, by defining:

r2+t2=1⇒∃θ∈[0,π],r=cos⁡(θ2) and t=sin⁡(θ2) with θ=arctan⁡(tr)

One could also have obtained:

BS3(θ,ϕr,ϕt,ϕ0)=ei.ϕ0(cos⁡(θ2).ei.ϕrsin⁡(θ2).e−i.ϕtsin⁡(θ2).ei.ϕt−cos⁡(θ2).e−i.ϕr)

In by using BS3 and considering ϕr=ϕt=ϕ0=0 and θ=π2, we obtain the expression of the Hadamard gate.

BS3(π2,0,0,0)=(cos⁡(π4)sin⁡(π4)sin⁡(π4)−cos⁡(π4))=(222222−22)=H

2.4 Quandela Conventions

Note that the conventions presented above are not exhaustive. Perceval, the library proposed by [27], offers by default three conventions.

Let us now consider, in each case, all phases set to zero.

Quandela Convention H:

BSH(θ∈[0;π],0,0,0)=(cos⁡(θ2)sin⁡(θ2)sin⁡(θ2)−cos⁡(θ2))

In particular, BSH(π2,0,0,0)=H defines the Hadamard gate. This convention can be obtained with the same definitions of ϕr,ϕt and ϕ0 as in the BS1 convention but by choosing:

θ∈[0,π],r=cos⁡(θ2)et t=sin⁡(θ2)

Quandela Convention Rx:

BSRx(θ∈[0;π],0,0,0)=(cos⁡(θ2)isin⁡(θ2)isin⁡(θ2)cos⁡(θ2))

The operator BSRx is obtained by again choosing θ∈[0,π],r=cos⁡(θ2) and t=sin⁡(θ2), in the previous calculations, along with the selections of ϕt,r,0 given in Table 1.

New parameters for ϕt,r,0 to define BSRx

ϕt=12(ϕad−ϕbc),

ϕr=12(ϕac−ϕbd),

ϕ0=12(ϕad+ϕbc−π)

which can be rewritten

ϕac=ϕ0+ϕr

ϕbd=ϕ0−ϕr

ϕad=ϕ0+ϕt+π2

ϕbc=ϕ0−ϕt+π2

Table 1 provides the correspondence between the old definitions of ϕt,r,0 and the new definitions.

Quandela Convention Ry:

BSRy(θ∈[0;π],0,0,0)=(cos⁡(θ2)−sin⁡(θ2)sin⁡(θ2)cos⁡(θ2))

This corresponds exactly to the RY gate according to IBM’s notation. This expression arises through similar transformations.

2.5 Symmetric Beam Splitter

A beam splitter is said to be symmetric if it is assigned r=12 and t=12. This corresponds to the photon having a “one in two chance of being reflected” and a “one in two chance of being transmitted”.

It means choosing θ=π4 in most conventions

In the first convention,

BS1(θ,ϕr,ϕt,ϕ0)=ei.ϕ0(sin⁡(θ).ei.ϕrcos⁡(θ).e−i.ϕtcos⁡(θ).ei.ϕt−sin⁡(θ).e−i.ϕr)

From a mathematical point of view, we choose ϕr=ϕt=ϕ0=0 with θ=π4 and we obtain, for example:

BS1(π4,0,0,0)=1(12.112.112.1−12.1)=12(111−1)=H

The values of the symmetric beam splitter in the Quandela convention (for θ=π2 and all phases equal to 0) are:

BSH(π2,0,0,0,0)=12(111−1)

which represents the Hadamard gate

BSRx(π2,0,0,0,0)=12(1ii1)

BSRy(π2,0,0,0,0)=12(1−111)

All these conventions are equivalent in the sense that one can be transformed into another.

Now, let us clarify a subtlety that was previously overlooked. The Perceval implementation is as follows:

BSRx(θ,Φtl,Φbl,Φtr,Φbr)=(ei.(Φtl+Φtr)cos⁡(θ2)i.ei.(Φbl+Φtr)sin⁡(θ2)i.ei.(Φbr+Φtl)sin⁡(θ2)ei.(Φbl+Φbr)cos⁡(θ2))

BSH(θ,Φtl,Φbl,Φtr,Φbr)=(ei.(Φtl+Φtr)cos⁡(θ2)ei.(Φbl+Φtr)sin⁡(θ2)ei.(Φbr+Φtl)sin⁡(θ2)−ei.(Φbl+Φbr)cos⁡(θ2))

BSRy(θ,Φtl,Φbl,Φtr,Φbr)=(ei.(Φti+Φtr)cos⁡(θ2)−ei.(Φbl+Φtr)sin⁡(θ2)ei.(Φbr+Φtl)sin⁡(θ2)ei.(Φbl+Φbr)cos⁡(θ2))

With the special cases:

BSH(θ,0,0,0,0)=BSH(θ)=12(cos⁡(θ2)sin⁡(θ2)sin⁡(θ2)−cos⁡(θ2))(13)

BSRx(θ,0,0,0,0)=BSRx(θ)=12(cos⁡(θ2)i.sin⁡(θ2)i.sin⁡(θ2)cos⁡(θ2))(14)

BSRy(θ,0,0,0,0)=BSRy(θ)=12(cos⁡(θ2)−sin⁡(θ2)sin⁡(θ2)cos⁡(θ2))(15)

It should be noted that the library allows passing 4 phases as parameters, which differs from the 3 phases described previously. Fig. 5 shows on the left the placement of these four phases the library uses, to be compared with the three parameters of the generic beam splitter on the right we described.

Figure 5: Differences between the phase parameters of the beam splitter provided by the Perceval library [27] and those considered previously.

The phases considered are those associated with each complex coefficient describing the action of the beam splitter, whereas Perceval adds a phase to each mode (inputs and outputs) that can be controlled using phase shifters (described in the next section). This allows them to obtain a more tunable beam splitter. By default, in the Quandela library, the parameters are θ=π2 and ϕ=0 for all ϕ, so that:

BSH(θ)=12(cos⁡(θ2)sin⁡(θ2)sin⁡(θ2)−cos⁡(θ2))

BSH()=BSH(π2,0,0,0,0)=12(111−1)

The notation BSH() corresponds to calling the procedure BSH() without specifying the parameters.

3  Photonic Circuit Representation and Two Other Optical Components

3.1 The Phase Shifter Component

Apart from the beam splitter, there is a second essential optical component that enables transformations on Fock states. One can also use phase shifters: a transparent plate acting on a single input mode and allowing the phase of the output mode to be changed if a photon is present. Physically, a phase shifter is a device that “slows down” the photon, creating a phase shift relative to the rest of the system. There are several ways to implement it. One way is simply to place a transparent plate in the photon’s path, as stressed in the figure below (Fig. 6).

Figure 6: Representation of a type of physical phase shifter.

The phase shifter acts on a single input mode and produces a single output mode. In other words, its action can be described in terms of Fock states. Recall that the basis states of qubits are encoded on two modes by

|0⟩=|10⟩01

And

|1⟩=|01⟩01

Let PSp(ϕ) denote the operator representing the application of a phase shifter with phase ϕ to mode p. Its action on Fock states can be described as follows:

PSp(ϕ)|0⟩p↦|0⟩p and PSp(ϕ)|1⟩p↦eiϕ|1⟩p

Its action on Fock states can be described as follows:

PS1(ϕ)|0⟩=PS1(|10⟩01)=|10⟩01=|0⟩

This time, applying it to the qubit |1⟩:

PS1(ϕ)|1⟩=PS1(ϕ)|01⟩01=eiϕ|01⟩01=eiϕ|1⟩

We thus have the action on qubits:

PS0(ϕ)=(ei.ϕ001) and PS1(ϕ)=(100eiϕ)

Reminders

According to the definition of Rz:

Rz(ϕ)=e−iZθ2=(e−iϕ200eiϕ2)