Simulating this circuit using PennyLane is easy; we can simply read off the gates from left to right, and convert it into a QNode. """ import numpy as np # set the random seed np.random.seed(42) # import PennyLane import pennylane as qml ###################################################################### # We must define the unitary matrix we would like to embed in the circuit. # We will use SciPy to generate a Haar-random unitary: from scipy.stats import unitary_group # define the linear interferometer U = unitary_group.rvs(4) print(U) ###################################################################### # We can now use this to construct the circuit, choosing a compatible # device. For the simulation, we can use the Strawberry Fields # Gaussian backend. This backend is perfectly suited for simulation of GBS, # as the initial states are Gaussian, and all gates transform Gaussian states to other # Gaussian states. n_wires = 4 cutoff = 10 dev = qml.device("strawberryfields.gaussian", wires=n_wires, cutoff_dim=cutoff) @qml.qnode(dev) def gbs_circuit(): # prepare the input squeezed states for i in range(n_wires): qml.Squeezing(1.0, 0.0, wires=i) # linear interferometer qml.InterferometerUnitary(U, wires=range(n_wires)) return qml.probs(wires=range(n_wires)) ###################################################################### # A couple of things to note in this particular example: # # 1. To prepare the input single mode squeezed vacuum state :math:`\ket{re^{i\phi}}`, # where :math:`r = 1` and :math:`\phi=0`, we # apply a squeezing gate (:class:`~pennylane.Squeezing`) to each of the wires (initially in # the vacuum state). # # 2. Next we apply the linear interferometer to all four wires using # :class:`~pennylane.Interferometer` and the unitary matrix ``U``. This operator # decomposes the unitary matrix representing the linear interferometer into single-mode # rotation gates (:class:`~pennylane.PhaseShift`) and two-mode beamsplitters # (:class:`~pennylane.Beamsplitter`). After applying the interferometer, we will denote the # output state by :math:`\ket{\psi'}`. # # 3. GBS takes place physically in an infinite-dimensional Hilbert space, # which is not practical for simulation. We need to set an upper limit on the maximum # number of photons we can detect. This is the # ``cutoff`` value we defined above; we will only be considering detection events # containing 0 to 9 photons per mode. # # We can now execute the QNode, and extract the resulting probability distribution: probs = gbs_circuit().reshape([cutoff] * n_wires) print(probs.shape) ###################################################################### # For example, element ``[1,2,0,1]`` represents the probability of # detecting 1 photon on wire # ``0`` and wire ``3``, and 2 photons at wire ``1``, i.e., the value # # .. math:: \text{prob}(1,2,0,1) = \left|\braketD{1,2,0,1}{\psi'}\right|^2. # # Let's extract and view the probabilities of measuring various Fock states. # Fock states to measure at output measure_states = [(0, 0, 0, 0), (1, 1, 0, 0), (0, 1, 0, 1), (1, 1, 1, 1), (2, 0, 0, 0)] # extract the probabilities of calculating several # different Fock states at the output, and print them out for i in measure_states: print(f"|{''.join(str(j) for j in i)}>: {probs[i]}") ###################################################################### # The GBS Distribution # -------------------- # # Hamilton et al. [#hamilton2017]_ showed that the probability of # measuring a final state containing only 0 or 1 photons per mode is given by # # .. math:: # # \left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 = # \frac{\left|\text{Haf}[(U(\bigoplus_i\tanh(r_i))U^T)]_{st}\right|^2}{\prod_{i=1}^N \cosh(r_i)} # # i.e., the sampled single-photon probability distribution is proportional to the **hafnian** of a # submatrix of :math:`U(\bigoplus_i\tanh(r_i))U^T`. # # .. note:: # # The hafnian of a matrix is defined by # # .. math:: \text{Haf}(A) = \sum_{\sigma \in \text{PMP}_{2N}}\prod_{i=1}^N A_{\sigma(2i-1)\sigma(2i)}, # # where :math:`\text{PMP}_{2N}` is the set of all perfect matching permutations of :math:`2N` elements. In graph theory, the # hafnian calculates the number of perfect `matchings #