This example demonstrates how to construct and simulate a photonic integrated circuit (PIC) representing a unitary matrix of arbitrary dimension and values using only directional couplers and phase shifters. A circuit is constructed to model the unitary evolution of a non-interacting Bose-Hubbard model, or 'hopping' Hamiltonian, and the probabilities of success for various measurements are simulated for three-photon input states representing both indistinguishable and distinguishable photons. Lastly, the circuit is simulated with variations in the gap widths of its directional couplers to characterize the impact of manufacturing imperfections on the measured probabilities of success.

## Overview

Understand the simulation workflow and key results

The simulations in this example are performed with the qINTERCONNECT solver. The simulation is set up and run through the Python API, which can either be run directly in INTERCONNECT or in an external Python environment. Since qINTERCONNECT requires an INTERCONNECT GUI license to run, at least two GUI licenses are required to run qINTERCONNECT directly from the INTERCONNECT environment. Documentation for the Python API is available here.

The Clements design is used to deconstruct the unitary matrix into a square mesh of 2x2 beamsplitters and a final array of phase shifters [1]. Each beamsplitter can then be constructed from two phase shifters and two 50:50 directional couplers. As an example, the unitary matrix representing a 'hopping' Hamiltonian is constructed and simulated [2]. As in [2], unitary matrices are generated for 6 different interaction strengths, representing different times. The hopping Hamiltonian can be expressed in matrix form as

0 & t & 0 & t \\

t & 0 & t & 0 \\

0 & t & 0 & t \\

t & 0 & t & 0 \\

\end{pmatrix} $$

from which the corresponding unitary matrix can be found by taking the matrix exponential

$$U_{NN}=e^{-iH_{NN}}.$$

### Step 1: Generate unitary matrix circuits

Six files are generated containing the photonic circuits corresponding to unitary matrices describing a hopping Hamiltonian at timesteps from 0.1 s to 5 s as in [2].

### Step 2: Perform simulations in qINTERCONNECT

The files generated in the previous step are used in a qINTERCONNECT simulation to measure the output of the circuit for indistinguishable and distinguishable photons.

### Step 3: Parameter sweep

The gap widths of the directional couplers inside the unitary matrix circuits are varied from their ideal values to characterize the impact of manufacturing imperfections.

## Run and results

Instructions for running the model and discussion of key results

### Step 1: Generate unitary matrix circuits

- Open the blank simulation file step1_generate_U.icp in INTERCONNECT. Then drag and drop the file lumfoundry_CML.cml. Close the project file (temporary step).
- Open step1_generate_U.icp again and run the script step1_generate_U.py

In this step the hopping Hamiltonian is calculated for 6 time steps and used to generate a photonic circuit for each time step labelled 'U_0.1.icp' to 'U_5.icp'.

### Step 2: Perform simulations in qINTERCONNECT

- Open U_0.1.icp and run step_2_1.py
- Open step2_2_U_Uinv_0.1.icp and run step_2_2.py
- Open step2_3_U_Uinv_QFT_0.1.icp and run step_2_3.py

In the first step, the circuits generated in the previous step are simulated

The script step_2_1.py loads each of the circuits generated in the previous step and performs multiple simulations for each time step. Firstly, an input state of \( \vert{1,1,1,0}\rangle\), consisting of one indistinguishable photon in each of channels 1, 2, and 3 and no photon in channel 4, is injected to the circuit and a measurement of the first channel is performed. The probabilities of success for these measurements are plotted for each photon number outcome as red dots in the figure below. Additionally, the simulations are repeated with distinguishable photons, plotted as blue dots.

The results are similar to those in Figure 3 of [2]. Next, additional simulations are performed with the input state of indistinguishable photons \( \vert{1,1,1,0}\rangle\) used above except this time all four output channels are measured. The probabilities of measuring each possible outcome are shown below.

In the second step, the simulation file contains a circuit for the unitary matrix above at *t = 0.1 s* followed by a circuit for the inverse of the unitary matrix. Under ideal conditions, these circuits combined represent the identity matrix, as shown below and in agreement with Figure 3 of [2].

In the third step, the simulation is similar to the one above except that the circuit following the unitary matrix now represents the inverse of the unitary matrix followed by a three-channel quantum Fourier transform.

### Step 3: Parameter sweep

- Open step3.icp and run step_3.py

The script step_3.py generates a new set of circuits with variations in the directional coupler gaps. Instead of an ideal value of 200 nm, gap values are taken from a normal distribution of gap values centered on 200 nm with a standard deviation of 10 nm. The simulation from step 2 where all channels are output measured is then repeated with the new circuits and compared to the ideal case, show below.

## Important model settings

Description of important objects and settings used in this model

The simulation parameters can be changed by modifying the parameter dictionary for each simulation, for instance:

### Common settings

**interconnect_file **: the name of the INTERCONNECT file "U_0.1.icp" containing the unitary matrix circuit.

**circuit, nchannel, mode per channel **: every circuit simulated in qINTERCONNECT must be contained in a compound element, so we provide the name of this compound element (circuit) as “Unitary Matrix”. The unitary matrix circuit contains four channels, so nchannel is set to 4 and mode_per_channel is set to 1.

**Input_ports, output_ports **: the names of the input and output ports as they are named in the INTERCONNECT circuit.

**input_dims **: the maximum number of photons to be simulated in each mode, which in this case is 3.

**frequency_points**: an optional parameter used for multi-frequency simulations.

**ket_in **: the input state is specified in the Fock basis under dictionary “ket_in”, where the key is a tuple corresponding to the number of photons in each channel, and the value is a list corresponding to the probability amplitude of the state in polar coordinates, ordered by the magnitude and then the angle. Here the input in the signal channel corresponds to the state \( \vert{1,1,1,0}\rangle \).

Distinguishable photons can be simulated by using the frequency degree of freedom. For instance, setting the number of frequency_points to 3, the state \( \vert{1,1,1,0}\rangle \) consisting of distinguishable photons can be represented as \( \vert{1,0,0,0}\rangle_1 \vert{0,1,0,0}\rangle_2 \vert{0,0,1,0}\rangle_3\) which corresponds to a ket_in value of (1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0): [1, 0].

**measurement_mask **: refers to what we are measuring at the output of each mode. A value of 100 means we do not directly measure the mode, but we want to know its state at the end of the simulation, whereas a value of 200 means the information in the mode is traced over, so that it doesn’t appear in the final density matrix. Any other value means we only keep the outcomes where such a number of photons is detected.

**ket_out **: the expected output state for the fidelity to be calculated against.

## Updating the model with your parameters

Instructions for updating the model based on your device parameters

### Changing circuit parameters

The directional coupler CML element can be replaced with a custom directional coupler by modifying the contents of the custom_builder.py script file. Furthermore, the file can be modified to support custom beamsplitter configurations, provided that they can be described by

$$U_{BS} = ie^{i\theta}e^{i\Phi_{BS}}\begin{pmatrix}

\text{sin}(\theta)e^{i\phi} & \text{cos}(\theta) \\

\text{cos}(\theta)e^{i\phi} & -\text{sin}(\theta) \\

\end{pmatrix} $$

up to a global phase. If using a different directional coupler or beamsplitter configuration, the global phase of the beamsplitter will be affected, and the value of this updated phase must be changed in the custom_builder.py script.

## Taking the model further

Information and tips for users that want to further customize the model

The UnitaryCircuit method supports the creation of arbitrary unitary matrices, and can be used to simulate other Hamiltonians such as the long-range Hamiltonian in Figure 3b of [2].

## Additional resources

Additional documentation, examples and training material

### Related publications

- W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer and I. A. Walmsley, "Optimal design for universal multiport interferometers", Optica Vol.
**3**, No. 12, p. 1460-1465, 2016. - F. H. B. Somhorst
*et al*, "Quantum simulation of thermodynamics in an integrated quantum photonic processor", Nat Commun, Vol. 14, No. 3895, 2023.