Implementation of a network of interconnected beamsplitters for coherent state evolution

Photonic Software Mr Mustard
1

Hi good morning dear colleagues, I am trying to evolve a system with six coherent states as input to three beamsplitters and I am facing the following errors described in the images below. How would I make this system evolve correctly?

# import numpy as np
from mrmustard.lab import *

vac = Vacuum(num_modes=2)        # 2-mode vacuum state
coh0 = Coherent(x=-0.1, y=-0.4)    # coh state |-alpha> with alpha = -0.1 - 0.4j
coh1 = Coherent(x=0.1, y=-0.4)    # coh state |alpha> with alpha = 0.1 - 0.4j
sq  = SqueezedVacuum(r=0.5)      # squeezed vacuum state
g   = Gaussian(num_modes=2)      # 2-mode Gaussian state with zero means
fock4 = Fock(4)                  # fock state |4>

D  = Dgate(x=1.0, y=-0.4)         # Displacement by 1.0 along x and -0.4 along y
S  = Sgate(r=0.5)                 # Squeezer with r=0.5

BS = BSgate(theta=np.pi/4)          # 50/50 beam splitter
S2 = S2gate(r=0.5)                  # two-mode squeezer
MZ = MZgate(phi_a=0.3, phi_b=0.1)   # Mach-Zehnder interferometer
I  = Interferometer(8)              # 8-mode interferometer
L  = Attenuator(0.5)                # pure lossy channel with 50% transmissivity
A  = Amplifier(gain=2.0, nbar=1.0)  # noisy amplifier with 200% gain

## Primeira bateria de BS

Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0,1]
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2,3]
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[4,5]

## Segunda bateria de BS
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[3,5]
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[1,6]
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2,4]

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

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 4
      1 ## Primeira bateria de BS
      3 Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0,1]
----> 4 Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2,3]
      5 Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[4,5]

File ~\anaconda3\lib\site-packages\mrmustard\lab\abstract\state.py:500, in State.__rshift__(self, other)
    496 if issubclass(other.__class__, State):
    497     raise TypeError(
    498         f"Cannot apply {other.__class__.__qualname__} to a state. Are you looking for the << operator?"
    499     )
--> 500 return other.primal(self)

File ~\anaconda3\lib\site-packages\mrmustard\lab\abstract\transformation.py:56, in Transformation.primal(self, state)
     47 r"""Applies ``self`` (a ``Transformation``) to other (a ``State``) and returns the transformed state.
     48 
     49 Args:
   (...)
     53     State: the transformed state
     54 """
     55 if state.is_gaussian:
---> 56     new_state = self.transform_gaussian(state, dual=False)
     57 else:
     58     new_state = self.transform_fock(state, dual=False)

File ~\anaconda3\lib\site-packages\mrmustard\lab\abstract\transformation.py:103, in Transformation.transform_gaussian(self, state, dual)
     93 r"""Transforms a Gaussian state into a Gaussian state.
     94 
     95 Args:
   (...)
    100     State: the transformed state
    101 """
    102 X, Y, d = self.XYd if not dual else self.XYd_dual
--> 103 cov, means = gaussian.CPTP(state.cov, state.means, X, Y, d, state.modes, self.modes)
    104 new_state = State(
    105     cov=cov, means=means, modes=state.modes, _norm=state.norm
    106 )  # NOTE: assumes modes don't change
    107 return new_state

File ~\anaconda3\lib\site-packages\mrmustard\physics\gaussian.py:432, in CPTP(cov, means, X, Y, d, state_modes, transf_modes)
    411 r"""Returns the cov matrix and means vector of a state after undergoing a CPTP channel.
    412 
    413 Computed as ``cov = X \cdot cov \cdot X^T + Y`` and ``d = X \cdot means + d``.
   (...)
    429     Tuple[Matrix, Vector]: the covariance matrix and the means vector of the state after the CPTP channel
    430 """
    431 if not set(transf_modes).issubset(state_modes):
--> 432     raise ValueError(
    433         f"The channel should act on a subset of the state modes ({transf_modes} is not a subset of {state_modes})"
    434     )
    435 # if single-mode channel, apply to all modes indicated in `modes`
    436 if X is not None and X.shape[-1] == 2:

ValueError: The channel should act on a subset of the state modes ([2, 3] is not a subset of [0, 1])

conda 23.3.1
mrmustard 0.3.0
numpy 1.23.5

2

Hey @Antonio_Aguiar! Welcome to the forum :rocket:

The issue here is that there are only two states indexed by 0 and 1. New state modes aren’t created after the first line you have :slightly_smiling_face:. You can fix this by doing BS[0, 1] instead of BS[2, 3]. This actually applies to every line after that, as well! Here’s a demonstration:

states = [
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2, 3]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[4, 5]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[3, 5]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[1, 6]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2, 4]",
]

for state in states:
    try:
        eval(state)
    except ValueError as e:
        print(state)
        print(e)
        print()

'''
Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2, 3]
The channel should act on a subset of the state modes ([2, 3] is not a subset of [0, 1])

Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[4, 5]
The channel should act on a subset of the state modes ([4, 5] is not a subset of [0, 1])

Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[3, 5]
The channel should act on a subset of the state modes ([3, 5] is not a subset of [0, 1])

Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[1, 6]
The channel should act on a subset of the state modes ([1, 6] is not a subset of [0, 1])

Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[2, 4]
The channel should act on a subset of the state modes ([2, 4] is not a subset of [0, 1])
'''

So, just change the indices you access in each BS to [0,1]:

states = [
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
    "Coherent(x=[1.0, 2.0], y=[-1.0, -2.0]) >> BS[0, 1]",
]
3

If your goal was to create a multi-mode state, then you need to compound those Coherent operations into one object like this:

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

print(state.num_modes) # 6

state >> BS[0, 1] >> BS[2, 3] >> BS[4, 5] # this works!

You can find out more about this here: Basic API Reference — Mr Mustard 0.4.1 documentation

1 Like
4

Thank you for your attention. The multimode solution suits me best. About the last code, when you defined the states, there are six coherent states that will then pass through the beam splitters, correct? My next step is to build a function that switches the inputs between |alpha> and |-alpha> so that I evaluate all 64 possible inputs and evolution of my system. Do you know a way to do this analysis?

5

Yep! In my code example, state is a 6-mode coherent state that then goes through three beam splitters that act on specific modes.

My next step is to build a function that switches the inputs between |alpha> and |-alpha> so that I evaluate all 64 possible inputs and evolution of my system.

Interesting! Well, another way that you can make a Coherent state is by applying the displacement operator to a vacuum state. You can do this iteratively as follows:

alpha = [0.1, 0.2] # alpha[0] = x, alpha[1] = y
num_modes = 4

state = Vacuum(num_modes=num_modes)

for m in range(num_modes):
    state >> Dgate(x=alpha[0], y=alpha[1])

Hope this helps!