Advanced time-domain photonic circuits

Authors: Fabian Laudenbach and Nicolas Quesada

In a previous tutorial we introduced the basics of time-domain photonic circuits — the same technology underpinning Xanadu’s quantum computational advantage hardware </demos/tutorial_borealis_beginner>, Borealis [1] — by translating a simple standard Strawberry Fields program into a strawberryfields.TDMProgram. Note that the examples considered previously are relatively simple, involving at most one spatial mode and two concurrent modes. In this tutorial we will go into more depth and show how one can take typical diagrams of experimental setups and from them write their correct representation as time-domain programs.

Note

We encourage you to read the introductory time-domain photonic circuits tutorial before jumping into this one!

To achieve the goal of this tutorial we will propose a set of general rules to translate diagrams of experimental setups into time-domain Strawberry Fields programs. We will apply these rules to two concrete experimental setups and then, after gaining some intuition, will justify their validity and correctness.

Time-domain rules

Now we state the rules to translate a photonic circuit diagram into a time-domain program:

1. Identify the number of spatial modes M in your circuit. This number is equal to the number of wires appearing in the diagram.

2. Obtain the number of concurrent modes in each spatial mode. This is the number of modes that are alive at any given time in each wire. For each spatial mode the number of concurrent modes is one more than the number of modes living in the delay lines in that wire (if there is any). The number of concurrent modes is specified as a list N where each entry is an integer greater than zero.

3. Once you know the number of concurrent modes n = sum(N), one needs to assign them into the different parts of the circuit. We start with the first spatial mode placing the first qumode q[0] at the end of the wire, where the detector is, and then we update its location by tracing its path from right to left until we can place it on the right most time-bin. If there are more concurrent modes in the first spatial mode because of a delay line we assign modes q[1] to q[delay1] until every concurrent mode is accounted for in the first spatial mode. The indices of the modes must always increase from right to left.

4. We now repeat step 3 for all the spatial modes. Assume we used delay1 modes for the first spatial mode. Then we place qumode q[delay1+1] at the end of the second wire, slide back until the first concurrent mode we meet and then start to assign more qumodes until all the concurrent modes of the second spatial mode (wire) are accounted for. This we continue until all spatial modes have been assigned.

5. Once the qumodes are assigned in a diagram read off which ones enter into the different gates. This step is more easily visualized once the previous steps are completed in a concrete example as we do below.

One dimensional cluster state

As a first example we consider the experiment from Yoshikawa et al. as described in Figure 1 of Ref. [2]. Below we reproduce a simplified version of their diagram:

Note that in this figure we removed some unnecessary details from theirs. We also make explicit that the delay line, the loop in the figure, supports only one time bin. Finally, we have added a rotation gate by \(\pi/2\) so that the squeezed states get entangled in the first beamsplitter. This is necessary to compensate for the slightly different beamsplitter convention we use here.

From the diagram above we can already start to apply the rules introduced above. We immediately conclude that there 2 spatial modes; we also easily see that the first spatial mode supports 1 concurrent mode while the second spatial mode supports 2 concurrent modes. This is because there is one concurrent mode living in the delay line. Thus we can already write

N = [1, 2]

This simple assignment completes steps 1 and 2. Our next step is to assign the labels for the first spatial mode. We first put q[0] at the end of the first spatial mode:

and then we slide them backwards until we hit a concurrent mode, represented by the red “pulse”

This completes step 3. For step 4 we repeat this process obtaining

With this final diagram in mind we can now move to step 5 and write the circuit.

import numpy as np
import strawberryfields as sf

sq_r = 4  # Set the value of the squeezing parameter
l = 100  # We will make a 100 temporal mode cluster state

# Measurement angles for the two spatial modes.
# Note the first half of the measurement are in the X quadrature
# and the second half is in the P quadrature
# This will be important when we process the samples below!
theta1 = [0] * int(l / 2) + [np.pi / 2] * int(l / 2)
theta2 = theta1

shots = 1
prog = sf.TDMProgram(N=[1, 2])

with prog.context(theta1, theta2) as (p, q):
    sf.ops.Sgate(sq_r) | q[0]
    sf.ops.Sgate(sq_r) | q[2]
    sf.ops.Rgate(np.pi / 2) | q[0]
    sf.ops.BSgate() | (q[0], q[2])
    sf.ops.BSgate() | (q[0], q[1])
    sf.ops.MeasureHomodyne(p[0]) | q[0]
    sf.ops.MeasureHomodyne(p[1]) | q[1]

eng = sf.Engine("gaussian")

result = eng.run(prog, shots=shots)
samples = result.samples
print("The shape of our samples is", samples.shape)

Out:

The shape of our samples is (1, 2, 100)

With these samples we can verify that we have generated an entangled state by computing the variances of the nullifiers of certain linear combinations of the quadrature operators. These nullifiers are given by (cf. Eq. 1 of Ref. [2])

\[\begin{split}X_k = x_{A,k}+x_{B,k}+x_{A,k+1} - x_{B,k+1}, \\ P_k = p_{A,k}+p_{B,k}-p_{A,k+1} - p_{B,k+1},\end{split}\]

where \(A\) and \(B\) refer to the two spatial modes and \(k\) is an index labeling the different temporal modes.

X_A = samples[0][0][: l // 2]  # X samples from first spatial mode
P_A = samples[0][0][l // 2 :]  # P samples from first spatial mode
X_B = samples[0][1][: l // 2]  # X samples from second spatial mode
P_B = samples[0][1][l // 2 :]  # P samples from second spatial mode

Having split the samples we easily calculate the nullifiers

ntot = len(X_A) - 1
nX = np.array([X_A[i] + X_B[i] + X_A[i + 1] - X_B[i + 1] for i in range(ntot)])
nP = np.array([P_A[i] + P_B[i] - P_A[i + 1] + P_B[i + 1] for i in range(ntot)])

nXvar = np.var(nX)
nPvar = np.var(nP)
print("The variance of the X nullifier is", nXvar)
print("The variance of the P nullifier is", nPvar)

Out:

The variance of the X nullifier is 0.0014269252205099667
The variance of the P nullifier is 0.0010979640031097433

We can compare these values with expected theoretical values

print(
    "The expected value of the variance of the nullifiers is",
    2 * sf.hbar * np.exp(-2 * sq_r),
)

Out:

The expected value of the variance of the nullifiers is 0.0013418505116100474

Two dimensional cluster state

As a second example we consider the experiment from Asavanant et al. as described in Figure 1.c of Ref. [3]. Below we reproduce a simplified version of their diagram:

As with the previous example we have simplified the diagram, added red pulses to make explicit how many concurrent modes are in each spatial mode and added two extra rotation gates to compensate for the way we define our beamsplitter. We have also assumed that the second longer loop supports 5 concurrent modes.

We easily identify the number of current modes in each spatial mode. The first and last spatial modes have no delay lines so only have one concurrent mode. The second and third spatial modes have delay lines supporting 1 and 5 concurrent modes from which we easily deduce that

N = [1, 2, 6, 1]

Now we start assigning qumodes to the concurrent modes. We start with the first spatial mode for which we put q[0] at the end of the wire and slide it back until we hit a concurrent mode to obtain

Now we move to the second spatial mode. For which we assign as follows

Next we look at the third spatial mode for which we need to assign 6 qumodes. The first one is easy

Then we consecutively assign the rest of the modes to the loop

And finally we assign the sixth qumode of this spatial mode

We are left with the last qumode which is trivial to assign and we obtain a fully labeled circuit

With this labeled circuit we can easily write the TDMProgram

sq_r = 5
n = 400  # number of timebins
theta_A = [0] * int(n / 2) + [np.pi / 2] * int(
    n / 2
)  # measurement angles for detector A
theta_B = theta_A  # measurement angles for detector B
theta_C = theta_A
theta_D = theta_A
shots = 1

prog = sf.TDMProgram(N)
with prog.context(theta_A, theta_B, theta_C, theta_D) as (p, q):

    sf.ops.Sgate(sq_r) | q[0]
    sf.ops.Sgate(sq_r) | q[2]
    sf.ops.Sgate(sq_r) | q[8]
    sf.ops.Sgate(sq_r) | q[9]

    sf.ops.Rgate(np.pi / 2) | q[0]
    sf.ops.Rgate(np.pi / 2) | q[8]

    sf.ops.BSgate() | (q[0], q[2])
    sf.ops.BSgate() | (q[8], q[9])
    sf.ops.BSgate() | (q[2], q[8])
    sf.ops.BSgate() | (q[0], q[1])
    sf.ops.BSgate() | (q[3], q[9])

    sf.ops.MeasureHomodyne(p[0]) | q[0]
    sf.ops.MeasureHomodyne(p[1]) | q[1]
    sf.ops.MeasureHomodyne(p[2]) | q[3]
    sf.ops.MeasureHomodyne(p[3]) | q[9]

eng = sf.Engine("gaussian")
result = eng.run(prog, shots=shots)
samples = result.samples

To verify the correctness of the simulation we can again look at the variances of the nullifiers. For this two-dimensional cluster the nullifiers are given by

\[\begin{split}X_k = x_{A,k} + x_{B,k} - \frac{1}{\sqrt{2}}\left(-x_{A,k+1}+x_{B,k+1}+x_{C,k+T}+x_{D,k+T} \right),\\ P_k = p_{A,k} + p_{B,k} + \frac{1}{\sqrt{2}}\left(-p_{A,k+1}+p_{B,k+1}+p_{C,k+T}+p_{D,k+T} \right),\\ Y_k = x_{C,k} - x_{D,k} - \frac{1}{\sqrt{2}}\left(-x_{A,k+1}+x_{B,k+1}-x_{C,k+T}-x_{D,k+T}\right),\\ Q_k = p_{C,k} - p_{D,k} + \frac{1}{\sqrt{2}}\left(-p_{A,k+1}+p_{B,k+1}-p_{C,k+T}-p_{D,k+T}\right),\end{split}\]

where we used \(A,B,C,D\) to label the four spatial modes and the index \(k\) to label temporal modes. The variances of the nullifiers should be given by \(2 \hbar e^{-2 r}\) where \(r\) is the squeezing parameter of the input squeezed states. With this information we numerically compute the nullifiers

X_A = samples[0][0][: n // 2]  # X samples from detector A
P_A = samples[0][0][n // 2 :]  # P samples from detector A
X_B = samples[0][1][: n // 2]  # X samples from detector B
P_B = samples[0][1][n // 2 :]  # P samples from detector B
X_C = samples[0][2][: n // 2]  # X samples from detector C
P_C = samples[0][2][n // 2 :]  # P samples from detector C
X_D = samples[0][3][: n // 2]  # X samples from detector D
P_D = samples[0][3][n // 2 :]  # P samples from detector D

T = 5  # Delay in third spatial mode
c = np.sqrt(1 / 2)
ntot = len(X_A) - T - 1
X = np.array(
    [
        X_A[k] + X_B[k] - c * (-X_A[k + 1] + X_B[k + 1] + X_C[k + T] + X_D[k + T])
        for k in range(ntot)
    ]
)
Y = np.array(
    [
        X_C[k] - X_D[k] - c * (-X_A[k + 1] + X_B[k + 1] - X_C[k + T] - X_D[k + T])
        for k in range(ntot)
    ]
)
P = np.array(
    [
        P_A[k] + P_B[k] + c * (-P_A[k + 1] + P_B[k + 1] + P_C[k + T] + P_D[k + T])
        for k in range(ntot)
    ]
)
Q = np.array(
    [
        P_C[k] - P_D[k] + c * (-P_A[k + 1] + P_B[k + 1] - P_C[k + T] - P_D[k + T])
        for k in range(ntot)
    ]
)

Xvar = np.var(X)
Yvar = np.var(Y)
Pvar = np.var(P)
Qvar = np.var(Q)

print("The variance of the X nullifier is", Xvar)
print("The variance of the Y nullifier is", Yvar)
print("The variance of the P nullifier is", Pvar)
print("The variance of the Q nullifier is", Qvar)
print(
    "The expected value of the variances is", 2 * sf.hbar * np.exp(-2 * sq_r),
)

Out:

The variance of the X nullifier is 0.0001889837390515925
The variance of the Y nullifier is 0.0001878866858368039
The variance of the P nullifier is 0.0001895953820689966
The variance of the Q nullifier is 0.0001836553442700422
The expected value of the variances is 0.00018159971904993942

We obtain an excellent agreement between our numerics and the theory.

Why does it work?

We conclude our tutorial by arguing why the rules we introduced earlier work. The first two rules should be fairly straightforward to understand: to simulate time-domain circuits we need to correctly account for how many spatial modes there are and also count correctly how many modes are alive in each wire at a given time. The third and fourth rules are a bit less obvious. Why do we insist that the qumodes with the smallest index are the ones closest to detectors? The reason is simple: it is because in our description of TDMPrograms we always shift the modes “downwards”. Recall from the basic time-domain tutorial that every strawberryfields.TDMProgram is equivalent to regular Strawberry Fields program where we apply a strawberryfields.shift_by() operation. This is why we always want to respect the chronology of the measurements. This is satisfied if and only if for every wire the qumodes are assigned consecutively and with indices decreasing from left to right.

References

1

Quantum computational advantage with a programmable photonic processor.

2(1,2)

J. Yoshikawa, S. Yokoyama, T. Kaji, C. Sornphiphatphong, Y. Shiozawa, K. Makino and A. Furusawa. Generation of one-million-mode continuous-variable cluster state by unlimited time-domain multiplexing. APL Photonics, 2016. doi:10.1063/1.4962732 .

3

W.Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J. Yoshikawa, N. C. Menicucci, H. Yonezawa and A. Furusawa. Generation of time-domain-multiplexed two-dimensional cluster state. Science, 2019. doi:10.1126/science.aay2645 .

demos/run_advanced_time_domain

 

Download Python script

 

Download Notebook