Time-domain photonic circuits

Authors: Fabian Laudenbach, Josh Izaac, Theodor Isacsson and Nicolas Quesada

In this tutorial we introduce the basic ideas of time-domain photonic systems together with the time-domain module of StrawberryFields, strawberryfields.tdm. Time-domain multiplexing allows for the creation of massive quantum systems having millions of entangled modes as shown in a number of recent experimental demonstrations [1], [2], [3], [4]. To motivate this architecture we follow Takeda et al. [1] and study the one loop setup shown in the picture below


An instruction set for the time-domain architecture

To introduce the time domain architecture we will simulate the setup in the figure. We can summarize the setup as the following set of intructions:

  1. The i+1 mode is squeezed by a certain amount r.

  2. The i+1 mode undergoes a beamsplitter transformation with its preceding mode, labeled by i, by a certain angle alpha[i]. This beamsplitter has the effect of switching in and out of the loop the modes i+1 and i respectively.

  3. The i+1 mode has entered the loop after the beamsplitter transformation and is rotated by angle phi[i].

  4. After interacting with the i+1 mode, the mode in the previous time bin i has left the loop and is homodyned by measurement angle theta[i].

We can use n modes to represent n timebins and implement the setup described by the figure using an standard StrawberryFields program.

import strawberryfields as sf
from strawberryfields.ops import Sgate, BSgate, Rgate, MeasureHomodyne
import numpy as np

Below we assume that the parameters of the squeezing gate \((S)\) and the angles of the beamsplitter \((BS)\) and rotation \(R\) do not change with the timebin.

np.random.seed(42) # We set the seed to facilitate later comparison
n = 20

r = 1.0
length =  num_modes - 1
alpha = [np.pi/4]*length
phi = [0]*length
theta = [0]*length

prog0 = sf.Program(n)
with prog0.context as q:
    for i in range(num_modes-1):
eng0 = sf.Engine('gaussian')
result0 = eng0.run(prog0)

# Every mode except the last two is measured.

In this program we follow each an every time-bin mode (a total of n=20 modes) as they progress through the different operations. However, at most only two modes are interacting at the same time. The rest of the modes either have not been squeezed (and thus are in the vacuum state) or have been measured (and thus reset to vacuum).

Instead of keeping track of the modes using labels telling us their time bin, we can label them in terms of the position that they occupy in the circuit above. A convenient choice of labels is to look at the two modes that enter into the beamsplitter. We can use 0 for the mode that enters in the top port of the beamsplitter and which has been circulating in the loop. We can use the label 1 for the mode that has just been squeezed and is about to enter the beamplitter in the bottom port. These two modes, which have fixed positions in the circuit, we call concurrent modes; they are, at all times, the only two modes that are not in vacuum in the circuit.

To make a closed loop of modes going in and out of the circuit one can recycle the mode that is measured (mode 0) and plug it back-in as a vacuum mode re-entering the circuit. As it turns out the Gaussian backend of Strawberry Fields automatically resets to vacuum a mode that has been measured thus the resetting is automatically taken care of.

To implement the recycling we will need a function that takes the mode that is measured and is at the beginning of the mode list q and that after measurement moves it back to the end. This functionality is provided by the function shift_by in the module strawberryfield.tdm.tdmprogram as illustrated below:

from strawberryfields.tdm.tdmprogram import shift_by
l = [0, 1, 2, 3, 4]
shift_by(l, 1)

We are now ready to simulate an n=20 mode circuit using only two concurrent modes:

np.random.seed(42) # We set the seed to facilitate comparison
prog1 = sf.Program(2)
num_modes = 20
with prog1.context as q:
    for i in range(num_modes-1):
        Sgate(r) | q[1]
        BSgate(alpha[i]) | (q[0], q[1])
        Rgate(phi[i]) | q[1]
        MeasureHomodyne(theta[i]) | q[0] # We always measure the register in the 0th position in q
        q = shift_by(q, 1)
eng1 = sf.Engine('gaussian')
result1 = eng1.run(prog1)

Eventhough we measured n-1 times we only get 2 values. This is because the samples attribute only contains the last ever measured value(s) in our modes; since we only have two concurrent modes we only get two numbers. To obtain the complete sample record we use the attribute all_samples


Note that we obtain the same samples as in result0.samples except that they are ordered in a slightly different arrangement. We can put them in the correct order by using reshape_samples

from strawberryfields.tdm.tdmprogram import reshape_samples
reshape_samples(result1.all_samples, [0,1], [2])
# The second argument [0,1] gives the labels of the concurrent modes were measurement happend,
# while the third argument [2] indicates, for this simple one loop program, that there are only two concurrent modes in a single loop.

Where now all the measurement outcomes are attached to the only mode that we ever measure, namely 0. Note that in result1.all_samples the samples are attached to different modes, this is an artifact of the recycling.

By recycling the modes we are now able to simulate very long circuits with just two modes. This motivates the introduction of the TDMProgram class. In this class one needs to specify the number of active modes (N=2 for the example above) and then the class takes care of automatically shifting the modes and reshaping the samples. Thus the program above could have been rewritten as

from strawberryfields.tdm import tdmprogram

np.random.seed(42) # We set the seed to facilitate comparison
N = 2 # Number of concurrent modes
prog2 = tdmprogram.TDMProgram(N=N)

with prog2.context(alpha, phi, theta) as (p, q):
    Sgate(r, 0) | q[1]
    BSgate(p[0]) | (q[0], q[1])
    Rgate(p[1]) | q[1]
    MeasureHomodyne(p[2]) | q[0]
eng2 = sf.Engine('gaussian')
result2 = eng2.run(prog2)

In a TDMProgram we need to specify the number of “concurrent” modes N, meaning the total number of modes that are “alive” at a given the time in the circuit. Then we pass to the context the sequences of the different parameters that will be used in the program, in our case the lists alpha, phi and theta. The three variables in the list p will cycle through the corresponding values of these lists, for example, p[0] will go through the values of alpha as the simulation progresses.

A single-loop architecture

We will now write a simple function that encapsulates the single loop architecture using a TDMprogram. Then we will use it to generate samples of interesting quantum states and investigate their correlations.

def singleloop(r, alpha, phi, theta, copies):
    """Single loop program.

        r (float): squeezing parameter
        alpha (Sequence[float]): beamsplitter angles
        phi (Sequence[float]): rotation angles
        theta (Sequence[float]): homodyne measurement angles
        copies (int)
        (list): homodyne samples from the single loop simulation

    prog = tdmprogram.TDMProgram(N=2)
    with prog.context(alpha, phi, theta, copies=copies) as (p, q):
        Sgate(r, 0) | q[1]
        BSgate(p[0]) | (q[0], q[1])
        Rgate(p[1]) | q[1]
        MeasureHomodyne(p[2]) | q[0]
    eng = sf.Engine("gaussian")
    result = eng.run(prog)
    # We only want the samples from mode 0 since it is the only one we ever measure.
    return result.samples[0]

We can use the function we just developed to reproduce for the fourth time the results of our very simple experiment:

np.random.seed(42) # We set the seed to facilitate comparison
singleloop(r, alpha, phi, theta, 1) # The last argument, 1, is used for the number of samples we want

Generation of EPR states

With the single loop architecture we can easily generate Einstein-Podolsky-Rosen (EPR) states [5]. In this section we will create this bipartite state and probe it using homodyne measurements.

r = 2.0
alpha = [np.pi/4, 0]
phi = [0, np.pi/2]
theta = [0,0] # In this setting we will loop at the positions of the two oscillators, thus we se the homodyne angles accordingly
copies = 2500

samplesEPRxx = singleloop(r, alpha, phi, theta, copies)

Recall that the samples are lumped together i.e. they have the format \(\text{[mode 0 first run, mode 1 first run, mode 0 second run, mode 1 second run, ...]}\) We will transform this into the ordering \(\text{[[mode 0 first run, mode 0 second run, ...], [mode 1 first run, mode 1 second run, ...]}\) by using the reshaping method and a transposition

samplesxx = samplesEPRxx.reshape([-1,2]).T
x0 = samplesxx[0]
x1 = samplesxx[1]

# We now want to look at the samples
import matplotlib.pyplot as plt

plt.figure(figsize = (5,5))
plt.plot(x0, x1,".")

We can put this observation in more quantitative terms by looking at the variance of \(x_1-x_0\):

sample_var_xx = np.var(x0-x1)

As it turns out the variance of the operator \(N_1 = x_1 - x_0\) is related to the amount of squeezing used to prepare the state. One can easily find that \(\text{Var}(N) = 2 e^{-2 r}\).

expected_var_xx = 2 * np.exp(-2*r)

In the literature, the operator \(N\) for which one finds a relation like the one above, is called a Nullifier of the state. A way to understand this name is that the variance of this operator approaches zero –nullifies– as the squeezing in the state grows [1].

The EPR state has a second nullifier, \(N_2 = p_1+p_2\). We can confirm this by running our circuit again, but this time measuring in the \(P\) quadrature

theta = [np.pi/2, np.pi/2] # Now we homodyne the p-quadratures by setting thr angle to pi/2
samplesEPRpp = singleloop(r, alpha, phi, theta, copies)

samplespp = samplesEPRpp.reshape([-1,2]).T # As before we reshape the samples
p0 = samplespp[0]
p1 = samplespp[1]

plt.figure(figsize = (5,5))
plt.plot(p0, p1,".")

# We now calculate the sample variance and the expected variance,
# except for the effect of finite size statistics, they are almost identical
np.var(p0+p1), 2*np.exp(-2*r)

Generating GHZ states

The GHZ states (named after Greenberger, Horne and Zeilinger [[#greenberger2007]_] can be thought as an \(n\)-partite generalization of the EPR states [7]. We will prepare these states using the single loop architecture encapsulated in the function singleloop defined before.

n = 10
vac_modes = 1 # This is an ancilla mode that is used at the beginning of the protocol
copies = 1
r = 4
alpha = [np.arccos(np.sqrt(1 / (n - i + 1))) if i != n + 1 else 0 for i in range(n + vac_modes)]
alpha[0] = 0.0
phi = [0] * (n + vac_modes)
phi[0] = np.pi / 2
theta = [0] * (n + vac_modes) # We will measure first all the states in the X quadrature

singleloop(r, alpha, phi, theta, copies)

Note that except for the first mode all the sampled positions are the same. We can also sample our state in the \(p\) basis

theta = [np.pi/2] * (n + vac_modes) # Now we measure in p quadrature

samplep = singleloop(r, alpha, phi, theta, copies)

There does not seem to be anything interesting going on when we look at the \(p\) quadratures, yet if we consider the sum of the sampled values of the momenta we find that the momenta add up to approximately zero:

np.sum(samplep[1:]) # Note that we exclude the first element, the "vacuum mode"

Indeed, like for the EPR state defined before one can introduce nullifiers for the GHZ states, these are given by

We now generate many samples and statistically verify that the variances of the nullifier decay exponentially with the amount of squeezing

# First collect x-samples
copies = 1000
theta = [0] * (n + vac_modes)
samplesGHZx = singleloop(r, alpha, phi, theta, copies)
# Then collect p-samples
theta = [np.pi/2] * (n + vac_modes)
samplesGHZp = singleloop(r, alpha, phi, theta, copies)

# As before we will reshape the samples so that they are in a rectangular array
reshaped_samplesGHZx = samplesGHZx.reshape([copies, n + vac_modes])
reshaped_samplesGHZp = samplesGHZp.reshape([copies, n + vac_modes])

nullifier_X = lambda sample: (sample - sample[-1])[vac_modes:-1]
val_nullifier_X = np.var([nullifier_X(x) for x in reshaped_samplesGHZx], axis=0)

# The expected value of this variance is

nullifier_P = lambda sample: np.sum(sample[vac_modes:])
val_nullifier_P = np.var([nullifier_P(p) for p in reshaped_samplesGHZp], axis=0)




S. Takeda, T. Kan and A. Furusawa. On-demand photonic entanglement synthesizer. Science Advances, 2019. doi:10.1126/sciadv.aaw4530 .


J. Yoshikawa, S. Yokoyama, T. Kaji, C. Sornphiphatphong, Y. Shiozawa, K. Makino1 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 .


M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen and U. L. Andersen. Deterministic generation of a two-dimensional cluster state. Science, 2019. doi:10.1126/science.aay4354 .


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 .


A. Einstein, B. Podolsky and N. Rosen. Can quantum-mechanical description of physical reality be considered complete? Physical Review, 1935. doi:10.1103/PhysRev.47.777 .


D. M. Greenberger, M. A. Horne and A. Zeilinger. Going beyond Bell’s Theorem arXiv:0712.0921, 2007.


P. van Loock, S. L. Braunstein. Multipartite Entanglement for Continuous Variables: A Quantum Teleportation Network Physical Review Letters, 2000. doi:10.1103/PhysRevLett.84.3482 .

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery