Beam Splitter (BS)

Hey, good morning.
I’m having difficulty implementing a beam splitter whose inputs receive alpha and -alpha coherent states, for example. How do I set the input state (alpha or -alpha) to interfere with the BS? Is it possible to place a phase shifter on one of the BS inputs? How do I measure output states?

import numpy as np
from mrmustard.lab import *

state = (Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]))

print(state.num_modes) # 2

state >> BSgate(theta=np.pi/4)[0, 1]

results = Gaussian(2) << PNRDetector(efficiency = 0.9, modes = [0])
results[0]
results[1]

If you want help with diagnosing an error, please put the full error message below:

C:\Users\anton\MrMustard\mrmustard\math\backend_numpy.py:117: RuntimeWarning: invalid value encountered in cast
  return np.array(array, dtype=dtype)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~\MrMustard\mrmustard\math\backend_manager.py:105, in BackendManager._apply(self, fn, args)
    104 try:
--> 105     attr = getattr(self.backend, fn)
    106 except AttributeError:

AttributeError: 'BackendNumpy' object has no attribute 'convolution'

During handling of the above exception, another exception occurred:

NotImplementedError                       Traceback (most recent call last)
Cell In[3], line 1
----> 1 results = Gaussian(2) << PNRDetector(efficiency = 0.9, modes = [0])
      2 results[0]
      3 results[1]

File ~\MrMustard\mrmustard\lab\detectors.py:97, in PNRDetector.__init__(self, efficiency, dark_counts, efficiency_trainable, dark_counts_trainable, efficiency_bounds, dark_counts_bounds, stochastic_channel, modes, cutoffs)
     90 self._add_parameter(
     91     make_parameter(efficiency_trainable, eff, "efficiency", efficiency_bounds)
     92 )
     93 self._add_parameter(
     94     make_parameter(dark_counts_trainable, dk, "dark_counts", dark_counts_bounds)
     95 )
---> 97 self.recompute_stochastic_channel()

File ~\MrMustard\mrmustard\lab\detectors.py:126, in PNRDetector.recompute_stochastic_channel(self, cutoffs)
    121 dark_prior = math.poisson(max_k=settings.PNR_INTERNAL_CUTOFF, rate=dc)
    122 condprob = math.binomial_conditional_prob(
    123     success_prob=qe, dim_in=settings.PNR_INTERNAL_CUTOFF, dim_out=c
    124 )
    125 self._internal_stochastic_channel.append(
--> 126     math.convolve_probs_1d(
    127         condprob, [dark_prior, math.eye(settings.PNR_INTERNAL_CUTOFF)[0]]
    128     )
    129 )

File ~\MrMustard\mrmustard\math\backend_manager.py:1535, in BackendManager.convolve_probs_1d(self, prob, other_probs)
   1532 for q_ in other_probs[1:]:
   1533     q = q[..., None] * q_[(None,) * q.ndim + (slice(None),)]
-> 1535 return self.convolve_probs(prob, q)

File ~\MrMustard\mrmustard\math\backend_manager.py:1550, in BackendManager.convolve_probs(self, prob, other)
   1548 prob_padded = self.pad(prob, [(s - 1, 0) for s in other.shape])
   1549 other_reversed = other[(slice(None, None, -1),) * other.ndim]
-> 1550 return self.convolution(
   1551     prob_padded[None, ..., None],
   1552     other_reversed[..., None, None],
   1553     data_format="N"
   1554     + ("HD"[: other.ndim - 1])[::-1]
   1555     + "WC",  # TODO: rewrite this to be more readable (do we need it?)
   1556 )[0, ..., 0]

File ~\MrMustard\mrmustard\math\backend_manager.py:449, in BackendManager.convolution(self, array, filters, padding, data_format)
    431 def convolution(
    432     self,
    433     array: Tensor,
   (...)
    436     data_format="NWC",
    437 ) -> Tensor:  # TODO: remove strides and data_format?
    438     r"""Performs a convolution on array with filters.
    439 
    440     Args:
   (...)
    447         The convolved array.
    448     """
--> 449     return self._apply("convolution", (array, filters, padding, data_format))

File ~\MrMustard\mrmustard\math\backend_manager.py:109, in BackendManager._apply(self, fn, args)
    107     msg = f"Function ``{fn}`` not implemented for backend ``{self.backend_name}``."
    108     # pylint: disable=raise-missing-from
--> 109     raise NotImplementedError(msg)
    110 return attr(*args)

NotImplementedError: Function ``convolution`` not implemented for backend ``numpy``.

Anaconda 3
strawberryfields 0.23

Hi @Antonio_Aguiar ,

There seems to be an issue with the Numpy backend. However you can use the Tensorflow backend instead.

By phase shifter do you mean applying a Displacement gate? Or a displacement + squeezing + rotation as shown in the first image here?

You can apply them all by running the code below.

import numpy as np
from mrmustard.lab import *
import mrmustard.math as math
math.change_backend("tensorflow")

# Create a separable two-mode coherent state
state = (Coherent(x=[1.0, -1.0], y=[2.0, -2.0]))

print(state.num_modes) # 2

# Apply 50/50 beam splitter
state = state >> BSgate(theta=np.pi/4)

# Apply squeezing, rotation, and displacement
state = state >> Sgate(r=1, phi=np.pi/2) >> Dgate(y=1.3)

# Perform a measurement
results = state << PNRDetector(efficiency = 0.9, modes = [0])
results[0]  # unnormalized density matrix of mode 1 conditioned on measuring 0 in mode 0
results[1]  # unnormalized density matrix of mode 1 conditioned on measuring 1 in mode 0

I hope this helps!

Firstly, thank you for responding. When I talk about a phase shifter I mean an e^(pi*j) acted on the input state, but if it is clear how I define a coherent state |alpha> and a |-alpha> at the input I won’t need the PS. As I would like to see interference between two coherent states (identical or distinct by a phase) passing through a beam splitter and there will be cases where there will be nothing at one of the outputs or something close to vacuum. I would like to know how I know if a coherent state has arrived or if nothing has arrived at the output. Considering that I am going to enter |alpha> or |-alpha>, what types of analysis should I do to detect the types of interference that occurred in the beam splitter?

Hi @Antonio_Aguiar,

Maybe for the phase shifting you could use a rotation gate (Rgate) if needed.

The PNR and Threshold detectors return an array of unnormalized measurement results, meaning that the elements of the array are the density matrices of the leftover systems, conditioned on the outcomes:

results = Gaussian(2) << PNRDetector(efficiency = 0.9, modes = [0])
results[0]  # unnormalized dm of mode 1 conditioned on measuring 0 in mode 0
results[1]  # unnormalized dm of mode 1 conditioned on measuring 1 in mode 0
results[2]  # unnormalized dm of mode 1 conditioned on measuring 2 in mode 0
# etc...

So for your two-mode state the following code gives you a 50x50 matrix

# Create a separable two-mode coherent state
state = (Coherent(x=[1.0, -1.0], y=[2.0, -2.0]))

# Apply 50/50 beam splitter
state = state >> BSgate(theta=np.pi/4)

# Perform a measurement
result = state << PNRDetector(efficiency = 0.9, modes = [0,1])
result

So results[0] gives you the unnormalized density matrix of mode 1 conditioned on measuring 0 in mode 0.

If needed you can compare this to the result of applying a beamsplitter to the vacuum.

You can find more information on measurements here. Note that Strawberry Fields works differently, but maybe the theory there can help you.

Hi @Antonio_Aguiar,

To complement my previous answer, for the Rgate you can insert it wherever you want it.

e.g. state = Coherent(x=[1.0, -1.0], y=[2.0, -2.0]) >> Rgate(phi=0.4)

Is there a possibility to add photon detectors in Strawberry Fields or Mr Mustard? How would I graph photon statistics?
How would I add a phase to the input in Strawberry Fields simulations?

Hi @Antonio_Aguiar ,

The PNRDetector is indeed a photon detector.

You can use matplotlib or some similar library to graph your results.

In this demo you can learn how to use a rotation gate (so adding a phase) in Strawberry Fields. Note that all of the code I’ve shared so far is for MrMustard, not Strawberry Fields so be careful not to get confused if you’re using both!

1 Like

How would I go about reconstructing the coherent output states of a beam splitter using Strawberry Fields? Would it be possible to simulate this experiment on a real machine?

Hi @Antonio_Aguiar ,

We have one device online but only upon request under limited availability.
I recommend that you take a look at this tutorial to learn what you can and cannot do with it.
You can learn more details about X8 here.

My recommendation would be to do everything in simulation instead. If you really need access to the hardware you can see availability times and the contact email to request access where it says contact us here.

Check out the details here to see the decomposition of a coherent state in the Fock basis.

The Strawberry Fields website and documentation have a lot of information on how to simulate circuits too.

Let me know if you have any further questions.

Hello how are you, thank you for your attention. I have no doubts about the theory of coherent states. My question is about how to reconstruct the output states after passing through a divider. Because I simply didn’t find anything in your (Strawberry Fields) documentation that talks about a topic as basic as this. I would just like to simulate the input of coherent states into a BS (50/50) and analyze the outputs and see if the behavior is similar to what is seen in quantum optics textbooks.

Hi Antonio, in Strawberry Fields in this case you should use a gaussian backend. If you use the Fock backend you won’t be able to recover this information.

Here’s how you can program your coherent states and measure them back as displacements.

# Program without beamsplitter

import strawberryfields as sf
from strawberryfields.ops import *
import numpy as np
import matplotlib.pyplot as plt

# set the random seed
np.random.seed(42)

# Set the info for the coherent states
alpha = 1+0.5j
r = np.abs(alpha)
phi1 = np.angle(alpha)
phi2 = np.angle(-alpha)

# Create the program
prog = sf.Program(2)
eng = sf.Engine("gaussian")

with prog.context as q:
    # prepare initial states
    Coherent(r,phi1) | q[0]
    Coherent(r,phi2) | q[1]

    #BSgate(np.pi/4,0) | (q[0], q[1]) # No beam splitter this time
    MeasureFock() | (q[0], q[1])

# Get and print the results
results1 = eng.run(prog)
print('Photon count: ',results1.samples)
state1 = results1.state
print('State: ', state1)
print('Discplacement for mode 0: ',state1.displacement([0]))
print('Discplacement for mode 1: ',state1.displacement([1]))

Here you can see these in nice plots

# Mode 0
fig = plt.figure()
X = np.linspace(-5, 5, 100)
P = np.linspace(-5, 5, 100)
Z = state1.wigner(0, X, P)
X, P = np.meshgrid(X, P)
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, P, Z, cmap="RdYlGn", lw=0.5, rstride=1, cstride=1)
fig.set_size_inches(4.8, 5)
ax.set_axis_off()

# Mode 1
fig = plt.figure()
X = np.linspace(-5, 5, 100)
P = np.linspace(-5, 5, 100)
Z = state1.wigner(1, X, P)
X, P = np.meshgrid(X, P)
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, P, Z, cmap="RdYlGn", lw=0.5, rstride=1, cstride=1)
fig.set_size_inches(4.8, 5)
ax.set_axis_off()

Now add the beam splitter and do the same

# With beam splitter

prog = sf.Program(2)
eng = sf.Engine("gaussian")

with prog.context as q:
    # prepare initial states
    Coherent(r,phi1) | q[0]
    Coherent(r,phi2) | q[1]
    
    BSgate(np.pi/4,0) | (q[0], q[1])
    MeasureFock() | (q[0], q[1])

results2 = eng.run(prog)
print('Photon count: ',results2.samples)
state2 = results2.state
print('State: ', state2)

# Get and print the results
print('Discplacement for mode 0: ',state2.displacement([0]))
print('Discplacement for mode 1: ',state2.displacement([1]))

See the difference

# Mode 0
fig = plt.figure()
X = np.linspace(-5, 5, 100)
P = np.linspace(-5, 5, 100)
Z = state2.wigner(0, X, P)
X, P = np.meshgrid(X, P)
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, P, Z, cmap="RdYlGn", lw=0.5, rstride=1, cstride=1)
fig.set_size_inches(4.8, 5)
ax.set_axis_off()

# Mode 1
fig = plt.figure()
X = np.linspace(-5, 5, 100)
P = np.linspace(-5, 5, 100)
Z = state2.wigner(1, X, P)
X, P = np.meshgrid(X, P)
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, P, Z, cmap="RdYlGn", lw=0.5, rstride=1, cstride=1)
fig.set_size_inches(4.8, 5)
ax.set_axis_off()

I hope this helps you

1 Like

Thank you very much, it helped me a lot. I have one more question, if I choose to now place cat states in the inputs of my BS, do I need to change the backend (due to the fact that it is a Gaussian state)? Would this change be easy to make?

Hi @Antonio_Aguiar,

As you can see in the docs for Catstate “Cat states are non-Gaussian, and thus can only be used in the Fock and Bosonic backends, not the Gaussian backend.”

We have three very good tutorials that go into a lot of detail on how to use the bosonic backend for studying cat states and GKP states. Check them out here, here, and here.

Note that this will not run on our X8 device. Also please note that the photon-counting measurement ( MeasureFock()) is not available on the Bosonic backend.

Hello, Thank you for your attention. Yes I understand the limitations of the Gaussian backend. I had already read the three tutorials because I knew the particularities of non-Gausian states. Is it possible to use “MeasureHD” on the Bosonic backend?

Hi @Antonio_Aguiar , yes you should be able to use MeasureHD on the bosonic backend. Let me know if you have any issues using it.