Repository: sequence-toolbox/SeQUeNCe (language: python)

Issue:
Memory simulation update
Adding all of the memory simulation work (absorptive memories and Fock state representation) to the main development branch


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index effbb1af..c9d07142 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -14,6 +14,11 @@ jobs:
       matrix:
         python-version: ['3.8', '3.9', '3.10']
         mpi: [mpich, openmpi]
+        exclude:
+          - python-version: '3.8'
+            mpi: mpich
+          - python-version: '3.9'
+            mpi: mpich
 
     steps:
     - uses: actions/checkout@v2
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b6bfe332..efe1641f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -145,3 +145,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Removed
 - residual `json5` references in examples
+
+## [0.6.0]
+### Added
+- Fock state density matrix encoding
+- Quantum manager to use fock state density matrix
+  - Add truncation attribute to quantum managers; Hilbert space dimension is truncation + 1
+  - Supports error channels based on kraus operators and direct operator action on state
+- Hardware models to use fock state density matrix
+  - Absorptive quantum memory built on AFC
+  - Interference and direct detection for state tomography
+
+### Changed
+- Moved around and updated some error checking
+- Temporarily removed mpich testing for python 3.8 and 3.9 (broken)
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 659726be..2fcbe3b2 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -22,7 +22,7 @@
 author = 'Xiaoliang Wu, Joaquin Chung, Alexander Kolar, Eugene Wang, Tian Zhong, Rajkumar Kettimuthu, Martin Suchara'
 
 # The full version, including alpha/beta/rc tags
-release = '0.5.4'
+release = '0.6.0'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/example/absorptive_experiment.py b/example/absorptive_experiment.py
index 53381184..613b4b6a 100644
--- a/example/absorptive_experiment.py
+++ b/example/absorptive_experiment.py
@@ -1,67 +1,163 @@
-from typing import List, Dict, Any
+"""The main script for simulating experiment of entanglement generation between two remote AFC absorptive quantum memories.
+
+There are 4 nodes involved: 2 memory nodes, 1 entangling node (for BSM) and 1 measurement node (for measurement of retrieved photonic state).
+Each memory node is connected with both entangling node and measurement node, but there is no direct connection between memory nodes.
+
+Each memory node contains an AFC memory instance and an SPDC source instance.
+The entangling node contains a QSDetectorFockInterference instance (BSM device with a beamsplitter and two photon detectors behind).
+The measurement node contians a QSDetectorFockDirect instance and a QSDetectorFockInterference instance, for measurement of 
+    diagonal and off-diagonal elements of the effective 4-d density matrix, respectively.
+"""
+
+from typing import List, Callable, TYPE_CHECKING
 from pathlib import Path
+from copy import copy
+
+if TYPE_CHECKING:
+    from sequence.components.photon import Photon
 
 from json import dump
 import numpy as np
 
-from sequence.components.bsm import make_bsm
+from sequence.kernel.event import Event
+from sequence.kernel.process import Process
+from sequence.kernel.timeline import Timeline
+from sequence.kernel.quantum_manager import FOCK_DENSITY_MATRIX_FORMALISM
 from sequence.components.detector import QSDetectorFockDirect, QSDetectorFockInterference
 from sequence.components.light_source import SPDCSource
 from sequence.components.memory import AbsorptiveMemory
+from sequence.components.optical_channel import QuantumChannel
 from sequence.components.photon import Photon
-from sequence.kernel.event import Event
-from sequence.kernel.process import Process
-from sequence.kernel.timeline import Timeline
 from sequence.topology.node import Node
-from sequence.topology.topology import Topology
 from sequence.protocol import Protocol
-from sequence.utils.encoding import absorptive
+from sequence.kernel.quantum_utils import *  # only for manual calculation and should not be used in simulation
+
+
+# define simulation constants
+
+# quantum manager
+TRUNCATION = 2  # truncation of Fock space (=dimension-1)
+
+# photon sources
+TELECOM_WAVELENGTH = 1436  # telecom band wavelength of SPDC source idler photon
+WAVELENGTH = 606  # wavelength of AFC memory resonant absorption, of SPDC source signal photon
+SPDC_FREQUENCY = 80e6  # frequency of both SPDC sources' photon creation (same as memory frequency and detector count rate)
+MEAN_PHOTON_NUM1 = 0.1  # mean photon number of SPDC source on node 1
+MEAN_PHOTON_NUM2 = 0.1  # mean photon number of SPDC source on node 2
+
+# detectors
+BSM_DET1_EFFICIENCY = 0.6  # efficiency of detector 1 of BSM
+BSM_DET2_EFFICIENCY = 0.6  # efficiency of detector 2 of BSM
+BSM_DET1_DARK = 150  # Dark count rate (Hz)
+BSM_DET2_DARK = 150
+MEAS_DET1_EFFICIENCY = 0.6  # efficiency of detector 1 of DM measurement
+MEAS_DET2_EFFICIENCY = 0.6  # efficiency of detector 2 of DM measurement
+MEAS_DET1_DARK = 150
+MEAS_DET2_DARK = 150
+
+# fibers
+DIST_ANL_ERC = 20  # distance between ANL and ERC, in km
+DIST_HC_ERC = 20  # distance between HC and ERC, in km
+ATTENUATION = 0.2  # attenuation rate of optical fibre (in dB/km)
+DELAY_CLASSICAL = 5e-3  # delay for classical communication between BSM node and memory nodes (in s)
+
+# memories
+MODE_NUM = 100  # number of temporal modes of AFC memory (same for both memories)
+MEMO_FREQUENCY1 = SPDC_FREQUENCY  # frequency of memory 1
+MEMO_FREQUENCY2 = SPDC_FREQUENCY  # frequency of memory 2
+ABS_EFFICIENCY1 = 0.35  # absorption efficiency of AFC memory 1
+ABS_EFFICIENCY2 = 0.35  # absorption efficiency of AFC memory 2
+PREPARE_TIME1 = 0  # time required for AFC structure preparation of memory 1
+PREPARE_TIME2 = 0  # time required for AFC structure preparation of memory 2
+COHERENCE_TIME1 = -1  # spin coherence time for AFC memory 1 spinwave storage, -1 means infinite time
+COHERENCE_TIME2 = -1  # spin coherence time for AFC memory 2 spinwave storage, -1 means infinite time
+AFC_LIFETIME1 = -1  # AFC structure lifetime of memory 1, -1 means infinite time
+AFC_LIFETIME2 = -1  # AFC structure lifetime of memory 2, -1 means infinite time
+DECAY_RATE1 = 4.3e-8  # retrieval efficiency decay rate for memory 1
+DECAY_RATE2 = 4.3e-8  # retrieval efficiency decay rate for memory 2
+
+# experiment settings
+time = int(1e12)
+calculate_fidelity_direct = True
+calculate_rate_direct = True
+num_direct_trials = 200
+num_bs_trials_per_phase = 50
+phase_settings = np.linspace(0, 2*np.pi, num=20, endpoint=False)
 
 
-# define constants
-FREQUENCY = 80e6
-ABS_EFFICIENCY = 1.0
-MODE_NUM = 10000
-COHERENCE_TIME = 1
-TELECOM_WAVELENGTH = 1000
-WAVELENGTH = 500
-OVERLAP_ERR = 0
+# function to generate standard pure Bell state for fidelity calculation
+def build_bell_state(truncation, sign, phase=0, formalism="dm"):
+    """Generate standard Bell state which is heralded in ideal BSM.
 
-# experiment settings
-num_direct_trials = 10
-phase_settings = np.linspace(0, 2*np.pi, num=10, endpoint=False)
+    For comparison with results from imperfect parameter choices.
+    """
+
+    basis0 = np.zeros(truncation+1)
+    basis0[0] = 1
+    basis1 = np.zeros(truncation+1)
+    basis1[1] = 1
+    basis10 = np.kron(basis1, basis0)
+    basis01 = np.kron(basis0, basis1)
+    
+    if sign == "plus":
+        ket = (basis10 + np.exp(1j*phase)*basis01)/np.sqrt(2)
+    elif sign == "minus":
+        ket = (basis10 - np.exp(1j*phase)*basis01)/np.sqrt(2)
+    else:
+        raise ValueError("Invalid Bell state sign type " + sign)
+
+    dm = np.outer(ket, ket.conj())
+
+    if formalism == "dm":
+        return dm
+    elif formalism == "ket":
+        return ket
+    else:
+        raise ValueError("Invalid quantum state formalism " + formalism)
 
 
-# for absorptive quantum memory
-def efficiency(_: int) -> float:
-    return 1.0
+# retrieval efficiency as function of storage time for absorptive quantum memory, using exponential decay model
+def efficiency1(t: int) -> float:
+    return np.exp(-t*DECAY_RATE1)
+
+
+def efficiency2(t: int) -> float:
+    return np.exp(-t*DECAY_RATE2)
+
+
+def add_channel(node1: Node, node2: Node, timeline: Timeline, **kwargs):
+    name = "_".join(["qc", node1.name, node2.name])
+    qc = QuantumChannel(name, timeline, **kwargs)
+    qc.set_ends(node1, node2.name)
+    return qc
 
 
 # protocol to control photon emission on end node
 class EmitProtocol(Protocol):
-    def __init__(self, own: "EndNode", name: str, other_node: str, num_qubits: int, source_name: str, memory_name: str):
+    def __init__(self, own: "EndNode", name: str, other_node: str, photon_pair_num: int,
+                 source_name: str, memory_name: str):
         """Constructor for Emission protocol.
 
         Args:
             own (EndNode): node on which the protocol is located.
             name (str): name of the protocol instance.
-            other_node (str): name of the other node entangling photons
-            num_qubits (int): number of qubits to send in one execution.
-            delay_time (int): time to wait before re-starting execution.
+            other_node (str): name of the other node to generate entanglement with
+            photon_pair_num (int): number of output photon pulses to send in one execution.
             source_name (str): name of the light source on the node.
             memory_name (str): name of the memory on the node.
         """
 
         super().__init__(own, name)
         self.other_node = other_node
-        self.num_qubits = num_qubits
+        self.num_output = photon_pair_num
         self.source_name = source_name
         self.memory_name = memory_name
 
     def start(self):
-        self.own.components[self.memory_name]._prepare_AFC()
+        if not self.own.components[self.memory_name].is_prepared:
+            self.own.components[self.memory_name]._prepare_AFC()
         
-        states = [None] * self.num_qubits  # TODO: rewrite spdc class?
+        states = [None] * self.num_output  # for Fock encoding only list length matters and list elements do not matter
         source = self.own.components[self.source_name]
         source.emit(states)
 
@@ -70,26 +166,28 @@ def received_message(self, src: str, msg):
 
 
 class EndNode(Node):
-    """Node for each end of the network.
+    """Node for each end of the network (the memory node).
 
     This node stores an SPDC photon source and a quantum memory.
+    The properties of attached devices are made customizable for each individual node.
     """
 
-    def __init__(self, name: str, timeline: "Timeline", other_node: str, bsm_node: str, measure_node: str):
+    def __init__(self, name: str, timeline: "Timeline", other_node: str, bsm_node: str, measure_node: str,
+                 mean_photon_num: float, spdc_frequency: float, memo_frequency: float, abs_effi: float,
+                 afc_efficiency: Callable, mode_number: int):
         super().__init__(name, timeline)
 
         self.bsm_name = bsm_node
         self.meas_name = measure_node
 
         # hardware setup
-        spdc_name = name + ".spdc_source"
-        memo_name = name + ".memory"
-        spdc = SPDCSource(spdc_name, timeline, wavelengths=[TELECOM_WAVELENGTH, WAVELENGTH],
-                          frequency=FREQUENCY, encoding_type=absorptive, mean_photon_num=0.2)
-        memory = AbsorptiveMemory(memo_name, timeline, fidelity=0.85, frequency=FREQUENCY,
-                                  absorption_efficiency=ABS_EFFICIENCY, mode_number=MODE_NUM,
-                                  coherence_time=COHERENCE_TIME, wavelength=WAVELENGTH, overlap_error=OVERLAP_ERR,
-                                  efficiency=efficiency, prepare_time=0, destination=measure_node)
+        self.spdc_name = name + ".spdc_source"
+        self.memo_name = name + ".memory"
+        spdc = SPDCSource(self.spdc_name, timeline, wavelengths=[TELECOM_WAVELENGTH, WAVELENGTH],
+                          frequency=spdc_frequency, mean_photon_num=mean_photon_num)
+        memory = AbsorptiveMemory(self.memo_name, timeline, frequency=memo_frequency,
+                                  absorption_efficiency=abs_effi, afc_efficiency=afc_efficiency,
+                                  mode_number=mode_number, wavelength=WAVELENGTH, destination=measure_node)
         self.add_component(spdc)
         self.add_component(memory)
         spdc.add_receiver(self)
@@ -97,7 +195,7 @@ def __init__(self, name: str, timeline: "Timeline", other_node: str, bsm_node: s
         memory.add_receiver(self)
 
         # protocols
-        self.emit_protocol = EmitProtocol(self, name + ".emit_protocol", other_node, MODE_NUM, spdc_name, memo_name)
+        self.emit_protocol = EmitProtocol(self, name + ".emit_protocol", other_node, mode_number, self.spdc_name, self.memo_name)
 
     def get(self, photon: "Photon", **kwargs):
         dst = kwargs.get("dst")
@@ -110,44 +208,54 @@ def get(self, photon: "Photon", **kwargs):
 
 
 class EntangleNode(Node):
-    def __init__(self, name: str, timeline: "Timeline"):
+    def __init__(self, name: str, timeline: "Timeline", src_list: List[str]):
         super().__init__(name, timeline)
 
         # hardware setup
-        bsm_name = name + ".bsm"
-        bsm = make_bsm(bsm_name, timeline, "absorptive")
+        self.bsm_name = name + ".bsm"
+        # assume no relative phase between two input optical paths
+        bsm = QSDetectorFockInterference(self.bsm_name, timeline, src_list)
         self.add_component(bsm)
         bsm.attach(self)
-        self.set_first_component(bsm_name)
+        self.set_first_component(self.bsm_name)
+        self.resolution = max([d.time_resolution for d in bsm.detectors])
 
-        self.resolution = bsm.resolution
-        self.bsm_times = [[], []]
+        # detector parameter setup
+        bsm.set_detector(0, efficiency=BSM_DET1_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=BSM_DET1_DARK)
+        bsm.set_detector(1, efficiency=BSM_DET2_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=BSM_DET1_DARK)
 
-    def bsm_update(self, bsm, info: Dict[str, Any]):
-        self.bsm_times[info['res']].append(info['time'])
+    def receive_qubit(self, src: str, qubit) -> None:
+        self.components[self.first_component_name].get(qubit, src=src)
 
-    def get_valid_bins(self, start_time: int, num_bins: int, frequency: float):
-        """Computes time bins containing a BSM measurement.
+    def get_detector_entries(self, detector_name: str, start_time: int, num_bins: int, frequency: float):
+        """Returns detection events for density matrix measurement. Used to determine BSM result.
 
         Args:
+            detector_name (str): name of detector to get measurements from.
             start_time (int): simulation start time of when photons received.
             num_bins (int): number of arrival bins
             frequency (float): frequency of photon arrival (in Hz).
 
         Returns:
-            List[int]: list of length num_bins containing 0/1.
-                List element is 0 for an unsuccessful measurement and 1 for successful.
+            List[int]: list of length (num_bins) with result for each time bin.
         """
 
-        return_bins = [0] * num_bins
+        trigger_times = self.components[detector_name].get_photon_times()
+        return_res = [0] * num_bins
+
+        for time in trigger_times[0]:
+            closest_bin = int(round((time - start_time) * frequency * 1e-12))
+            expected_time = (float(closest_bin) * 1e12 / frequency) + start_time
+            if abs(expected_time - time) < self.resolution and 0 <= closest_bin < num_bins:
+                return_res[closest_bin] += 1
 
-        for time in self.bsm_times[0] + self.bsm_times[1]:
+        for time in trigger_times[1]:
             closest_bin = int(round((time - start_time) * frequency * 1e-12))
             expected_time = (float(closest_bin) * 1e12 / frequency) + start_time
             if abs(expected_time - time) < self.resolution and 0 <= closest_bin < num_bins:
-                return_bins[closest_bin] = not return_bins[closest_bin]
+                return_res[closest_bin] += 2
 
-        return return_bins
+        return return_res
 
 
 class MeasureNode(Node):
@@ -158,16 +266,23 @@ def __init__(self, name: str, timeline: "Timeline", other_nodes: List[str]):
         direct_detector = QSDetectorFockDirect(self.direct_detector_name, timeline, other_nodes)
         self.add_component(direct_detector)
         direct_detector.attach(self)
-        self.set_first_component(self.direct_detector_name)
 
         self.bs_detector_name = name + ".bs"
-        bs_detector = QSDetectorFockInterference(self.bs_detector_name, tl)
+        bs_detector = QSDetectorFockInterference(self.bs_detector_name, timeline, other_nodes)
         self.add_component(bs_detector)
         bs_detector.add_receiver(self)
 
+        self.set_first_component(self.direct_detector_name)
+
         # time resolution of SPDs
         self.resolution = max([d.time_resolution for d in direct_detector.detectors + bs_detector.detectors])
 
+        # detector parameter setup
+        direct_detector.set_detector(0, efficiency=MEAS_DET1_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=MEAS_DET1_DARK)
+        direct_detector.set_detector(1, efficiency=MEAS_DET2_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=MEAS_DET2_DARK)
+        bs_detector.set_detector(0, efficiency=MEAS_DET1_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=MEAS_DET1_DARK)
+        bs_detector.set_detector(1, efficiency=MEAS_DET2_EFFICIENCY, count_rate=SPDC_FREQUENCY, dark_count=MEAS_DET2_DARK)
+
     def receive_qubit(self, src: str, qubit) -> None:
         self.components[self.first_component_name].get(qubit, src=src)
 
@@ -175,7 +290,7 @@ def set_phase(self, phase: float):
         self.components[self.bs_detector_name].set_phase(phase)
 
     def get_detector_entries(self, detector_name: str, start_time: int, num_bins: int, frequency: float):
-        """Computes distribution of detection events for density matrix.
+        """Returns detection events for density matrix measurement.
 
         Args:
             detector_name (str): name of detector to get measurements from.
@@ -184,7 +299,7 @@ def get_detector_entries(self, detector_name: str, start_time: int, num_bins: in
             frequency (float): frequency of photon arrival (in Hz).
 
         Returns:
-            List[int]: list of length (duration * 1e-12 * frequency) with result for each time bin.
+            List[int]: list of length (num_bins) with result for each time bin.
         """
 
         trigger_times = self.components[detector_name].get_photon_times()
@@ -206,72 +321,321 @@ def get_detector_entries(self, detector_name: str, start_time: int, num_bins: in
 
 
 if __name__ == "__main__":
-    # open file to store experiment results
-    Path("results").mkdir(parents=True, exist_ok=True)
-    filename = "results/absorptive.json"
-    fh = open(filename, 'w')
 
-    """Run Simulation"""
+    """Setup Simulation"""
+
+    """ For generating effective fidelity 2-d plot
+    dim = 100 
+    
+    mean_num1_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+    mean_num2_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+
+    fidelity_mat = np.zeros((dim,dim))
+    idx1 = -1
+    idx2 = -1
+    for mean_num1 in mean_num1_list:
+        idx1 += 1
+
+        for mean_num2 in mean_num2_list:
+            idx2 += 1
+            
+            tl = Timeline(time, formalism=FOCK_DENSITY_MATRIX_FORMALISM, truncation=TRUNCATION)
+
+            anl_name = "Argonne"
+            hc_name = "Harper Court"
+            erc_name = "Eckhardt Research Center BSM"
+            erc_2_name = "Eckhardt Research Center Measurement"
+            seeds = [1, 2, 3, 4]
+            src_list = [anl_name, hc_name]  # the list of sources, note the order
+        
+            anl = EndNode(anl_name, tl, hc_name, erc_name, erc_2_name, mean_photon_num=mean_num1,
+                          spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY1, abs_effi=ABS_EFFICIENCY1,
+                          afc_efficiency=efficiency1, mode_number=MODE_NUM)
+            hc = EndNode(hc_name, tl, anl_name, erc_name, erc_2_name, mean_photon_num=mean_num2,
+                         spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY2, abs_effi=ABS_EFFICIENCY2,
+                         afc_efficiency=efficiency2, mode_number=MODE_NUM)
+            erc = EntangleNode(erc_name, tl, src_list)
+            erc_2 = MeasureNode(erc_2_name, tl, src_list)
+        
+            for seed, node in zip(seeds, [anl, hc, erc, erc_2]):
+                node.set_seed(seed)
+        
+            # extend fiber lengths to be equivalent
+            fiber_length = max(DIST_ANL_ERC, DIST_HC_ERC)
+        
+            qc1 = add_channel(anl, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+            qc2 = add_channel(hc, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+            qc3 = add_channel(anl, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+            qc4 = add_channel(hc, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        
+            tl.init()
+        
+            # Pre-simulation explicit calculation of entanglement fidelity upon successful BSM
+        
+            # use non-transmitted Photon as interface with existing methods in SeQUeNCe
+            spdc_anl = anl.components[anl.spdc_name]
+            spdc_hc = hc.components[hc.spdc_name]
+            memo_anl = anl.components[anl.memo_name]
+            memo_hc = hc.components[hc.memo_name]
+            channel_anl = anl.qchannels[erc_name]
+            channel_hc = hc.qchannels[erc_name]
+            bsm = erc.components[erc.bsm_name]
+        
+            # photon0: idler, photon1: signal
+            photon0_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[0], location=spdc_anl, encoding_type=spdc_anl.encoding_type, use_qm=True)
+            photon1_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[1], location=spdc_anl, encoding_type=spdc_anl.encoding_type, use_qm=True)
+            # set shared state to squeezed state
+            state_spdc_anl = spdc_anl._generate_tmsv_state()
+            keys = [photon0_anl.quantum_state, photon1_anl.quantum_state]
+            tl.quantum_manager.set(keys, state_spdc_anl)
+        
+            photon0_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[0], location=spdc_hc, encoding_type=spdc_hc.encoding_type, use_qm=True)
+            photon1_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[1], location=spdc_hc, encoding_type=spdc_hc.encoding_type, use_qm=True)
+            # set shared state to squeezed state
+            state_spdc_hc = spdc_hc._generate_tmsv_state()
+            keys = [photon0_hc.quantum_state, photon1_hc.quantum_state]
+            tl.quantum_manager.set(keys, state_spdc_hc)
+        
+            # photon loss upon absorption by memories
+            key_anl_memo = photon1_anl.quantum_state
+            loss_anl_memo = 1 - memo_anl.absorption_efficiency
+            tl.quantum_manager.add_loss(key_anl_memo, loss_anl_memo)
+            key_hc_memo = photon1_hc.quantum_state
+            loss_hc_memo = 1 - memo_hc.absorption_efficiency
+            tl.quantum_manager.add_loss(key_hc_memo, loss_hc_memo)
+    
+            # transmission loss through optical fibres (overwrites previous variables)
+            key_anl_pho = photon0_anl.quantum_state
+            loss_anl_pho = channel_anl.loss
+            tl.quantum_manager.add_loss(key_anl_pho, loss_anl_pho)
+            key_hc_pho = photon0_hc.quantum_state
+            loss_hc_pho = channel_anl.loss
+            tl.quantum_manager.add_loss(key_hc_pho, loss_hc_pho)
+        
+            # QSDetector measurement and remaining state after partial trace
+            povms = bsm.povms
+            povm_tuple = tuple([tuple(map(tuple, povm)) for povm in povms])
+            keys = [photon0_anl.quantum_state, photon0_hc.quantum_state]
+            new_state, all_keys = tl.quantum_manager._prepare_state(keys)
+            indices = tuple([all_keys.index(key) for key in keys])
+            state_tuple = tuple(map(tuple, new_state))
+            states, probs = measure_multiple_with_cache_fock_density(state_tuple, indices, len(all_keys), povm_tuple, tl.quantum_manager.truncation)
+            state_plus, state_minus = states[1], states[2]
+                    
+            # assign remaining state (minus as example)
+            if len(keys) < len(all_keys):
+                indices = tuple([all_keys.index(key) for key in keys])
+                new_state_tuple = tuple(map(tuple, state_plus))
+                remaining_state = density_partial_trace(new_state_tuple, indices, len(all_keys), tl.quantum_manager.truncation)
+                
+            # effective Bell state generated 
+            def effective_state(state):
+                state[0][0] = 0
+                state = state/np.trace(state)
+                
+                return state
+            
+            remaining_state_eff = effective_state(remaining_state)
+            
+            # calculate the fidelity with reference Bell state
+            bell_plus = build_bell_state(tl.quantum_manager.truncation, "plus")
+            bell_minus = build_bell_state(tl.quantum_manager.truncation, "minus")
+            fidelity = np.trace(remaining_state_eff.dot(bell_minus)).real
+            
+            fidelity_mat[idx1][idx2] = fidelity
+                
+        idx2 = -1
+            
+    # plot the 2-d matrix
+    
+    num1, num2 = np.meshgrid(mean_num1_list, mean_num2_list)
+    
+    fig = plt.figure(figsize=(8,6))
+    im = plt.imshow(fidelity_mat, cmap=plt.cm.RdBu, extent=(0.02, 0.1, 0.1, 0.02))
+    plt.colorbar(im)  
+    plt.title("Effective Fidelity with Mean Photon Number")
+    plt.xlabel("mean photon number of source 1 $\mu_1$")
+    plt.ylabel("mean photon number of source 2 $\mu_2$")
+    plt.show()
+    """
 
-    tl = Timeline(1e12, 'density_matrix')
-    tl.seed(0)
+    tl = Timeline(time, formalism=FOCK_DENSITY_MATRIX_FORMALISM, truncation=TRUNCATION)
 
     anl_name = "Argonne"
     hc_name = "Harper Court"
-    erc_name = "Eckhardt Research Center"
-    erc_2_name = "Eckhardt Research Center 2"
-
-    anl = EndNode(anl_name, tl, hc_name, erc_name, erc_2_name)
-    hc = EndNode(hc_name, tl, anl_name, erc_name, erc_2_name)
-    erc = EntangleNode(erc_name, tl)
-    erc_2 = MeasureNode(erc_2_name, tl, [anl_name, hc_name])
-
-    topo = Topology("Experiment Topo", tl)
-    for node in [anl, hc, erc, erc_2]:
-        topo.add_node(node)
-    topo.add_quantum_channel(anl_name, erc_name, distance=20, attenuation=0)
-    topo.add_quantum_channel(hc_name, erc_name, distance=20, attenuation=0)
-    topo.add_quantum_channel(anl_name, erc_2_name, distance=20, attenuation=0)
-    topo.add_quantum_channel(hc_name, erc_2_name, distance=20, attenuation=0)
+    erc_name = "Eckhardt Research Center BSM"
+    erc_2_name = "Eckhardt Research Center Measurement"
+    seeds = [1, 2, 3, 4]
+    src_list = [anl_name, hc_name]  # the list of sources, note the order
+
+    anl = EndNode(anl_name, tl, hc_name, erc_name, erc_2_name, mean_photon_num=MEAN_PHOTON_NUM1,
+                  spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY1, abs_effi=ABS_EFFICIENCY1,
+                  afc_efficiency=efficiency1, mode_number=MODE_NUM)
+    hc = EndNode(hc_name, tl, anl_name, erc_name, erc_2_name, mean_photon_num=MEAN_PHOTON_NUM2,
+                 spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY2, abs_effi=ABS_EFFICIENCY2,
+                 afc_efficiency=efficiency2, mode_number=MODE_NUM)
+    erc = EntangleNode(erc_name, tl, src_list)
+    erc_2 = MeasureNode(erc_2_name, tl, src_list)
+
+    for seed, node in zip(seeds, [anl, hc, erc, erc_2]):
+        node.set_seed(seed)
+
+    # extend fiber lengths to be equivalent
+    fiber_length = max(DIST_ANL_ERC, DIST_HC_ERC)
+
+    qc1 = add_channel(anl, erc, tl, distance=fiber_length, attenuation=ATTENUATION)
+    qc2 = add_channel(hc, erc, tl, distance=fiber_length, attenuation=ATTENUATION)
+    qc3 = add_channel(anl, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION)
+    qc4 = add_channel(hc, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION)
 
     tl.init()
 
-    # calculations for when to start protocol
-    delay_anl = topo.nodes[anl_name].qchannels[erc_name].delay
-    delay_hc = topo.nodes[hc_name].qchannels[erc_name].delay
-    time_anl = max(delay_anl, delay_hc) - delay_anl
-    time_hc = max(delay_anl, delay_hc) - delay_hc
+    # calculate start time for protocol
+    # since fiber lengths equal, both start at 0
+    start_time_anl = start_time_hc = 0
 
     # calculations for when to start recording measurements
-    start_time_bsm = time_anl + delay_anl
-    mem = anl.get_components_by_type("AbsorptiveMemory")[0]
+    delay_anl = anl.qchannels[erc_2_name].delay
+    delay_hc = hc.qchannels[erc_2_name].delay
+    assert delay_anl == delay_hc
+    start_time_bsm = start_time_anl + delay_anl
+    mem = anl.components[anl.memo_name]
     total_time = mem.total_time
-    start_time_meas = time_anl + total_time + delay_anl
+    start_time_meas = start_time_anl + total_time + delay_anl
 
     results_direct_measurement = []
-    results_bs_measurement = []
+    results_bs_measurement = [[] for _ in phase_settings]
+
+    """Pre-simulation explicit calculation of effective entanglement fidelity upon successful BSM"""
+
+    if calculate_fidelity_direct:
+        # use non-transmitted Photon as interface with existing methods in SeQUeNCe
+        spdc_anl = anl.components[anl.spdc_name]
+        spdc_hc = hc.components[hc.spdc_name]
+        memo_anl = anl.components[anl.memo_name]
+        memo_hc = hc.components[hc.memo_name]
+        channel_anl = anl.qchannels[erc_name]
+        channel_hc = hc.qchannels[erc_name]
+        bsm = erc.components[erc.bsm_name]
+
+        # photon0: idler, photon1: signal
+        photon0_anl = Photon("", tl, wavelength=spdc_anl.wavelengths[0], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        photon1_anl = Photon("", tl, wavelength=spdc_anl.wavelengths[1], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_anl = spdc_anl._generate_tmsv_state()
+        keys = [photon0_anl.quantum_state, photon1_anl.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_anl)
+
+        photon0_hc = Photon("", tl, wavelength=spdc_hc.wavelengths[0], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        photon1_hc = Photon("", tl, wavelength=spdc_hc.wavelengths[1], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_hc = spdc_hc._generate_tmsv_state()
+        keys = [photon0_hc.quantum_state, photon1_hc.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_hc)
+
+        # photon loss upon absorption by memories
+        key_anl = photon1_anl.quantum_state
+        loss_anl = 1 - memo_anl.absorption_efficiency
+        tl.quantum_manager.add_loss(key_anl, loss_anl)
+        key_hc = photon1_hc.quantum_state
+        loss_hc = 1 - memo_hc.absorption_efficiency
+        tl.quantum_manager.add_loss(key_hc, loss_hc)
+
+        # transmission loss through optical fibres
+        key_anl = photon0_anl.quantum_state
+        loss_anl = channel_anl.loss
+        tl.quantum_manager.add_loss(key_anl, loss_anl)
+        key_hc = photon0_hc.quantum_state
+        loss_hc = channel_anl.loss
+        tl.quantum_manager.add_loss(key_hc, loss_hc)
+
+        # QSDetector measurement
+        povms = bsm.povms
+        povm_tuple = tuple([tuple(map(tuple, povm)) for povm in povms])
+        keys = [photon0_anl.quantum_state, photon0_hc.quantum_state]
+        new_state, all_keys = tl.quantum_manager._prepare_state(keys)
+        indices = tuple([all_keys.index(key) for key in keys])
+        state_tuple = tuple(map(tuple, new_state))
+        states, probs = measure_multiple_with_cache_fock_density(state_tuple, indices, len(all_keys), povm_tuple,
+                                                                 tl.quantum_manager.truncation)
+        state_plus, state_minus = states[1], states[2]
+    
+        # calculate remaining state
+        indices = tuple([all_keys.index(key) for key in keys])
+        new_state_tuple = tuple(map(tuple, state_plus))
+        remaining_state = density_partial_trace(new_state_tuple, indices, len(all_keys),
+                                                tl.quantum_manager.truncation)
+        remaining_keys = [key for key in all_keys if key not in keys]
+            
+        # effective Bell state generated 
+        def effective_state(state):
+            state_copy = copy(state)
+            state_copy[0][0] = 0
+            state_copy = state_copy/np.trace(state_copy)
+            
+            return state_copy
+        
+        remaining_state_eff = effective_state(remaining_state)
+        
+        # calculate the fidelity with reference Bell state
+        bell_plus = build_bell_state(tl.quantum_manager.truncation, "plus")
+        bell_minus = build_bell_state(tl.quantum_manager.truncation, "minus")
+        fidelity = np.trace(remaining_state_eff.dot(bell_minus)).real
+
+        print("Directly calculated effective fidelity:", fidelity)
+
+        """Pre-simulation explicit calculation of entanglement generation rate based on calculation above"""
+
+        if calculate_rate_direct:
+            duration_photon = MODE_NUM / SPDC_FREQUENCY  # duration of emitted photon train from SPDC source
+            delay_fiber_anl = anl.qchannels[erc_name].delay
+            delay_fiber_hc = hc.qchannels[erc_name].delay
+            assert delay_fiber_anl == delay_fiber_hc
+            delay_fiber = delay_fiber_anl  # time for photon to travel from SPDC source to BSM device
+            delay_classical = DELAY_CLASSICAL  # delay for classical communication between BSM node and memory nodes
+
+            # total duration from first photon emitted to last photon's detection result communicated back
+            duration_tot = duration_photon + delay_fiber + delay_classical
+
+            prob_herald = probs[1] + probs[2]  # calculate heralding probability
+            num_generated_avg = MODE_NUM * prob_herald  # average number of entangled pairs generated in one emission cycle
+
+            rate = num_generated_avg / duration_tot
+
+            print("Directly calculated entanglement generation:", rate)
+
+    """Run Simulation"""
 
     for i in range(num_direct_trials):
         # start protocol for emitting
         process = Process(anl.emit_protocol, "start", [])
-        event = Event(time_anl, process)
+        event = Event(start_time_anl, process)
         tl.schedule(event)
         process = Process(hc.emit_protocol, "start", [])
-        event = Event(time_hc, process)
+        event = Event(start_time_hc, process)
         tl.schedule(event)
 
         tl.run()
         print("finished direct measurement trial {} out of {}".format(i+1, num_direct_trials))
 
         # collect data
-        bsm_res = erc.get_valid_bins(start_time_bsm, MODE_NUM, FREQUENCY)
-        meas_res = erc_2.get_detector_entries(erc_2.direct_detector_name, start_time_meas, MODE_NUM, FREQUENCY)
-        num_bsm_res = sum(bsm_res)
-        meas_res_valid = [m for m, b in zip(meas_res, bsm_res) if b]
-        probs = [0.0] * 4
+
+        # BSM results determine relative sign of reference Bell state and herald successful entanglement
+        bsm_res = erc.get_detector_entries(erc.bsm_name, start_time_bsm, MODE_NUM, SPDC_FREQUENCY)
+        bsm_success_indices = [i for i, res in enumerate(bsm_res) if res == 1 or res == 2]
+        meas_res = erc_2.get_detector_entries(erc_2.direct_detector_name, start_time_meas, MODE_NUM, SPDC_FREQUENCY)
+
+        num_bsm_res = len(bsm_success_indices)
+        meas_res_valid = [meas_res[i] for i in bsm_success_indices]
+        counts_diag = [0] * 4
         for j in range(4):
-            probs[j] = meas_res_valid.count(j) / num_bsm_res
-        results_direct_measurement.append(probs)
+            counts_diag[j] = meas_res_valid.count(j)
+        res_diag = {"counts": counts_diag, "total_count": num_bsm_res}
+        results_direct_measurement.append(res_diag)
 
         # reset timeline
         tl.time = 0
@@ -282,35 +646,59 @@ def get_detector_entries(self, detector_name: str, start_time: int, num_bins: in
     for i, phase in enumerate(phase_settings):
         erc_2.set_phase(phase)
 
-        # start protocol for emitting
-        process = Process(anl.emit_protocol, "start", [])
-        event = Event(time_anl, process)
-        tl.schedule(event)
-        process = Process(hc.emit_protocol, "start", [])
-        event = Event(time_hc, process)
-        tl.schedule(event)
-
-        tl.run()
-        print("finished interference measurement trial {} out of {}".format(i+1, len(phase_settings)))
-
-        # collect data
-        bsm_res = erc.get_valid_bins(start_time_bsm, MODE_NUM, FREQUENCY)
-        meas_res = erc_2.get_detector_entries(erc_2.bs_detector_name, start_time_meas, MODE_NUM, FREQUENCY)
-        meas_res_valid = [m for m, b in zip(meas_res, bsm_res) if b]
-        num_detector_0 = meas_res.count(1) + meas_res_valid.count(3)
-        num_detector_1 = meas_res.count(2) + meas_res_valid.count(3)
-        freqs = [num_detector_0/(total_time * 1e-12), num_detector_1/(total_time * 1e-12)]
-        results_bs_measurement.append(freqs)
-
-        # reset timeline
-        tl.time = 0
-        tl.init()
+        for j in range(num_bs_trials_per_phase):
+            # start protocol for emitting
+            process = Process(anl.emit_protocol, "start", [])
+            event = Event(start_time_anl, process)
+            tl.schedule(event)
+            process = Process(hc.emit_protocol, "start", [])
+            event = Event(start_time_hc, process)
+            tl.schedule(event)
+
+            tl.run()
+            print("finished interference measurement trial {} out of {} for phase {} out ouf {}".format(
+                j+1, num_bs_trials_per_phase, i+1, len(phase_settings)))
+
+            # collect data
+
+            # relative sign should influence interference pattern
+            bsm_res = erc.get_detector_entries(erc.bsm_name, start_time_bsm, MODE_NUM, SPDC_FREQUENCY)
+            bsm_success_indices_1 = [i for i, res in enumerate(bsm_res) if res == 1]
+            bsm_success_indices_2 = [i for i, res in enumerate(bsm_res) if res == 2]
+            meas_res = erc_2.get_detector_entries(erc_2.bs_detector_name, start_time_meas, MODE_NUM, SPDC_FREQUENCY)
+            res_interference = {}
+
+            # detector 1
+            num_bsm_res = len(bsm_success_indices_1)
+            meas_res_valid = [meas_res[i] for i in bsm_success_indices_1]
+            num_detector_0 = meas_res_valid.count(1) + meas_res_valid.count(3)
+            num_detector_1 = meas_res_valid.count(2) + meas_res_valid.count(3)
+            counts_interfere = [num_detector_0, num_detector_1]
+            res_interference["counts1"] = counts_interfere
+            res_interference["total_count1"] = num_bsm_res
+
+            # detector 2
+            num_bsm_res = len(bsm_success_indices_2)
+            meas_res_valid = [meas_res[i] for i in bsm_success_indices_2]
+            num_detector_0 = meas_res_valid.count(1) + meas_res_valid.count(3)
+            num_detector_1 = meas_res_valid.count(2) + meas_res_valid.count(3)
+            counts_interfere = [num_detector_0, num_detector_1]
+            res_interference["counts2"] = counts_interfere
+            res_interference["total_count2"] = num_bsm_res
+
+            results_bs_measurement[i].append(res_interference)
+
+            # reset timeline
+            tl.time = 0
+            tl.init()
 
     """Store results"""
 
-    info = {"direct results": results_direct_measurement, "bs results": results_bs_measurement}
+    # open file to store experiment results
+    Path("results").mkdir(parents=True, exist_ok=True)
+    filename = "results/absorptive.json"
+    fh = open(filename, 'w')
+    info = {"num_direct_trials": num_direct_trials, "num_bs_trials": num_bs_trials_per_phase,
+            "num_phase": len(phase_settings),
+            "direct results": results_direct_measurement, "bs results": results_bs_measurement}
     dump(info, fh)
-
-    # plt.bar(list(range(4)), counts)
-    # plt.yscale('log')
-    # plt.show()
diff --git a/example/absorptive_fidelity.py b/example/absorptive_fidelity.py
new file mode 100644
index 00000000..eef52e23
--- /dev/null
+++ b/example/absorptive_fidelity.py
@@ -0,0 +1,151 @@
+from absorptive_experiment import *
+from json import JSONEncoder
+from matplotlib import pyplot as plt
+
+dim = 100
+mean_num1_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+mean_num2_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+fidelity_mat = np.zeros((dim, dim))
+
+
+# effective Bell state generated
+def effective_state(state):
+    state[0][0] = 0
+    state = state / np.trace(state)
+    return state
+
+
+# for storage
+class NumpyEncoder(JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        return JSONEncoder.default(self, obj)
+
+
+for idx1, mean_num1 in enumerate(mean_num1_list):
+    for idx2, mean_num2 in enumerate(mean_num2_list):
+
+        tl = Timeline(time, formalism=FOCK_DENSITY_MATRIX_FORMALISM, truncation=TRUNCATION)
+
+        anl_name = "Argonne"
+        hc_name = "Harper Court"
+        erc_name = "Eckhardt Research Center BSM"
+        erc_2_name = "Eckhardt Research Center Measurement"
+        seeds = [1, 2, 3, 4]
+        src_list = [anl_name, hc_name]  # the list of sources, note the order
+
+        anl = EndNode(anl_name, tl, hc_name, erc_name, erc_2_name, mean_photon_num=mean_num1,
+                      spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY1, abs_effi=ABS_EFFICIENCY1,
+                      afc_efficiency=efficiency1, mode_number=MODE_NUM)
+        hc = EndNode(hc_name, tl, anl_name, erc_name, erc_2_name, mean_photon_num=mean_num2,
+                     spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY2, abs_effi=ABS_EFFICIENCY2,
+                     afc_efficiency=efficiency2, mode_number=MODE_NUM)
+        erc = EntangleNode(erc_name, tl, src_list)
+        erc_2 = MeasureNode(erc_2_name, tl, src_list)
+
+        for seed, node in zip(seeds, [anl, hc, erc, erc_2]):
+            node.set_seed(seed)
+
+        # extend fiber lengths to be equivalent
+        fiber_length = max(DIST_ANL_ERC, DIST_HC_ERC)
+
+        qc1 = add_channel(anl, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc2 = add_channel(hc, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc3 = add_channel(anl, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc4 = add_channel(hc, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+
+        tl.init()
+
+        # Pre-simulation explicit calculation of entanglement fidelity upon successful BSM
+
+        # use non-transmitted Photon as interface with existing methods in SeQUeNCe
+        spdc_anl = anl.components[anl.spdc_name]
+        spdc_hc = hc.components[hc.spdc_name]
+        memo_anl = anl.components[anl.memo_name]
+        memo_hc = hc.components[hc.memo_name]
+        channel_anl = anl.qchannels[erc_name]
+        channel_hc = hc.qchannels[erc_name]
+        bsm = erc.components[erc.bsm_name]
+
+        # photon0: idler, photon1: signal
+        photon0_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[0], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        photon1_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[1], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_anl = spdc_anl._generate_tmsv_state()
+        keys = [photon0_anl.quantum_state, photon1_anl.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_anl)
+
+        photon0_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[0], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        photon1_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[1], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_hc = spdc_hc._generate_tmsv_state()
+        keys = [photon0_hc.quantum_state, photon1_hc.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_hc)
+
+        # photon loss upon absorption by memories
+        key_anl_memo = photon1_anl.quantum_state
+        loss_anl_memo = 1 - memo_anl.absorption_efficiency
+        tl.quantum_manager.add_loss(key_anl_memo, loss_anl_memo)
+        key_hc_memo = photon1_hc.quantum_state
+        loss_hc_memo = 1 - memo_hc.absorption_efficiency
+        tl.quantum_manager.add_loss(key_hc_memo, loss_hc_memo)
+
+        # transmission loss through optical fibres (overwrites previous variables)
+        key_anl_pho = photon0_anl.quantum_state
+        loss_anl_pho = channel_anl.loss
+        tl.quantum_manager.add_loss(key_anl_pho, loss_anl_pho)
+        key_hc_pho = photon0_hc.quantum_state
+        loss_hc_pho = channel_anl.loss
+        tl.quantum_manager.add_loss(key_hc_pho, loss_hc_pho)
+
+        # QSDetector measurement and remaining state after partial trace
+        povms = bsm.povms
+        povm_tuple = tuple([tuple(map(tuple, povm)) for povm in povms])
+        keys = [photon0_anl.quantum_state, photon0_hc.quantum_state]
+        new_state, all_keys = tl.quantum_manager._prepare_state(keys)
+        indices = tuple([all_keys.index(key) for key in keys])
+        state_tuple = tuple(map(tuple, new_state))
+        states, probs = measure_multiple_with_cache_fock_density(state_tuple, indices, len(all_keys), povm_tuple,
+                                                                 tl.quantum_manager.truncation)
+        state_plus, state_minus = states[1], states[2]
+
+        indices = tuple([all_keys.index(key) for key in keys])
+        new_state_tuple = tuple(map(tuple, state_plus))
+        remaining_state = density_partial_trace(new_state_tuple, indices, len(all_keys),
+                                                tl.quantum_manager.truncation)
+
+
+        remaining_state_eff = effective_state(remaining_state)
+
+        # calculate the fidelity with reference Bell state
+        bell_plus = build_bell_state(tl.quantum_manager.truncation, "plus")
+        bell_minus = build_bell_state(tl.quantum_manager.truncation, "minus")
+        fidelity = np.trace(remaining_state_eff.dot(bell_minus)).real
+
+        fidelity_mat[idx1][idx2] = fidelity
+
+    print("finished mean num {} out of {}".format(idx1 + 1, dim))
+
+"""Store results"""
+
+# open file to store experiment results
+Path("results").mkdir(parents=True, exist_ok=True)
+filename = "results/absorptive_fidelity.json"
+fh = open(filename, 'w')
+info = {"fidelity matrix": fidelity_mat}
+dump(info, fh, cls=NumpyEncoder)
+
+# plot the 2-d matrix
+num1, num2 = np.meshgrid(mean_num1_list, mean_num2_list)
+fig = plt.figure(figsize=(8, 6))
+im = plt.imshow(fidelity_mat, cmap=plt.cm.RdBu, extent=(0.02, 0.1, 0.1, 0.02))
+plt.colorbar(im)
+plt.title("Effective Fidelity with Mean Photon Number")
+plt.xlabel("mean photon number of source 1 $\mu_1$")
+plt.ylabel("mean photon number of source 2 $\mu_2$")
+plt.show()
diff --git a/example/absorptive_fidelity_graph.py b/example/absorptive_fidelity_graph.py
new file mode 100644
index 00000000..420e7f20
--- /dev/null
+++ b/example/absorptive_fidelity_graph.py
@@ -0,0 +1,41 @@
+from json5 import load
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+filename = "results/absorptive_fidelity.json"
+data = load(open(filename))
+fidelity_mat = np.asarray(data["fidelity matrix"])
+
+dim = 10
+mean_num1_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+mean_num2_list = np.linspace(0.02, 0.1, num=dim, endpoint=True)
+num1, num2 = np.meshgrid(mean_num1_list, mean_num2_list)
+
+# plot the 2-d matrix
+plt.rc('font', size=10)
+plt.rc('axes', titlesize=18)
+plt.rc('axes', labelsize=15)
+plt.rc('xtick', labelsize=15)
+plt.rc('ytick', labelsize=15)
+plt.rc('legend', fontsize=15)
+
+plt.rcParams['axes.titley'] = 1.05
+plt.rcParams['axes.titlepad'] = 0
+
+fig = plt.figure(figsize=(8, 6))
+ax = fig.subplots()
+divider = make_axes_locatable(ax)
+cax = divider.append_axes('right', size='5%', pad=0.1)
+
+im = ax.imshow(fidelity_mat, cmap=plt.cm.RdBu, extent=(0.02, 0.1, 0.1, 0.02))
+fig.colorbar(im, cax=cax, orientation='vertical')
+ax.invert_yaxis()
+ax.set_title("Effective Fidelity versus Mean Photon Number")
+ax.set_xlabel(r'Mean Photon Number of Source 1 $\mu_1$')
+ax.set_ylabel(r'Mean Photon Number of Source 2 $\mu_2$')
+ax.yaxis.set_ticks(np.linspace(0.02, 0.1, 5))
+fig.tight_layout()
+plt.savefig('fidelity.png')
+plt.show()
diff --git a/example/absorptive_graph.py b/example/absorptive_graph.py
index d673c080..8834d63f 100644
--- a/example/absorptive_graph.py
+++ b/example/absorptive_graph.py
@@ -1,5 +1,9 @@
 from json import load
+import numpy as np
+from scipy.optimize import curve_fit
 import matplotlib.pyplot as plt
+import matplotlib.ticker as tck
+from matplotlib.colors import LightSource
 
 
 filename = "results/absorptive.json"
@@ -8,9 +12,105 @@
 direct_results = data["direct results"]
 bs_results = data["bs results"]
 
-# plotting bs results
-freq_0 = [res[0] for res in bs_results]
-freq_1 = [res[1] for res in bs_results]
-plt.plot(freq_0)
-plt.plot(freq_1)
+# acquire direct detection results
+bins = np.zeros(4)
+for trial in direct_results:
+    bins += np.array(trial["counts"])
+bins *= (1 / sum(bins))
+
+# acquire bs results
+num = data["num_phase"]
+phases = np.linspace(0, 2*np.pi, num=num)
+freq_0 = []
+freq_1 = []
+st_dev_0 = []
+st_dev_1 = []
+for res in bs_results:
+    counts = []
+    total = 0
+    for trial in res:
+        counts.append(trial["counts1"])
+        total += trial["total_count1"]
+    counts = np.array(counts, dtype=float)
+    st_dev = np.std(counts, axis=0)
+    # average
+    counts = np.sum(counts, axis=0)
+    counts *= (1 / total)
+    st_dev *= (1 / total)
+    freq_0.append(counts[0])
+    freq_1.append(counts[1])
+    st_dev_0.append(2 * st_dev[0])  # mult. by 2, so have 1 on top and bottom of point
+    st_dev_1.append(2 * st_dev[1])
+
+# curve fitting
+def sin(x, A, phi, z): return A*np.sin(x + phi) + z
+params_0, _ = curve_fit(sin, phases, freq_0)
+params_1, _ = curve_fit(sin, phases, freq_1)
+vis_0 = abs(params_0[0]) / abs(params_0[2])
+vis_1 = abs(params_1[0]) / abs(params_1[2])
+
+# off-diagonal density matrix
+vis_total = (vis_0 + vis_1) / 2
+off_diag = vis_total * (bins[1] + bins[2]) / 2
+
+# plotting density matrix
+# plot the 2-d matrix
+plt.rc('font', size=10)
+plt.rc('axes', titlesize=18)
+plt.rc('axes', labelsize=15)
+plt.rc('xtick', labelsize=15)
+plt.rc('ytick', labelsize=15)
+plt.rc('legend', fontsize=15)
+
+fig = plt.figure(figsize=(8, 8))
+ax = fig.add_subplot(projection="3d")
+ax.view_init(azim=-30, elev=20)
+
+offset = 3
+total = list(bins) + [off_diag, off_diag]
+height = np.log10(total) + offset
+width = depth = 0.8
+x = np.array([0, 1, 2, 3, 1, 2]) - (width/2)
+y = np.array([0, 1, 2, 3, 2, 1]) - (width/2)
+color = (['tab:blue'] * 4) + (['tab:orange'] * 2)
+ls = LightSource(100, 60)
+ax.bar3d(x, y, np.zeros(6), width, depth, height,
+         color=color, edgecolor='black', shade=True, lightsource=ls)
+
+def log_tick_formatter(val, pos=None): return f"$10^{{{int(val-offset)}}}$"
+ax.zaxis.set_major_formatter(tck.FuncFormatter(log_tick_formatter))
+ax.zaxis.set_major_locator(tck.MaxNLocator(integer=True))
+
+# ax.set_title(r'Reconstructed $|\tilde{\rho}|$')
+ax.set_xticks([0, 1, 2, 3])
+ax.set_yticks([0, 1, 2, 3])
+ax.set_xticklabels([r'$|00\rangle$', r'$|01\rangle$', r'$|10\rangle$', r'$|11\rangle$'])
+ax.set_yticklabels([r'$\langle00|$', r'$\langle01|$', r'$\langle10|$', r'$\langle11|$'])
+plt.savefig('density.png', bbox_inches='tight')
+plt.show()
+
+# plotting interference
+fig = plt.figure(figsize=(8, 5))
+ax = fig.add_subplot()
+
+ax.errorbar(phases/np.pi, freq_0, yerr=st_dev_0, marker='o', ls='', capsize=5, label=r'$p_{01}$')
+ax.errorbar(phases/np.pi, freq_1, yerr=st_dev_1, marker='o', ls='', capsize=5, label=r'$p_{10}$')
+ax.plot(phases/np.pi, np.zeros(phases.shape), ls='--', color='gray')
+plt.plot(phases/np.pi, sin(phases, *params_0), color='tab:blue', ls='--')
+plt.plot(phases/np.pi, sin(phases, *params_1), color='tab:orange', ls='--')
+
+ax.set_title("Photon Interference")
+ax.xaxis.set_major_formatter(tck.FormatStrFormatter(r'%g $\pi$'))
+ax.xaxis.set_major_locator(tck.MultipleLocator(base=0.5))
+ax.set_xlabel("Relative Phase")
+ax.set_ylabel("Detection Rate")
+ax.legend()
+fig.tight_layout()
+plt.savefig('interference.png')
 plt.show()
+
+# output
+print("Detector 1 visibility:", vis_1)
+print("Detector 2 visibility:", vis_0)
+print("Combined:", vis_total)
+print("Off-diagonal:", off_diag)
diff --git a/example/absorptive_rate.py b/example/absorptive_rate.py
new file mode 100644
index 00000000..9fd17d09
--- /dev/null
+++ b/example/absorptive_rate.py
@@ -0,0 +1,139 @@
+from absorptive_experiment import *
+from matplotlib import pyplot as plt
+
+num = 4
+mean_num_list = np.linspace(0.02, 0.05, num=num, endpoint=True)
+mode_num_list = np.arange(10, 105, step=5)
+rate_list = []
+
+for mean_num in mean_num_list:
+    mean_num_list = []
+
+    for mode_num in mode_num_list:
+        tl = Timeline(time, formalism=FOCK_DENSITY_MATRIX_FORMALISM, truncation=TRUNCATION)
+
+        anl_name = "Argonne"
+        hc_name = "Harper Court"
+        erc_name = "Eckhardt Research Center BSM"
+        erc_2_name = "Eckhardt Research Center Measurement"
+        seeds = [1, 2, 3, 4]
+        src_list = [anl_name, hc_name]  # the list of sources, note the order
+
+        anl = EndNode(anl_name, tl, hc_name, erc_name, erc_2_name, mean_photon_num=mean_num,
+                      spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY1, abs_effi=ABS_EFFICIENCY1,
+                      afc_efficiency=efficiency1, mode_number=mode_num)
+        hc = EndNode(hc_name, tl, anl_name, erc_name, erc_2_name, mean_photon_num=mean_num,
+                     spdc_frequency=SPDC_FREQUENCY, memo_frequency=MEMO_FREQUENCY2, abs_effi=ABS_EFFICIENCY2,
+                     afc_efficiency=efficiency2, mode_number=mode_num)
+        erc = EntangleNode(erc_name, tl, src_list)
+        erc_2 = MeasureNode(erc_2_name, tl, src_list)
+
+        for seed, node in zip(seeds, [anl, hc, erc, erc_2]):
+            node.set_seed(seed)
+
+        # extend fiber lengths to be equivalent
+        fiber_length = max(DIST_ANL_ERC, DIST_HC_ERC)
+
+        qc1 = add_channel(anl, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc2 = add_channel(hc, erc, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc3 = add_channel(anl, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+        qc4 = add_channel(hc, erc_2, tl, distance=fiber_length, attenuation=ATTENUATION, frequency=10e7)
+
+        tl.init()
+
+        # Pre-simulation explicit calculation of entanglement fidelity upon successful BSM
+
+        # use non-transmitted Photon as interface with existing methods in SeQUeNCe
+        spdc_anl = anl.components[anl.spdc_name]
+        spdc_hc = hc.components[hc.spdc_name]
+        memo_anl = anl.components[anl.memo_name]
+        memo_hc = hc.components[hc.memo_name]
+        channel_anl = anl.qchannels[erc_name]
+        channel_hc = hc.qchannels[erc_name]
+        bsm = erc.components[erc.bsm_name]
+
+        # photon0: idler, photon1: signal
+        photon0_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[0], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        photon1_anl = Photon("", spdc_anl.timeline, wavelength=spdc_anl.wavelengths[1], location=spdc_anl,
+                             encoding_type=spdc_anl.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_anl = spdc_anl._generate_tmsv_state()
+        keys = [photon0_anl.quantum_state, photon1_anl.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_anl)
+
+        photon0_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[0], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        photon1_hc = Photon("", spdc_hc.timeline, wavelength=spdc_hc.wavelengths[1], location=spdc_hc,
+                            encoding_type=spdc_hc.encoding_type, use_qm=True)
+        # set shared state to squeezed state
+        state_spdc_hc = spdc_hc._generate_tmsv_state()
+        keys = [photon0_hc.quantum_state, photon1_hc.quantum_state]
+        tl.quantum_manager.set(keys, state_spdc_hc)
+
+        # photon loss upon absorption by memories
+        key_anl_memo = photon1_anl.quantum_state
+        loss_anl_memo = 1 - memo_anl.absorption_efficiency
+        tl.quantum_manager.add_loss(key_anl_memo, loss_anl_memo)
+        key_hc_memo = photon1_hc.quantum_state
+        loss_hc_memo = 1 - memo_hc.absorption_efficiency
+        tl.quantum_manager.add_loss(key_hc_memo, loss_hc_memo)
+
+        # transmission loss through optical fibres (overwrites previous variables)
+        key_anl_pho = photon0_anl.quantum_state
+        loss_anl_pho = channel_anl.loss
+        tl.quantum_manager.add_loss(key_anl_pho, loss_anl_pho)
+        key_hc_pho = photon0_hc.quantum_state
+        loss_hc_pho = channel_anl.loss
+        tl.quantum_manager.add_loss(key_hc_pho, loss_hc_pho)
+
+        # QSDetector measurement and remaining state after partial trace
+        povms = bsm.povms
+        povm_tuple = tuple([tuple(map(tuple, povm)) for povm in povms])
+        keys = [photon0_anl.quantum_state, photon0_hc.quantum_state]
+        new_state, all_keys = tl.quantum_manager._prepare_state(keys)
+        indices = tuple([all_keys.index(key) for key in keys])
+        state_tuple = tuple(map(tuple, new_state))
+        states, probs = measure_multiple_with_cache_fock_density(state_tuple, indices, len(all_keys), povm_tuple,
+                                                                 tl.quantum_manager.truncation)
+
+        # Pre-simulation explicit calculation of entanglement generation rate based on calculation above
+
+        duration_photon = mode_num / SPDC_FREQUENCY * 1e12  # duration of emitted photon train from SPDC source
+        delay_fiber_anl = anl.qchannels[erc_name].delay
+        delay_fiber_hc = hc.qchannels[erc_name].delay
+        assert delay_fiber_anl == delay_fiber_hc
+        delay_fiber = delay_fiber_anl  # time for photon to travel from SPDC source to BSM device
+        delay_classical = DELAY_CLASSICAL * 1e12  # delay for classical communication between BSM node and memory nodes
+
+        # total duration from first photon emitted to last photon's detection result communicated back
+        duration_tot = duration_photon + delay_fiber + delay_classical
+
+        prob_herald = probs[1] + probs[2]  # calculate heralding probability
+        num_generated_avg = mode_num * prob_herald  # average number of entangled pairs generated in one emission cycle
+
+        rate = num_generated_avg / duration_tot * 1e12
+
+        mean_num_list.append(rate)
+    rate_list.append(mean_num_list)
+
+"""Store results"""
+
+# open file to store experiment results
+Path("results").mkdir(parents=True, exist_ok=True)
+filename = "results/absorptive_rate.json"
+fh = open(filename, 'w')
+info = {"rate_list": rate_list}
+dump(info, fh)
+
+list1, list2, list3, list4 = rate_list
+fig = plt.figure(figsize=(8, 6))
+plt.plot(mode_num_list, list1, label='$\mu=0.02$')
+plt.plot(mode_num_list, list2, label='$\mu=0.03$')
+plt.plot(mode_num_list, list3, label='$\mu=0.04$')
+plt.plot(mode_num_list, list4, label='$\mu=0.05$')
+plt.title("Entanglement Generation Rate with Mode Number")
+plt.xlabel("mode number $M$")
+plt.ylabel("entanglement generation rate $R$ (Hz)")
+plt.legend()
+plt.show()
diff --git a/example/absorptive_rate_graph.py b/example/absorptive_rate_graph.py
new file mode 100644
index 00000000..833dc1b7
--- /dev/null
+++ b/example/absorptive_rate_graph.py
@@ -0,0 +1,34 @@
+from json5 import load
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+filename = "results/absorptive_rate.json"
+data = load(open(filename))
+
+list1, list2, list3, list4 = data["rate_list"]
+mode_num_list = np.arange(10, 105, step=5)
+
+plt.rc('font', size=10)
+plt.rc('axes', titlesize=18)
+plt.rc('axes', labelsize=15)
+plt.rc('xtick', labelsize=15)
+plt.rc('ytick', labelsize=15)
+plt.rc('legend', fontsize=15)
+
+fig = plt.figure(figsize=(8, 5))
+ax = fig.subplots()
+
+ax.plot(mode_num_list, list1, label=r'$\mu=0.02$')
+ax.plot(mode_num_list, list2, label=r'$\mu=0.03$')
+ax.plot(mode_num_list, list3, label=r'$\mu=0.04$')
+ax.plot(mode_num_list, list4, label=r'$\mu=0.05$')
+
+ax.set_title("Entanglement Generation Rate with Mode Number")
+ax.set_xlabel("Mode Number $M$")
+ax.set_ylabel("Entanglement Generation Rate $R$ (Hz)")
+ax.legend()
+fig.tight_layout()
+plt.savefig('rate.png')
+plt.show()
+plt.show()
diff --git a/example/starlight_experiments.py b/example/starlight_experiments.py
index d3935373..5f2fc9c6 100644
--- a/example/starlight_experiments.py
+++ b/example/starlight_experiments.py
@@ -8,7 +8,7 @@
 from sequence.topology.router_net_topo import RouterNetTopo
 
 
-def get_component(node: Node, component_type: str):
+def get_component(node: "Node", component_type: str):
     for comp in node.components.values():
         if type(comp).__name__ == component_type:
             return comp
@@ -21,6 +21,7 @@ def get_component(node: Node, component_type: str):
     network_config_file = "example/starlight.json"
     network_topo = RouterNetTopo(network_config_file)
     tl = network_topo.get_timeline()
+    tl.show_progress = True
     routers = network_topo.get_nodes_by_type(RouterNetTopo.QUANTUM_ROUTER)
     bsm_nodes = network_topo.get_nodes_by_type(RouterNetTopo.BSM_NODE)
 
diff --git a/requirements.txt b/requirements.txt
index 1b7afc7a..d212c2e5 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,5 @@
 numpy>=1.22
+scipy
 matplotlib
 dash>=1.20.0
 dash-core-components
diff --git a/setup.py b/setup.py
index 375a2d7d..4e6a50c7 100644
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@
 
 setup(
     name="sequence",
-    version="0.5.4",
+    version="0.6.0",
     author="Xiaoliang Wu, Joaquin Chung, Alexander Kolar, Alexander Kiefer, Eugene Wang, Tian Zhong, Rajkumar Kettimuthu, Martin Suchara",
     author_email="xwu64@hawk.iit.edu, chungmiranda@anl.gov, akolar@anl.gov, akiefer@iu.edu, eugenewang@yahoo.com, tzh@uchicago.edu, kettimut@mcs.anl.gov, msuchara@anl.gov",
     description="Simulator of Quantum Network Communication: SEQUENCE-Python is a prototype version of the official SEQUENCE release.",
diff --git a/src/components/detector.py b/src/components/detector.py
index 4a8453ed..2256a23e 100644
--- a/src/components/detector.py
+++ b/src/components/detector.py
@@ -2,20 +2,23 @@
 
 This module models a single photon detector (SPD) for measurement of individual photons.
 It also defines a QSDetector class, which combines models of different hardware devices to measure photon states in different bases.
-QSDetector is defined as an abstract template and as implementaions for polarization and time bin qubits.
+QSDetector is defined as an abstract template and as implementations for polarization and time bin qubits.
 """
 
 from abc import ABC, abstractmethod
 from typing import TYPE_CHECKING, Any, Dict, List
+from numpy import eye, kron, exp, sqrt
+from scipy.linalg import fractional_matrix_power
+from math import factorial
 
 if TYPE_CHECKING:
     from ..kernel.timeline import Timeline
 
-from ..components.photon import Photon
-from ..components.beam_splitter import BeamSplitter, FockBeamSplitter
-from ..components.switch import Switch
-from ..components.interferometer import Interferometer
-from ..components.circuit import Circuit
+from .photon import Photon
+from .beam_splitter import BeamSplitter
+from .switch import Switch
+from .interferometer import Interferometer
+from .circuit import Circuit
 from ..kernel.entity import Entity
 from ..kernel.event import Event
 from ..kernel.process import Process
@@ -55,7 +58,8 @@ def init(self):
         """Implementation of Entity interface (see base class)."""
         self.next_detection_time = -1
         self.photon_counter = 0
-        self.add_dark_count()
+        if self.dark_count > 0:
+            self.add_dark_count()
 
     def get(self, photon=None, **kwargs) -> None:
         """Method to receive a photon for measurement.
@@ -63,31 +67,23 @@ def get(self, photon=None, **kwargs) -> None:
         Args:
             photon (Photon): photon to detect (currently unused)
 
-        Keyword Args:
-            dark_get (bool): Signifies if the call is the result of a false positive dark count event.
-                If true, will ignore probability calculations (default false).
-
         Side Effects:
             May notify upper entities of a detection event.
         """
 
-        dark_get = kwargs.get("dark_get", False)
-
         self.photon_counter += 1
-        now = self.timeline.now()
-        time = round(now / self.time_resolution) * self.time_resolution
 
         # if get a photon and it has single_atom encoding, measure
         if photon and photon.encoding_type["name"] == "single_atom":
             key = photon.quantum_state
             res = self.timeline.quantum_manager.run_circuit(Detector._meas_circuit, [key],
                                                             self.get_generator().random())
+            # if we measure |0>, return (do not record detection)
             if not res[key]:
                 return
 
-        if (self.get_generator().random() < self.efficiency or dark_get) and now > self.next_detection_time:
-            self.notify({'time': time})
-            self.next_detection_time = now + (1e12 / self.count_rate)  # period in ps
+        if self.get_generator().random() < self.efficiency:
+            self.record_detection()
 
     def add_dark_count(self) -> None:
         """Method to schedule false positive detection events.
@@ -99,20 +95,35 @@ def add_dark_count(self) -> None:
             May schedule future calls to self.
         """
 
-        if self.dark_count > 0:
-            time_to_next = int(self.get_generator().exponential(
+        assert self.dark_count > 0, "Detector().add_dark_count called with 0 dark count rate"
+        time_to_next = int(self.get_generator().exponential(
                 1 / self.dark_count) * 1e12)  # time to next dark count
-            time = time_to_next + self.timeline.now()  # time of next dark count
+        time = time_to_next + self.timeline.now()  # time of next dark count
+
+        process1 = Process(self, "add_dark_count", [])  # schedule photon detection and dark count add in future
+        process2 = Process(self, "record_detection", [])
+        event1 = Event(time, process1)
+        event2 = Event(time, process2)
+        self.timeline.schedule(event1)
+        self.timeline.schedule(event2)
 
-            process1 = Process(self, "add_dark_count", [])  # schedule photon detection and dark count add in future
-            process2 = Process(self, "get", [None], {"dark_get": True})
-            event1 = Event(time, process1)
-            event2 = Event(time, process2)
-            self.timeline.schedule(event1)
-            self.timeline.schedule(event2)
+    def record_detection(self):
+        """Method to record a detection event.
+
+        Will calculate if detection succeeds (by checking if we have passed `next_detection_time`)
+        and will notify observers with the detection time (rounded to the nearest multiple of detection frequency).
+        """
+
+        now = self.timeline.now()
+
+        if now > self.next_detection_time:
+            time = round(now / self.time_resolution) * self.time_resolution
+            self.notify({'time': time})
+            self.next_detection_time = now + (1e12 / self.count_rate)  # period in ps
 
     def notify(self, info: Dict[str, Any]):
         """Custom notify function (calls `trigger` method)."""
+
         for observer in self._observers:
             observer.trigger(self, info)
 
@@ -125,7 +136,7 @@ class QSDetector(Entity, ABC):
     Attributes:
         name (str): label for QSDetector instance.
         timeline (Timeline): timeline for simulation.
-        componennts (List[entity]): list of all aggregated hardware components.
+        components (List[entity]): list of all aggregated hardware components.
         detectors (List[Detector]): list of attached detectors.
         trigger_times (List[List[int]]): tracks simulation time of detection events for each detector.
     """
@@ -155,6 +166,21 @@ def trigger(self, detector: Detector, info: Dict[str, Any]) -> None:
         detector_index = self.detectors.index(detector)
         self.trigger_times[detector_index].append(info['time'])
 
+    def set_detector(self, idx: int,  efficiency=0.9, dark_count=0, count_rate=int(25e6), time_resolution=150):
+        """Method to set the properties of an attached detector.
+
+        Args:
+            idx (int): the index of attached detector whose properties are going to be set.
+            For other parameters see the `Detector` class. Default values are same as in `Detector` class.
+        """
+        assert 0 <= idx < len(self.detectors), "`idx` must be a valid index of attached detector."
+
+        detector = self.detectors[idx]
+        detector.efficiency = efficiency
+        detector.dark_count = dark_count
+        detector.count_rate = count_rate
+        detector.time_resolution = time_resolution
+
     def get_photon_times(self):
         return self.trigger_times
 
@@ -285,14 +311,15 @@ def update_interferometer_params(self, arg_name: str, value: Any) -> None:
 class QSDetectorFockDirect(QSDetector):
     """QSDetector to directly measure photons in Fock state.
 
-    Used to calculate diagonal density matrix elements.
+    Usage: to measure diagonal elements of effective density matrix.
 
     Attributes:
         name (str): label for QSDetector instance.
         timeline (Timeline): timeline for simulation.
-        src_list (List[int]):
+        src_list (List[str]): list of two sources which send photons to this detector (length 2).
         detectors (List[Detector]): list of attached detectors (length 2).
         trigger_times (List[List[int]]): tracks simulation time of detection events for each detector.
+        arrival_times (List[List[int]]): tracks simulation time of Photon arrival at each input port
     """
 
     def __init__(self, name: str, timeline: "Timeline", src_list: List[str]):
@@ -303,18 +330,64 @@ def __init__(self, name: str, timeline: "Timeline", src_list: List[str]):
         for i in range(2):
             d = Detector(name + ".detector" + str(i), timeline)
             self.detectors.append(d)
+        self.components = self.detectors
 
         self.trigger_times = [[], []]
+        self.arrival_times = [[], []]
 
-        self.components = self.detectors
+        self.povms = [None] * 4
+
+    def init(self):
+        self._generate_povms()
+        super().init()
+
+    def _generate_povms(self):
+        """Method to generate POVM operators corresponding to photon detector having 0 and 1 click
+        Will be used to generated outcome probability distribution.
+        """
+
+        # assume using Fock quantum manager
+        truncation = self.timeline.quantum_manager.truncation
+        create, destroy = self.timeline.quantum_manager.build_ladder()
+
+        create0 = create * sqrt(self.detectors[0].efficiency)
+        destroy0 = destroy * sqrt(self.detectors[0].efficiency)
+        series_elem_list = [((-1)**i) * fractional_matrix_power(create0, i+1).dot(
+            fractional_matrix_power(destroy0, i+1)) / factorial(i+1) for i in range(truncation)]
+        povm0_1 = sum(series_elem_list)
+        povm0_0 = eye(truncation+1) - povm0_1
+
+        create1 = create * sqrt(self.detectors[1].efficiency)
+        destroy1 = destroy * sqrt(self.detectors[1].efficiency)
+        series_elem_list = [((-1) ** i) * fractional_matrix_power(create1, i + 1).dot(
+            fractional_matrix_power(destroy1, i + 1)) / factorial(i + 1) for i in range(truncation)]
+        povm1_1 = sum(series_elem_list)
+        povm1_0 = eye(truncation + 1) - povm0_1
+
+        self.povms = [povm0_0, povm0_1, povm1_0, povm1_1]
 
     def get(self, photon: "Photon", **kwargs):
         src = kwargs["src"]
-        # measure (0/1 determines existence of photons in encoding)
-        res = Photon.measure(None, photon, self.get_generator())
-        if res:
-            detector_num = self.src_list.index(src)
-            self.detectors[detector_num].get()
+        assert photon.encoding_type["name"] == "fock", "Photon must be in Fock representation."
+        input_port = self.src_list.index(src)  # determine at which input the Photon arrives, an index
+
+        # record arrival time
+        arrival_time = self.timeline.now()
+        self.arrival_times[input_port].append(arrival_time)
+
+        key = photon.quantum_state  # the photon's key pointing to the quantum state in quantum manager
+        samp = self.get_generator().random()  # random measurement sample
+        if input_port == 0:
+            result = self.timeline.quantum_manager.measure([key], self.povms[0:2], samp)
+        elif input_port == 1:
+            result = self.timeline.quantum_manager.measure([key], self.povms[2:4], samp)
+        else:
+            raise Exception("too many input ports for QSDFockDirect {}".format(self.name))
+
+        assert result in list(range(len(self.povms))), "The measurement outcome is not valid."
+        if result == 1:
+            # trigger time recording will be done by SPD
+            self.detectors[input_port].record_detection()
 
     def get_photon_times(self) -> List[List[int]]:
         trigger_times = self.trigger_times
@@ -327,25 +400,183 @@ def set_basis_list(self, basis_list: List[int], start_time: int, frequency: int)
 
 
 class QSDetectorFockInterference(QSDetector):
-    """WIP"""
+    """QSDetector with two input ports and two photon detectors behind beamsplitter.
 
-    def __init__(self, name: str, timeline: "Timeline"):
+    The detectors will physically measure the two beamsplitter output  photonic modes' Fock states, respectively.
+    POVM operators which apply to pre-beamsplitter photonic state are used.
+    NOTE: in the current implementation, to realize interference, we require that Photons arrive at both input ports
+    simultaneously, and at most 1 Photon instance can be input at an input port at a time.
+
+    Usage: to realize Bell state measurement (BSM) and to measure off-diagonal elements of the effective density matrix.
+
+    Attributes:
+        name (str): label for QSDetector instance.
+        timeline (Timeline): timeline for simulation.
+        src_list (List[str]): list of two sources which send photons to this detector (length 2).
+        detectors (List[Detector]): list of attached detectors (length 2).
+        phase (float): relative phase between two input optical paths.
+        trigger_times (List[List[int]]): tracks simulation time of detection events for each detector.
+        detect_info (List[List[Dict]]): tracks detection information, including simulation time of detection events
+            and detection outcome for each detector.
+        arrival_times (List[List[int]]): tracks simulation time of arrival of photons at each input mode.
+
+        temporary_photon_info (List[Dict]): temporary list of information of Photon arriving at each input port.
+            Specific to current implementation. At most 1 Photon's information will be recorded in a dictionary.
+            When there are 2 non-empty dictionaries,
+            e.g. [{"photon":Photon1, "time":arrival_time1}, {"photon":Photon2, "time":arrival_time2}],
+            the entangling measurement will be carried out. After measurement, the temporary list will be reset.
+    """
+
+    def __init__(self, name: str, timeline: "Timeline", src_list: List[str], phase: float = 0):
         super().__init__(name, timeline)
+        assert len(src_list) == 2
+        self.src_list = src_list
+        self.phase = phase
 
-        self.beamsplitter = FockBeamSplitter(name + ".beamsplitter", timeline)
         for i in range(2):
             d = Detector(name + ".detector" + str(i), timeline)
-            self.beamsplitter.add_receiver(d)
             self.detectors.append(d)
-
-        # default circuit: no phase applied
-        self._circuit = Circuit(2)
+        self.components = self.detectors
 
         self.trigger_times = [[], []]
+        self.detect_info = [[], []]
+        self.arrival_times = [[], []]
+        self.temporary_photon_info = [{}, {}]
+
+        self.povms = [None] * 4
+
+    def init(self):
+        self._generate_povms()
+        super().init()
+
+    def _generate_transformed_ladders(self):
+        """Method to generate transformed creation/annihilation operators by the beamsplitter.
+
+        Will be used to construct POVM operators.
+        """
+
+        truncation = self.timeline.quantum_manager.truncation
+        identity = eye(truncation + 1)
+        create, destroy = self.timeline.quantum_manager.build_ladder()
+        phase = self.phase
+        efficiency1 = sqrt(self.detectors[0].efficiency)
+        efficiency2 = sqrt(self.detectors[1].efficiency)
+
+        # Modified mode operators in Heisenberg picture by beamsplitter transformation
+        # considering inefficiency and ignoring relative phase
+        create1 = (kron(efficiency1*create, identity) + exp(1j*phase)*kron(identity, efficiency2*create)) / sqrt(2)
+        destroy1 = create1.conj().T
+        create2 = (kron(efficiency1*create, identity) - exp(1j*phase)*kron(identity, efficiency2*create)) / sqrt(2)
+        destroy2 = create2.conj().T
+
+        return create1, destroy1, create2, destroy2
 
-        self.components = self.detectors + [self.beamsplitter]
+    def _generate_povms(self):
+        """Method to generate POVM operators corresponding to photon detector having 00, 01, 10 and 11 click(s).
+
+        Will be used to generated outcome probability distribution.
+        """
+
+        # assume using Fock quantum manager
+        truncation = self.timeline.quantum_manager.truncation
+        create1, destroy1, create2, destroy2 = self._generate_transformed_ladders()
+
+        # for detector1 (index 0)
+        series_elem_list1 = [(-1)**i * fractional_matrix_power(create1, i+1).dot(
+            fractional_matrix_power(destroy1, i+1)) / factorial(i+1) for i in range(truncation)]
+        povm1_1 = sum(series_elem_list1)
+        povm0_1 = eye((truncation+1) ** 2) - povm1_1
+        # for detector2 (index 1)
+        series_elem_list2 = [(-1)**i * fractional_matrix_power(create2, i+1).dot(
+            fractional_matrix_power(destroy2,i+1)) / factorial(i+1) for i in range(truncation)]
+        povm1_2 = sum(series_elem_list2)
+        povm0_2 = eye((truncation+1) ** 2) - povm1_2
+
+        # POVM operators for 4 possible outcomes
+        # Note: povm01 and povm10 are relevant to BSM
+        povm00 = povm0_1 @ povm0_2
+        povm01 = povm0_1 @ povm1_2
+        povm10 = povm1_1 @ povm0_2
+        povm11 = povm1_1 @ povm1_2
+
+        self.povms = [povm00, povm01, povm10, povm11]
 
     def get(self, photon, **kwargs):
+        src = kwargs["src"]
+        assert photon.encoding_type["name"] == "fock", "Photon must be in Fock representation."
+        input_port = self.src_list.index(src)  # determine at which input the Photon arrives, an index
+        # record arrival time
+        arrival_time = self.timeline.now()
+        self.arrival_times[input_port].append(arrival_time)
+        # record in temporary photon list
+        assert not self.temporary_photon_info[input_port], \
+            "At most 1 Photon instance should arrive at an input port at a time."
+        self.temporary_photon_info[input_port]["photon"] = photon
+        self.temporary_photon_info[input_port]["time"] = arrival_time
+
+        # judge if there have already been two input Photons arriving simultaneously
+        dict0 = self.temporary_photon_info[0]
+        dict1 = self.temporary_photon_info[1]
+        # if both two dictionaries are non-empty
+        if dict0 and dict1:
+            assert dict0["time"] == dict1["time"], "To realize interference photons must arrive simultaneously."
+            photon0 = dict0["photon"]
+            photon1 = dict1["photon"]
+            key0 = photon0.quantum_state
+            key1 = photon1.quantum_state
+
+            # determine the outcome
+            samp = self.get_generator().random()  # random measurement sample
+            result = self.timeline.quantum_manager.measure([key0, key1], self.povms, samp)
+
+            assert result in list(range(len(self.povms))), "The measurement outcome is not valid."
+            if result == 0:
+                # no click for either detector, but still record the zero outcome
+                # record detection information
+                detection_time = self.timeline.now()
+                info = {"time": detection_time, "outcome": 0}
+                self.detect_info[0].append(info)
+                self.detect_info[1].append(info)
+
+            elif result == 1:
+                # detector 1 has a click
+                # trigger time recording will be done by SPD
+                self.detectors[1].record_detection()
+                # record detection information
+                detection_time = self.timeline.now()
+                info0 = {"time": detection_time, "outcome": 0}
+                info1 = {"time": detection_time, "outcome": 1}
+                self.detect_info[0].append(info0)
+                self.detect_info[1].append(info1)
+
+            elif result == 2:
+                # detector 0 has a click
+                # trigger time recording will be done by SPD
+                self.detectors[0].record_detection()
+                # record detection information
+                detection_time = self.timeline.now()
+                info0 = {"time": detection_time, "outcome": 1}
+                info1 = {"time": detection_time, "outcome": 0}
+                self.detect_info[0].append(info0)
+                self.detect_info[1].append(info1)
+
+            elif result == 3:
+                # both detectors have a click
+                # trigger time recording will be done by SPD
+                self.detectors[0].record_detection()
+                self.detectors[1].record_detection()
+                # record detection information
+                detection_time = self.timeline.now()
+                info = {"time": detection_time, "outcome": 1}
+                self.detect_info[0].append(info)
+                self.detect_info[1].append(info)
+
+            self.temporary_photon_info = [{}, {}]
+
+        else:
+            pass
+
+        """
         # check if we have non-null photon
         if not photon.is_null:
             state = self.timeline.quantum_manager.get(photon.quantum_state)
@@ -355,10 +586,17 @@ def get(self, photon, **kwargs):
                 self.timeline.quantum_manager.run_circuit(self._circuit, state.keys)
 
             self.beamsplitter.get(photon)
+        """
 
     def get_photon_times(self) -> List[List[int]]:
+        """Method to get detector trigger times and detection information.
+        Will clear `trigger_times` and `detect_info`.
+        """
         trigger_times = self.trigger_times
+        # detect_info = self.detect_info
         self.trigger_times = [[], []]
+        self.detect_info = [[], []]
+        # return trigger_times, detect_info
         return trigger_times
 
     # does nothing for this class
@@ -366,5 +604,5 @@ def set_basis_list(self, basis_list: List[int], start_time: int, frequency: floa
         pass
 
     def set_phase(self, phase: float):
-        self._circuit = Circuit(2)
-        self._circuit.phase(0, phase)
+        self.phase = phase
+        self._generate_povms()
diff --git a/src/components/interferometer.py b/src/components/interferometer.py
index 562df745..74cc4ad6 100644
--- a/src/components/interferometer.py
+++ b/src/components/interferometer.py
@@ -11,7 +11,6 @@
 if TYPE_CHECKING:
     from ..kernel.timeline import Timeline
     from ..components.photon import Photon
-    from ..components.detector import Detector
 
 from ..kernel.process import Process
 from ..kernel.entity import Entity
diff --git a/src/components/light_source.py b/src/components/light_source.py
index d7232661..cd7f62bd 100644
--- a/src/components/light_source.py
+++ b/src/components/light_source.py
@@ -6,13 +6,13 @@
 
 from typing import List
 
-from numpy import multiply
+from numpy import multiply, sqrt, zeros, kron, outer
 
 from .photon import Photon
 from ..kernel.entity import Entity
 from ..kernel.event import Event
 from ..kernel.process import Process
-from ..utils.encoding import polarization
+from ..utils.encoding import polarization, fock
 from ..utils import log
 
 
@@ -110,22 +110,50 @@ class SPDCSource(LightSource):
         timeline (Timeline): timeline for simulation
         frequency (float): frequency (in Hz) of photon creation.
         wavelengths (List[float]): wavelengths (in nm) of emitted entangled photons.
+            If a list is given, it should contain two elements (corresponding to two modes).
         linewidth (float): st. dev. in photon wavelength (in nm) (currently unused).
         mean_photon_num (float): mean number of photons emitted each period.
         encoding_type (Dict): encoding scheme of emitted photons (as defined in the encoding module).
         phase_error (float): phase error applied to qubits.
     """
 
-    def __init__(self, name, timeline, wavelengths=None, frequency=8e7, bandwidth=0, mean_photon_num=0.1,
-                 encoding_type=polarization, phase_error=0):
+    def __init__(self, name, timeline, wavelengths=None, frequency=8e7, mean_photon_num=0.1,
+                 encoding_type=fock, phase_error=0, bandwidth=0):
         super().__init__(name, timeline, frequency, 0, bandwidth, mean_photon_num, encoding_type, phase_error)
-        if wavelengths is None:
-            wavelengths = [1550, 1550]
         self.wavelengths = wavelengths
+        if self.wavelengths is None or len(self.wavelengths) != 2:
+            self.set_wavelength()
 
     def init(self):
         assert len(self._receivers) == 2, "SPDC source must connect to 2 receivers."
 
+    def _generate_tmsv_state(self):
+        """Method to generate two-mode squeezed vacuum state of two output photonic modes
+
+        Returns:
+            array: generated state.
+        """
+
+        mean_num = self.mean_photon_num
+        truncation = self.timeline.quantum_manager.truncation
+
+        # create state component amplitudes list
+        amp_list = [(sqrt(mean_num / (mean_num + 1)) ** m) / sqrt(mean_num + 1) for m in range(truncation)]
+        amp_square_list = [amp ** 2 for amp in amp_list]
+        amp_list.append(sqrt(1 - sum(amp_square_list)))
+
+        # create two-mode state vector
+        state_vec = zeros((truncation+1) ** 2)
+
+        for i in range(truncation+1):
+            amp = amp_list[i]
+            basis = zeros(truncation+1)
+            basis[i] = 1
+            basis = kron(basis,basis)
+            state_vec += amp * basis
+
+        return state_vec
+
     def emit(self, state_list):
         """Method to emit photons.
 
@@ -134,14 +162,40 @@ def emit(self, state_list):
 
         Arguments:
             state_list (List[List[complex]]): list of complex coefficient arrays to send as photon-encoded qubits.
-                This is ignored for the absorptive encoding type.
+                This is ignored for absorptive and Fock encoding types.
+                For these encoding types only the length of list matters and elements can be arbitrary.
         """
 
         log.logger.info("SPDC sourcee {} emitting {} photons".format(self.name, len(state_list)))
 
         time = self.timeline.now()
 
-        if self.encoding_type["name"] == "absorptive":
+        if self.encoding_type["name"] == "fock":
+            # Use Fock encoding.
+            # The two generated photons should be entangled and should have keys pointing to same Fock state.
+            for _ in state_list:
+                # generate two new photons
+                new_photon0 = Photon("", self.timeline,
+                                     wavelength=self.wavelengths[0],
+                                     location=self,
+                                     encoding_type=self.encoding_type,
+                                     use_qm=True)
+                new_photon1 = Photon("", self.timeline,
+                                     wavelength=self.wavelengths[1],
+                                     location=self,
+                                     encoding_type=self.encoding_type,
+                                     use_qm=True)
+
+                # set shared state to squeezed state
+                state = self._generate_tmsv_state()
+                keys = [new_photon0.quantum_state, new_photon1.quantum_state]
+                self.timeline.quantum_manager.set(keys, state)
+
+                self.send_photons(time, [new_photon0, new_photon1])
+                self.photon_counter += 1
+                time += 1e12 / self.frequency
+
+        elif self.encoding_type["name"] == "absorptive":
             for _ in state_list:
                 num_photon_pairs = self.get_generator().poisson(self.mean_photon_num)
 
@@ -218,3 +272,7 @@ def send_photons(self, time, photons: List["Photon"]):
             process = Process(dst, "get", [photon])
             event = Event(int(round(time)), process)
             self.timeline.schedule(event)
+
+    def set_wavelength(self, wavelength1=1550, wavelength2=1550):
+        """Method to set the wavelengths of photons emitted in two output modes."""
+        self.wavelengths = [wavelength1, wavelength2]
\ No newline at end of file
diff --git a/src/components/memory.py b/src/components/memory.py
index 8e4f0065..8345d129 100644
--- a/src/components/memory.py
+++ b/src/components/memory.py
@@ -7,7 +7,7 @@
 
 from copy import copy
 from math import inf
-from typing import Any, List, TYPE_CHECKING, Dict, Callable
+from typing import Any, List, TYPE_CHECKING, Dict, Callable, Union
 
 from scipy import stats
 
@@ -22,6 +22,11 @@
 from ..utils.encoding import single_atom
 
 
+def const(t):
+    """Constant function thal always returns 1. For AFC memory default spin efficiency."""
+    return 1
+
+
 # array of single atom memories
 class MemoryArray(Entity):
     """Aggregator for Memory objects.
@@ -290,25 +295,31 @@ class AbsorptiveMemory(Entity):
     """Atomic ensemble absorptive memory.
 
     This class models an AFC(-spinwave) absorptive memory, where the quantum state is stored as collective excitation of atomic ensemble.
-    Retrieved photon sequence might be reversed (only for AFC spinwave), which is physically determine by RF pulses used during spinwave.
-    This class will not support qubit state manipulation.
-    Before invoking methods like "get" and "retrieve", must call "prepare" first to prepare the AFC, will take finite simulation time.
-    This class models an AFC absorptive memory, where the quantum state is stored as collective excitation of an atomic ensemble.
+    Rephasing time (predetermined storage time for AFC type) is given by temporal mode number and length of each temporal mode bin.
     This class does not support qubit state manipulation, individual photons should be manipulated instead.
+
     Before invoking methods like "get" and "retrieve", need to call "prepare" first to prepare the AFC, will take finite simulation time.
-    Rephasing time (predetermined storage time for AFC type) is given by temporal mode number and length of each temporal mode bin.
+    Retrieved photon sequence might be reversed (only for AFC-spinwave), which is physically determined by RF pulses used during spinwave.
     Note that for AFC (not spinwave) type, is_reversed must be False.
 
+    Note that the memory is reusable (available for other rounds of absorption) as long as AFC structure is still usable.
+    Once AFC structure expires the memory also expires, and the memory needs to be re-prepared so stored photons cannot be retrieved.
+    Spinwave coherence time determines whether a certain cycle of absorption and re-emission will be successful.
+    Eg. If on-demand storage time is longer than spinwave coherence time, the information stored in this cycle will be cleared.
+        However, if the AFC structure retains usability (in this example AFC_lifetime > coherence_time), the memory itself is not expired.
+
     Attributes:
         name (str): label for memory instance.
         timeline (Timeline): timeline for simulation.
         fidelity (float): (current) fidelity of memory's entanglement.
         frequency (float): maximum frequency of absorption for memory (total frequency bandwidth of AFC memory) (in Hz).
         absorption_efficiency (float): probability of absorbing a photon when arriving at the memory.
-        efficiency (Callable): probability of emitting a photon as a function of storage time.
+        afc_efficiency (Callable): probability of emitting a photon as a function of AFC re-emission time of optical AFC.
+        spin_efficiency (Callable): effeciency of spinwave storage as a function of storage time.
         mode_number (int): number of temporal modes available for storing photons, i.e. number of peaks in Atomic Frequency Comb.
         mode_bin (int):
-        coherence_time (float): average usable lifetime of memory (in s).
+        afc_lifetime (float): average usable lifetime of AFC structure (in s), 0 means infinite lifetime.
+        coherence_time (float): average usable lifetime of spinwave storage (spinwave transition coherence time) (in s), 0 means infinite coherence time.
         wavelength (float): wavelength (in nm) of absorbed and emitted photons.
         total_time (float): AFC re-phasing time (in ps)
         overlap_error (float): error due to photon overlap in one temporal mode, will degrade fidelity.
@@ -322,12 +333,13 @@ class AbsorptiveMemory(Entity):
         memory_array (MemoryArray): memory array aggregating current memory.
         destination (str): name of predetermined re-emission destination node, default None.
         entangled_memory (Dict[str, Any]): tracks entanglement state of memory with a memory.
-        stored_photons (Dict[str, Any]): photons stored in memory temporal modes.
+        stored_photons (List[Dict]): photons stored in memory temporal modes.
     """
 
-    def __init__(self, name: str, timeline: "Timeline", fidelity: float, frequency: float, absorption_efficiency: float,
-                 efficiency: Callable, mode_number: int, coherence_time: int, wavelength: int, overlap_error: float,
-                 prepare_time: int, is_spinwave=False, is_reversed=False, destination=None):
+    def __init__(self, name: str, timeline: "Timeline", frequency: float, absorption_efficiency: float,
+                 afc_efficiency: Callable, mode_number: int, wavelength: int, prepare_time: int = 0,
+                 afc_lifetime: float = -1, coherence_time: float = -1, fidelity: float = 1, overlap_error: float = 0,
+                 is_spinwave=False, is_reversed=False, destination=None, spin_efficiency=const):
         """Constructor for the AbsorptiveMemory class.
 
         Args:
@@ -336,9 +348,11 @@ def __init__(self, name: str, timeline: "Timeline", fidelity: float, frequency:
             fidelity (float): fidelity of memory.
             frequency (float): maximum frequency of absorption for memory (total frequency bandwidth of AFC memory).
             absorption_efficiency (float): probability of absorbing a photon when arriving at the memory.
-            efficiency (Callable): probability of emitting a photon as a function of storage time.
+            afc_efficiency (Callable): probability of emitting a photon as a function of AFC re-emission time of optical AFC.
+            spin_efficiency (Callable): effeciency of spinwave storage as a function of storage time. Default contant unity.
             mode_number (int): number of modes supported for storing photons.
-            coherence_time (float): average time (in s) that memory state is valid.
+            afc_lifetime (float): average usable lifetime of AFC structure (in s), -1 means infinite lifetime.
+            coherence_time (float): average usable lifetime of spinwave storage (spinwave transition coherence time) (in s), -1 means infinite coherence time.
             wavelength (int): wavelength (in nm) of photons emitted by memories.
             overlap_error (float): error due to photon overlap in one temporal mode.
             prepare_time (float): time to prepare AFC (in ps).
@@ -355,12 +369,14 @@ def __init__(self, name: str, timeline: "Timeline", fidelity: float, frequency:
         self.raw_fidelity = fidelity
         self.frequency = frequency
         self.absorption_efficiency = absorption_efficiency
-        self.efficiency = efficiency
+        self.afc_efficiency = afc_efficiency
+        self.spin_efficiency = spin_efficiency
         self.mode_number = mode_number
-        self.coherence_time = coherence_time  # coherence time in seconds
+        self.afc_lifetime = afc_lifetime  # AFC lifetime in seconds
+        self.coherence_time = coherence_time  # spinwave coherencetime in seconds
         self.wavelength = wavelength
-        self.mode_bin = 1e12 / self.frequency  # time bin for each separate temporal mode
-        self.total_time = self.mode_number * self.mode_bin  # AFC rephasing time
+        self.mode_bin = 1e12 / self.frequency  # time bin for each separate temporal mode in ps
+        self.total_time = self.mode_number * self.mode_bin  # AFC rephasing time in ps
         self.overlap_error = overlap_error
         self.prepare_time = prepare_time
 
@@ -388,12 +404,12 @@ def __init__(self, name: str, timeline: "Timeline", fidelity: float, frequency:
         self.excited_photons = []
 
         # initialization of stored_photons dictionary
-        self.stored_photons = [None] * self.mode_number
+        self.stored_photons: List[Union[None, Dict]] = [None] * self.mode_number
 
     def init(self):
         """Implementation of Entity interface (see base class)."""
 
-        self.photon_counter = 0
+        pass
 
     def set_memory_array(self, memory_array: MemoryArray):
         """Method to set the memory array to which the memory belongs
@@ -425,7 +441,10 @@ def _prepare_AFC(self):
             raise Exception("AFC has already been prepared")
         else:
             self.is_prepared = True
-            self.photon_counter = 0
+
+            # schedule AFC expiration once it is prepared, if finite AFC lifetime
+            if self.afc_lifetime > 0:
+                self._schedule_expiration()
 
     def get(self, photon: "Photon", **kwargs):
         """Method to receive a photon to store in the absorptive memory."""
@@ -434,11 +453,44 @@ def get(self, photon: "Photon", **kwargs):
         if not self.is_prepared:
             raise Exception("AFC is not prepared yet.")
 
-        now = -1
+        self.photon_counter += 1
+        now = self.timeline.now()
+
+        # determine which temporal mode the photon is stored in
+        absorb_time = now - self.absorb_start_time
+        index = int(absorb_time / self.mode_bin)
+        if index < 0 or index >= self.mode_number:
+            return
+
         # require resonant absorption of photons
-        if photon.wavelength == self.wavelength and self.get_generator().random() < self.absorption_efficiency:
-            self.photon_counter += 1
-            now = self.timeline.now()
+        # if photon uses Fock representation, inefficiency will be reflected with loss channel
+        if photon.encoding_type["name"] == "fock" and photon.wavelength == self.wavelength:
+            # invoke loss channel due to absorption inefficiency
+            key = photon.quantum_state  # if using Fock representation, the `quantum_state` field is the state key.
+            loss = 1 - self.absorption_efficiency  # loss rate due to absorption inefficiency
+            # apply loss channel on photonic state and return a new state
+            self.timeline.quantum_manager.add_loss(key, loss)
+
+            # store photon information, NOTE the difference between information recorded for different encodings
+            # in this case, store more than one photon
+            if self.stored_photons[index] is None:
+                self.stored_photons[index] = {"photons": [photon], "times": [absorb_time]}
+                self.excited_photons.append(photon)
+            else:
+                self.stored_photons[index]["photons"].append(photon)
+                self.stored_photons[index]["times"].append(absorb_time)
+
+        # otherwise, use random counter w/ efficiency
+        elif photon.wavelength == self.wavelength and self.get_generator().random() < self.absorption_efficiency:
+            # keep one photon per mode since most hardware cannot resolve photon number
+            # photon_counter might be larger than mode_number, multi-photon events counted by "number"
+            # if "overlap" is True, memory fidelity will be corrected by overlap_error
+            if self.stored_photons[index] is None:
+                self.stored_photons[index] = {"photon": photon, "time": absorb_time, "number": 1, "overlap": False}
+                self.excited_photons.append(photon)
+            else:
+                self.stored_photons[index]["number"] += 1
+                self.stored_photons[index]["overlap"] = True
 
         # determine absorb_start_time
         if self.photon_counter == 1:
@@ -449,36 +501,21 @@ def get(self, photon: "Photon", **kwargs):
                 process = Process(self, "retrieve", [])
                 event = Event(self.absorb_start_time + self.total_time, process)
                 self.timeline.schedule(event)
-
-        # schedule expiration at absorb_start_time
-        if self.coherence_time > 0 and self.absorb_start_time == now:
-            self._schedule_expiration()
-
-        # determine which temporal mode the photon is stored in
-        absorb_time = now - self.absorb_start_time
-        index = int(absorb_time / self.mode_bin)
-        if index < 0 or index >= self.mode_number:
-            return
-
-        # keep one photon per mode since most hardware cannot resolve photon number
-        # photon_counter might be larger than mode_number, multi-photon events counted by "number"
-        # if "degradation" is True, memory fidelity will be corrected by overlap_error
-        if self.stored_photons[index] is None:
-            self.stored_photons[index] = {"photon": photon, "time": absorb_time, "number": 1, "degradation": False}
-            self.excited_photons.append(photon)
-        else:
-            self.stored_photons[index]["number"] += 1
-            self.stored_photons[index]["degradation"] = True
+            else:
+                # if spinwave type, and if finite spin coherence time, schedule spinwave decoherence induced storage resetting
+                if self.coherence_time>0:
+                    self._schedule_storage_reset()
 
     def retrieve(self, dst=""):
-        """Method to re-emit all stored photons in a reverse sequence on demand.
+        """Method to re-emit all stored photons in normal/reverse sequence on demand.
 
         Efficiency is a function of time.
         """
 
         # AFC needs to be prepared first
+        # for simplicity, do not allow memory to re-emit if AFC expires
         if not self.is_prepared:
-            raise Exception("AFC is not prepared yet.")
+            raise Exception("AFC is not prepared yet or has expired.")
 
         # do nothing if there are no photons stored
         if len(self.excited_photons) == 0:
@@ -487,15 +524,43 @@ def retrieve(self, dst=""):
         now = self.timeline.now()
         store_time = now - self.absorb_start_time - self.total_time
 
-        # TODO: clear locally stored photons
         for index in range(self.mode_number):
+
             if self.stored_photons[index] is not None:
-                if self.get_generator().random() < self.efficiency(store_time):
-                    photon = self.stored_photons[index]["photon"]
-                    absorb_time = self.stored_photons[index]["time"]
+                stored_photons = self.stored_photons[index]
+
+                # check if using Fock representation (will include
+                # if photon uses Fock representation, no need for random number judgement
+                # inefficiency will be reflected with loss channel
+                if "photons" in stored_photons:
+                    for photon, absorb_time in zip(stored_photons["photons"], stored_photons["times"]):
+
+                        key = photon.quantum_state
+                        # loss rate due to emission inefficiency
+                        loss = 1 - self.afc_efficiency(self.total_time) * self.spin_efficiency(store_time)
+                        # apply loss channel on photonic state and return a new state
+                        self.timeline.quantum_manager.add_loss(key, loss)
+
+                        if self.is_reversed:
+                            if not self.is_spinwave:
+                                raise Exception("AFC memory can only have normal order of re-emission")
+                            emit_time = self.total_time - self.mode_bin - absorb_time  # reversed order of re-emission
+                        else:
+                            emit_time = absorb_time  # normal order of re-emission
+
+                        if self.destination is not None:
+                            dst = self.destination
+
+                        process = Process(self._receivers[0], "get", [photon], {"dst": dst})
+                        event = Event(self.timeline.now() + emit_time, process)
+                        self.timeline.schedule(event)
+
+                elif self.get_generator().random() < self.afc_efficiency(store_time):
+                    photon = stored_photons["photon"]
+                    absorb_time = stored_photons["time"]
 
                     if self.is_reversed:
-                        if self.is_spinwave:
+                        if not self.is_spinwave:
                             raise Exception("AFC memory can only have normal order of re-emission")
                         emit_time = self.total_time - self.mode_bin - absorb_time  # reversed order of re-emission
                     else:
@@ -504,27 +569,28 @@ def retrieve(self, dst=""):
                     if self.destination is not None:
                         dst = self.destination
 
-                    # process = Process(self.owner, "send_qubit", [dst, photon])
-                    # self._receivers[0].get(photon, dst)
                     process = Process(self._receivers[0], "get", [photon], {"dst": dst})
                     event = Event(self.timeline.now() + emit_time, process)
                     self.timeline.schedule(event)
 
-        # TODO: make sure if after each round of re-emission AFC needs to be prepared before next storage task
-        self.is_prepared = False
+        # clear entanglement and storage information after re-emission
+        # retrieval will re-emit all stored photons and information should no longer be stored
+        # clearance unified as storage_reset method
+        self.storage_reset()
 
     def expire(self) -> None:
-        """Method to handle memory expiration.
+        """Method to handle memory expiration due to AFC expiration.
 
         Side Effects:
             Will notify upper entities of expiration via the `pop` interface.
-            Will modify the quantum state of the memory.
+            Will modify information stored by the memory.
         """
 
         # AFC needs to be prepared first
         if not self.is_prepared:
             raise Exception("AFC is not prepared yet.")
 
+        # if stored photons have not been re-emitted when memory expires
         if self.excited_photons:
             for i in range(len(self.excited_photons)):
                 self.excited_photons[i].is_null = True
@@ -537,33 +603,56 @@ def reset(self) -> None:
         """Method to clear quantum memory.
 
         Will reset memory state to no photon stored and will clear entanglement information.
+        Will reset AFC to unprepared state.
 
         Side Effects:
             Will modify internal parameters and photon storage information.
         """
 
-        self.fidelity = 0
-        self.entangled_memory = {'node_id': None, 'memo_id': None}
-        self.photon_counter = 0
-        self.absorb_start_time = 0
-        self.excited_photons = []
+        self.storage_reset()
         self.is_prepared = False
 
         if self.expiration_event is not None:
             self.timeline.remove_event(self.expiration_event)
             self.expiration_event = None
 
+    def storage_reset(self) -> None:
+        """Method to clear information stored in one cycle of photon storage."""
+
+        self.fidelity = 0
+        self.entangled_memory = {'node_id': None, 'memo_id': None}
+        self.photon_counter = 0
+        self.absorb_start_time = 0
+        # no photon should be retrieved after storage reset
+        # re-emitted photons' is_null attribute should not be modified when memory expires
+        self.excited_photons = []
+        self.stored_photons = [None] * self.mode_number
+
     def _schedule_expiration(self) -> None:
+        """Schedule the expiration of memory according to AFC lifetime."""
+
         if self.expiration_event is not None:
             self.timeline.remove_event(self.expiration_event)
 
-        decay_time = self.timeline.now() + int(self.coherence_time * 1e12)
+        decay_time = self.timeline.now() + int(self.afc_lifetime * 1e12)
         process = Process(self, "expire", [])
         event = Event(decay_time, process)
         self.timeline.schedule(event)
 
         self.expiration_event = event
 
+    def _schedule_storage_reset(self) -> None:
+        """Schedule the resetting of a cycle of photon storage in spinwave.
+
+        Only invoked when the memory type is AFC-spinwave.
+        Storage resetting time is determined by AFC rephasing time and spinwave coherence time.
+        """
+
+        reset_time = self.timeline.now() + int(self.total_time + self.coherence_time * 1e12)
+        process = Process(self, "storage_reset", [])
+        event = Event(reset_time, process)
+        self.timeline.schedule(event)
+
     def update_expire_time(self, time: int):
         """Method to change time of expiration.
 
diff --git a/src/components/optical_channel.py b/src/components/optical_channel.py
index 2dc1c611..13dadaff 100644
--- a/src/components/optical_channel.py
+++ b/src/components/optical_channel.py
@@ -153,8 +153,20 @@ def transmit(self, qubit: "Photon", source: "Node") -> None:
                 time = int(time_bin * (1e12 / self.frequency))
             assert time == self.timeline.now(), "qc {} transmit method called at invalid time".format(self.name)
 
-        # check if photon kept
-        if (self.sender.get_generator().random() > self.loss) or qubit.is_null:
+        # check if photon state using Fock representation
+        if qubit.encoding_type["name"] == "fock":
+            key = qubit.quantum_state  # if using Fock representation, the `quantum_state` field is the state key.
+            # apply loss channel on photonic statex
+            self.timeline.quantum_manager.add_loss(key, self.loss)
+
+            # schedule receiving node to receive photon at future time determined by light speed
+            future_time = self.timeline.now() + self.delay
+            process = Process(self.receiver, "receive_qubit", [source.name, qubit])
+            event = Event(future_time, process)
+            self.timeline.schedule(event)
+
+        # if not using Fock representation, check if photon kept
+        elif (self.sender.get_generator().random() > self.loss) or qubit.is_null:
             if self._receiver_on_other_tl():
                 self.timeline.quantum_manager.move_manage_to_server(
                     qubit.quantum_state)
@@ -173,7 +185,7 @@ def transmit(self, qubit: "Photon", source: "Node") -> None:
             event = Event(future_time, process)
             self.timeline.schedule(event)
 
-        # if photon lost, exit
+        # if not using Fock representation, if photon lost, exit
         else:
             pass
 
diff --git a/src/kernel/quantum_manager.py b/src/kernel/quantum_manager.py
index c4d441d6..4cf4e8d9 100644
--- a/src/kernel/quantum_manager.py
+++ b/src/kernel/quantum_manager.py
@@ -17,13 +17,16 @@
 
 from qutip.qip.circuit import QubitCircuit, Gate
 from qutip.qip.operations import gate_sequence_product
-from numpy import log2
+from numpy import log, array, cumsum, base_repr, zeros
+from scipy.sparse import csr_matrix
+from scipy.special import binom
 
 from .quantum_state import KetState, DensityState
 from .quantum_utils import *
 
 KET_STATE_FORMALISM = "ket_vector"
 DENSITY_MATRIX_FORMALISM = "density_matrix"
+FOCK_DENSITY_MATRIX_FORMALISM = "fock_density"
 
 
 class QuantumManager:
@@ -33,13 +36,17 @@ class QuantumManager:
 
     Attributes:
         states (Dict[int, State]): mapping of state keys to quantum state objects.
-        formalism (str): formalism used for local quantum state representation.
+        truncation (int): maximally allowed number of excited states for elementary subsystems.
+                Default is 1 for qubit.
+        dim (int): subsystem Hilbert space dimension. dim = truncation + 1
     """
 
-    def __init__(self, formalism: str):
+    def __init__(self, formalism: str, truncation: int = 1):
         self.states: Dict[int, State] = {}
         self._least_available: int = 0
         self.formalism: str = formalism
+        self.truncation = truncation
+        self.dim = self.truncation + 1
 
     @abstractmethod
     def new(self, state: any) -> int:
@@ -131,10 +138,22 @@ def set(self, keys: List[int], amplitudes: any) -> None:
             amplitudes (any): Amplitudes to set state to, type determined by type of subclass.
         """
 
-        num_qubits = log2(len(amplitudes))
-        assert num_qubits.is_integer(), "Length of amplitudes should be 2 ** n, where n is the number of keys"
-        num_qubits = int(num_qubits)
-        assert num_qubits == len(keys), "Length of amplitudes should be 2 ** n, where n is the number of keys"
+        # num_subsystems = log(len(amplitudes)) / log(self.dim)
+        # assert self.dim ** int(round(num_subsystems)) == len(amplitudes),\
+        #     "Length of amplitudes should be d ** n, " \
+        #     "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+        #     "Actual amplitude length: {}, dim: {}, num subsystems: {}".format(
+        #         len(amplitudes), self.dim, num_subsystems
+        #     )
+        # num_subsystems = int(round(num_subsystems))
+        # assert num_subsystems == len(keys),\
+        #     "Length of amplitudes should be d ** n, " \
+        #     "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+        #     "Amplitude length: {}, expected subsystems: {}, num keys: {}".format(
+        #         len(amplitudes), num_subsystems, len(keys)
+        #     )
+
+        pass
 
     def remove(self, key: int) -> None:
         """Method to remove state stored at key."""
@@ -328,6 +347,7 @@ def _measure(self, state: List[List[complex]], keys: List[int],
             state (List[complex]): state to measure.
             keys (List[int]): list of keys to measure.
             all_keys (List[int]): list of all keys corresponding to state.
+            meas_samp (float): random sample used for measurement result.
 
         Returns:
             Dict[int, int]: mapping of measured keys to measurement results.
@@ -385,3 +405,301 @@ def _measure(self, state: List[List[complex]], keys: List[int],
             self.states[key] = new_state_obj
 
         return dict(zip(keys, result_digits))
+
+
+class QuantumManagerDensityFock(QuantumManager):
+    """Class to track and manage Fock states with the density matrix formalism."""
+
+    def __init__(self, truncation: int = 1):
+        # default truncation is 1 for 2-d Fock space.
+        super().__init__(DENSITY_MATRIX_FORMALISM, truncation=truncation)
+
+    def new(self, state='gnd') -> int:
+        """Method to create a new state with key
+
+        Args:
+            state (Union[str, List[complex], List[List[complex]]]): amplitudes of new state.
+                Default value is 'gnd': create zero-excitation state with current truncation.
+                Other inputs are passed to the constructor of `DensityState`.
+        """
+
+        key = self._least_available
+        self._least_available += 1
+        if state == 'gnd':
+            gnd = [1] + [0]*self.truncation
+            self.states[key] = DensityState(gnd, [key], truncation=self.truncation)
+        else:
+            self.states[key] = DensityState(state, [key], truncation=self.truncation)
+
+        return key
+
+    def run_circuit(self, circuit: Circuit, keys: List[int], meas_samp=None) -> Dict[int, int]:
+        """Currently the Fock states do not support quantum circuits.
+        This method is only to implement abstract method of parent class and SHOULD NOT be called after instantiation.
+        """
+        raise Exception("run_circuit method of class QuantumManagerDensityFock called")
+
+    def _generate_swap_operator(self, num_systems: int, i: int, j: int):
+        """Helper function to generate swapping unitary.
+
+        Args:
+            num_systems (int): number of subsystems in state
+            i (int): index of first subsystem to swap
+            j (int): index of second subsystem to swap
+
+        Returns:
+            Array[int]: unitary swapping operator
+        """
+
+        size = self.dim ** num_systems
+        swap_unitary = zeros((size, size))
+
+        for old_index in range(size):
+            old_str = base_repr(old_index, self.dim)
+            old_str = old_str.zfill(num_systems)
+            new_str = ''.join((old_str[:i], old_str[j], old_str[i+1:j], old_str[i], old_str[j+1:]))
+            new_index = int(new_str, base=self.dim)
+            swap_unitary[new_index, old_index] = 1
+
+        return swap_unitary
+
+    def _prepare_state(self, keys: List[int]):
+        """Function to prepare states at given keys for operator application.
+
+        Will take composite quantum state and swap subsystems to correspond with listed keys.
+        Should not be called directly, but from method to apply operator or measure state.
+
+        Args:
+            keys (List[int]): keys for states to apply operator to.
+
+        Returns:
+            Tuple(List[List[complex]], List[int]): Tuple containing:
+                1. new state to apply operator to, with keys swapped to be consecutive.
+                2. list of keys corresponding to new state.
+        """
+
+        old_states = []
+        all_keys = []
+
+        # go through keys and get all unique qstate objects
+        for key in keys:
+            qstate = self.states[key]
+            if qstate.keys[0] not in all_keys:
+                old_states.append(qstate.state)
+                all_keys += qstate.keys
+
+        # construct compound state
+        new_state = [1]
+        for state in old_states:
+            new_state = kron(new_state, state)
+
+        # apply any necessary swaps to order keys
+        if len(keys) > 1:
+
+            # generate desired key order
+            start_idx = all_keys.index(keys[0])
+            if start_idx + len(keys) > len(all_keys):
+                start_idx = len(all_keys) - len(keys)
+
+            for i, key in enumerate(keys):
+                i = i + start_idx
+                j = all_keys.index(key)
+                if j != i:
+                    swap_unitary = self._generate_swap_operator(len(all_keys), i, j)
+                    new_state = swap_unitary @ new_state @ swap_unitary.T
+                    all_keys[i], all_keys[j] = all_keys[j], all_keys[i]
+
+        return new_state, all_keys
+
+    def _prepare_operator(self, all_keys: List[int], keys: List[int], operator) -> array:
+        # pad operator with identity
+        left_dim = self.dim ** all_keys.index(keys[0])
+        right_dim = self.dim ** (len(all_keys) - all_keys.index(keys[-1]) - 1)
+        prepared_operator = operator
+
+        if left_dim > 0:
+            prepared_operator = kron(identity(left_dim), prepared_operator)
+        if right_dim > 0:
+            prepared_operator = kron(prepared_operator, identity(right_dim))
+
+        return prepared_operator
+
+    def apply_operator(self, operator: array, keys: List[int]):
+        prepared_state, all_keys = self._prepare_state(keys)
+        prepared_operator = self._prepare_operator(all_keys, keys, operator)
+        new_state = prepared_operator @ prepared_state @ prepared_operator.conj().T
+        self.set(all_keys, new_state)
+
+    def set(self, keys: List[int], state: List[List[complex]]) -> None:
+        """Method to set the quantum state at the given keys.
+
+        The `state` argument should be passed as List[List[complex]], where each internal list is a row.
+        However, the `state` may also be given as a one-dimensional pure state.
+        If the list is one-dimensional, will be converted to matrix with the outer product operation.
+
+        Args:
+            keys (List[int]): list of quantum manager keys to modify.
+            state: quantum state to set input keys to.
+        """
+
+        super().set(keys, state)
+        new_state = DensityState(state, keys, truncation=self.truncation)
+        for key in keys:
+            self.states[key] = new_state
+
+    def set_to_zero(self, key: int):
+        """set the state to ground (zero) state."""
+        gnd = [1] + [0]*self.truncation
+        self.set([key], gnd)
+
+    def build_ladder(self):
+        """Generate matrix of creation and annihilation (ladder) operators on truncated Hilbert space."""
+        truncation = self.truncation
+        data = array([sqrt(i+1) for i in range(truncation)])  # elements in create/annihilation operator matrix
+        row = array([i+1 for i in range(truncation)])
+        col = array([i for i in range(truncation)])
+        create = csr_matrix((data, (row, col)), shape=(truncation+1, truncation+1)).toarray()
+        destroy = create.conj().T
+
+        return create, destroy
+
+    def measure(self, keys: List[int], povms: List[array], meas_samp: float) -> int:
+        """Method to measure subsystems at given keys in POVM formalism.
+
+        Serves as wrapper for private `_measure` method, performing quantum manager specific operations.
+
+        Args:
+            keys (List[int]): list of keys to measure.
+            povms: (List[array]): list of POVM operators to use for measurement.
+            meas_samp (float): random measurement sample to use for computing resultant state.
+
+        Returns:
+            int: measurement as index of matching POVM in supplied tuple.
+        """
+
+        new_state, all_keys = self._prepare_state(keys)
+        return self._measure(new_state, keys, all_keys, povms, meas_samp)
+
+    def _measure(self, state: List[List[complex]], keys: List[int],
+                 all_keys: List[int], povms: List[array], meas_samp: float) -> int:
+        """Method to measure subsystems at given keys in POVM formalism.
+
+        Modifies quantum state of all qubits given by all_keys, post-measurement operator determined
+        by measurement operators which are chosen as square root of POVM operators.
+
+        Args:
+            state (List[List[complex]]): state to measure.
+            keys (List[int]): list of keys to measure.
+            all_keys (List[int]): list of all keys corresponding to state.
+            povms: (List[array]): list of POVM operators to use for measurement.
+            meas_samp (float): random measurement sample to use for computing resultant state.
+
+        Returns:
+            int: measurement as index of matching POVM in supplied tuple.
+        """
+
+        state_tuple = tuple(map(tuple, state))
+        povm_tuple = tuple([tuple(map(tuple, povm)) for povm in povms])
+        new_state = None
+        result = 0
+
+        # calculate meas probabilities and projected states
+        if len(keys) == 1:
+            if len(all_keys) == 1:
+                states, probs = measure_state_with_cache_fock_density(state_tuple, povm_tuple)
+
+            else:
+                key = keys[0]
+                num_states = len(all_keys)
+                state_index = all_keys.index(key)
+                states, probs = \
+                    measure_entangled_state_with_cache_fock_density(state_tuple, state_index, num_states, povm_tuple,
+                                                                    self.truncation)
+
+        else:
+            indices = tuple([all_keys.index(key) for key in keys])
+            states, probs = \
+                measure_multiple_with_cache_fock_density(state_tuple, indices, len(all_keys), povm_tuple,
+                                                         self.truncation)
+
+        # calculate result based on measurement sample.
+        prob_sum = cumsum(probs)
+        for i, (output_state, p) in enumerate(zip(states, prob_sum)):
+            if meas_samp < p:
+                new_state = output_state
+                result = i
+                break
+
+        """
+        # for potential future work
+        result_digits = [int(x) for x in base_repr(result, base=self.dim)[2:]]
+        while len(result_digits) < len(keys):
+            result_digits.insert(0, 0)
+
+        # assign measured states
+        for key, result in zip(keys, result_digits):
+            state = [0] * self.dim
+            state[result] = 1
+            self.set([key], state)
+        """
+
+        for key in keys:
+            self.states[key] = None  # clear the stored state at key (particle destructively measured)
+
+        # assign remaining state
+        if len(keys) < len(all_keys):
+            indices = tuple([all_keys.index(key) for key in keys])
+            new_state_tuple = tuple(map(tuple, new_state))
+            remaining_state = density_partial_trace(new_state_tuple, indices, len(all_keys), self.truncation)
+            remaining_keys = [key for key in all_keys if key not in keys]
+            self.set(remaining_keys, remaining_state)
+
+        return result
+
+    def _build_loss_kraus_operators(self, loss_rate: float, all_keys: List[int], key: int) -> List[array]:
+        """Method to build Kraus operators of a generalized amplitude damping channel.
+
+        This represents the effect of photon loss.
+
+        Args:
+            loss_rate (float): loss rate for the quantum channel.
+            all_keys (List[int]): list of all keys in affected state.
+            key (int): key for subsystem experiencing loss.
+
+        Returns:
+            List[array]: list of generated Kraus operators.
+        """
+
+        assert 0 <= loss_rate <= 1
+        kraus_ops = []
+
+        for k in range(self.dim):
+            total_kraus_op = zeros((self.dim ** len(all_keys), self.dim ** len(all_keys)))
+
+            for n in range(k, self.dim):
+                coeff = sqrt(binom(n, k)) * sqrt(((1-loss_rate) ** (n-k)) * (loss_rate ** k))
+                single_op = zeros((self.dim, self.dim))
+                single_op[n-k, n] = 1
+                total_op = self._prepare_operator(all_keys, [key], single_op)
+                total_kraus_op += coeff * total_op
+
+            kraus_ops.append(total_kraus_op)
+
+        return kraus_ops
+
+    def add_loss(self, key, loss_rate):
+        """Method to apply generalized amplitude damping channel on a *single* subspace corresponding to `key`.
+
+        Args:
+            key (int): key for the subspace experiencing loss.
+            loss_rate (float): loss rate for the quantum channel.
+        """
+
+        prepared_state, all_keys = self._prepare_state([key])
+        kraus_ops = self._build_loss_kraus_operators(loss_rate, all_keys, key)
+        output_state = zeros(prepared_state.shape, dtype=complex)
+
+        for kraus_op in kraus_ops:
+            output_state += kraus_op @ prepared_state @ kraus_op.conj().T
+
+        self.set(all_keys, output_state)
diff --git a/src/kernel/quantum_state.py b/src/kernel/quantum_state.py
index 4ea77be4..a237e49a 100644
--- a/src/kernel/quantum_state.py
+++ b/src/kernel/quantum_state.py
@@ -11,7 +11,7 @@
 from abc import ABC
 from typing import Tuple, Dict, List
 
-from numpy import pi, cos, sin, arange, log2
+from numpy import pi, cos, sin, arange, log, log2
 from numpy.random import Generator
 
 from .quantum_utils import *
@@ -38,7 +38,9 @@ class State(ABC):
         keys (List[int]): list of keys pointing to the state, for use with a quantum manager.
     """
 
-    def __init__(self):
+    def __init__(self, **kwargs):
+        # potential key word arguments for derived classes, e.g. truncation = d-1 for qudit
+
         super().__init__()
 
         self.state = None
@@ -76,19 +78,43 @@ class KetState(State):
     """Class to represent an individual quantum state as a ket vector.
 
     Attributes:
-        state (np.array): state vector. Should be of length 2 ** len(keys).
-        keys (List[int]): list of keys (qubits) associated with this state.
+        state (np.array): state vector. Should be of length d ** len(keys), where d is dimension of elementary
+            Hilbert space. Default is 2 for qubits.
+        keys (List[int]): list of keys (subsystems) associated with this state.
+        truncation (int): maximally allowed number of excited states for elementary subsystems.
+                Default is 1 for qubit. dim = truncation + 1
     """
 
-    def __init__(self, amplitudes: List[complex], keys: List[int]):
+    def __init__(self, amplitudes: List[complex], keys: List[int], truncation: int = 1):
+        """Constructor for ket state class.
+
+        Args:
+            amplitudes
+            truncation (int): maximally allowed number of excited states for elementary subsystems.
+                Default is 1 for qubit. dim = truncation + 1
+        """
         super().__init__()
+        self.truncation = truncation
+        dim = self.truncation + 1  # dimension of element Hilbert space
 
         # check formatting
         assert all([abs(a) <= 1.01 for a in amplitudes]), "Illegal value with abs > 1 in ket vector"
         assert abs(sum([abs(a) ** 2 for a in amplitudes]) - 1) < 1e-5, "Squared amplitudes do not sum to 1"
-        num_qubits = log2(len(amplitudes))
-        assert num_qubits.is_integer(), "Length of amplitudes should be 2 ** n, where n is the number of qubits"
-        assert num_qubits == len(keys), "Length of amplitudes should be 2 ** n, where n is the number of qubits"
+
+        num_subsystems = log(len(amplitudes)) / log(dim)
+        assert dim ** int(round(num_subsystems)) == len(amplitudes),\
+            "Length of amplitudes should be d ** n, " \
+            "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+            "Actual amplitude length: {}, dim: {}, num subsystems: {}".format(
+                len(amplitudes), dim, num_subsystems
+            )
+        num_subsystems = int(round(num_subsystems))
+        assert num_subsystems == len(keys),\
+            "Length of amplitudes should be d ** n, " \
+            "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+            "Amplitude length: {}, expected subsystems: {}, num keys: {}".format(
+                len(amplitudes), num_subsystems, len(keys)
+            )
 
         self.state = array(amplitudes, dtype=complex)
         self.keys = keys
@@ -98,32 +124,51 @@ class DensityState(State):
     """Class to represent an individual quantum state as a density matrix.
 
     Attributes:
-        state (np.array): density matrix values. NxN matrix with N = 2 ** len(keys).
-        keys (List[int]): list of keys (qubits) associated with this state.
+        state (np.array): density matrix values. NxN matrix with N = d ** len(keys), where d is dimension of elementary
+            Hilbert space. Default is d = 2 for qubits.
+        keys (List[int]): list of keys (subsystems) associated with this state.
+        truncation (int): maximally allowed number of excited states for elementary subsystems.
+            Default is 1 for qubit. dim = truncation + 1
     """
 
-    def __init__(self, state: List[List[complex]], keys: List[int]):
+    def __init__(self, state: List[List[complex]], keys: List[int], truncation: int = 1):
         """Constructor for density state class.
 
         Args:
             state (List[List[complex]]): density matrix elements given as a list.
                 If the list is one-dimensional, will be converted to matrix with outer product operation.
             keys (List[int]): list of keys to this state in quantum manager.
+            truncation (int): maximally allowed number of excited states for elementary subsystems.
+                Default is 1 for qubit. dim = truncation + 1
         """
 
         super().__init__()
+        self.truncation = truncation
+        dim = self.truncation + 1  # dimension of element Hilbert space
 
         state = array(state, dtype=complex)
         if state.ndim == 1:
-            state = outer(state.conj(), state)
+            state = outer(state, state.conj())
 
         # check formatting
-        assert abs(trace(array(state)) - 1) < 0.1, "density matrix trace must be 1"
+        assert abs(trace(array(state)) - 1) < 0.01, "density matrix trace must be 1"
         for row in state:
             assert len(state) == len(row), "density matrix must be square"
-        num_qubits = log2(len(state))
-        assert num_qubits.is_integer(), "Dimensions of density matrix should be 2 ** n, where n is the number of qubits"
-        assert num_qubits == len(keys), "Dimensions of density matrix should be 2 ** n, where n is the number of qubits"
+
+        num_subsystems = log(len(state)) / log(dim)
+        assert dim ** int(round(num_subsystems)) == len(state), \
+            "Length of amplitudes should be d ** n, " \
+            "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+            "Actual amplitude length: {}, dim: {}, num subsystems: {}".format(
+                len(state), dim, num_subsystems
+            )
+        num_subsystems = int(round(num_subsystems))
+        assert num_subsystems == len(keys), \
+            "Length of amplitudes should be d ** n, " \
+            "where d is subsystem Hilbert space dimension and n is the number of subsystems. " \
+            "Amplitude length: {}, expected subsystems: {}, num keys: {}".format(
+                len(state), num_subsystems, len(keys)
+            )
 
         self.state = state
         self.keys = keys
@@ -196,10 +241,19 @@ def set_state(self, state: Tuple[complex]):
         # check formatting of state
         assert all([abs(a) <= 1.01 for a in state]), "Illegal value with abs > 1 in quantum state"
         assert abs(sum([abs(a) ** 2 for a in state]) - 1) < 1e-5, "Squared amplitudes do not sum to 1"
+
         num_qubits = log2(len(state))
-        assert num_qubits.is_integer(), "Length of amplitudes should be 2 ** n, where n is the number of qubits"
+        assert 2 ** int(round(num_qubits)) == len(state), \
+            "Length of amplitudes should be 2 ** n, where n is the number of qubits. " \
+            "Actual amplitude length: {}, num qubits: {}".format(
+                len(state), num_qubits
+            )
+        num_qubits = int(round(num_qubits))
         assert num_qubits == len(self.entangled_states), \
-            "Length of amplitudes should be 2 ** n, where n is the number of qubits"
+            "Length of amplitudes should be 2 ** n, where n is the number of qubits. " \
+            "Num qubits in state: {}, num qubits in object: {}".format(
+                num_qubits, len(self.entangled_states)
+            )
 
         for qs in self.entangled_states:
             qs.state = state
@@ -227,7 +281,9 @@ def measure(self, basis: Tuple[Tuple[complex]], rng: Generator) -> int:
         """Method to measure a single quantum state.
 
         Args:
-            basis (Tuple[Tuple[complex]]): measurement basis, given as list of states (that are themselves lists of complex coefficients).
+            basis (Tuple[Tuple[complex]]): measurement basis, given as list of states
+                (that are themselves lists of complex coefficients).
+            rng (Generator): random number generator for measurement
 
         Returns:
             int: 0/1 measurement result, corresponding to one basis vector.
@@ -276,6 +332,7 @@ def measure_multiple(basis, states, rng: Generator):
         Args:
             basis (List[List[complex]]): list of basis vectors.
             states (List[QuantumState]): list of quantum state objects to measure.
+            rng (Generator): random number generator for measurement
 
         Returns:
             int: measurement result in given basis.
diff --git a/src/kernel/quantum_utils.py b/src/kernel/quantum_utils.py
index 197b0a9b..d6601608 100644
--- a/src/kernel/quantum_utils.py
+++ b/src/kernel/quantum_utils.py
@@ -5,10 +5,11 @@
 """
 
 from functools import lru_cache
-from typing import Tuple
+from typing import List, Tuple
 from math import sqrt
 
 from numpy import array, kron, identity, zeros, trace, outer, eye
+from scipy.linalg import sqrtm
 
 
 a = array([[0, 1], [0, 0]])
@@ -20,6 +21,7 @@
 
 @lru_cache(maxsize=1000)
 def measure_state_with_cache(state: Tuple[complex, complex], basis: Tuple[Tuple[complex]]) -> float:
+
     state = array(state)
     u = array(basis[0], dtype=complex)
     # measurement operator
@@ -31,9 +33,10 @@ def measure_state_with_cache(state: Tuple[complex, complex], basis: Tuple[Tuple[
 
 
 @lru_cache(maxsize=1000)
-def measure_entangled_state_with_cache(state: Tuple[complex], basis: Tuple[Tuple[complex]],
-                                        state_index: int, num_states: int) -> Tuple[
-    Tuple[complex], Tuple[complex], float]:
+def measure_entangled_state_with_cache(state: Tuple[complex], basis: Tuple[Tuple[complex]], state_index: int,
+                                       num_states: int) -> \
+        Tuple[array, array, float]:
+
     state = array(state)
     u = array(basis[0], dtype=complex)
     v = array(basis[1], dtype=complex)
@@ -70,7 +73,8 @@ def measure_entangled_state_with_cache(state: Tuple[complex], basis: Tuple[Tuple
 
 @lru_cache(maxsize=1000)
 def measure_multiple_with_cache(state: Tuple[complex], basis: Tuple[Tuple[complex]], length_diff: int) \
-        -> Tuple[Tuple[Tuple[complex]], Tuple[float]]:
+        -> Tuple[List[array], List[float]]:
+
     state = array(state)
     # construct measurement operators, projectors, and probabilities of measurement
     projectors = [None] * len(basis)
@@ -88,14 +92,14 @@ def measure_multiple_with_cache(state: Tuple[complex], basis: Tuple[Tuple[comple
         # project to new state
         if probabilities[i] > 0:
             new_state = (proj @ state) / sqrt(probabilities[i])
-            new_state = tuple(new_state)
             return_states[i] = new_state
 
-    return tuple(return_states), tuple(probabilities)
+    return return_states, probabilities
 
 
 @lru_cache(maxsize=1000)
 def measure_state_with_cache_ket(state: Tuple[complex, complex]) -> float:
+
     state = array(state)
     M0 = array([[1, 0], [0, 0]], dtype=complex)
 
@@ -105,8 +109,9 @@ def measure_state_with_cache_ket(state: Tuple[complex, complex]) -> float:
 
 
 @lru_cache(maxsize=1000)
-def measure_entangled_state_with_cache_ket(state: Tuple[complex], state_index: int, num_states: int) -> Tuple[
-        Tuple[complex], Tuple[complex], float]:
+def measure_entangled_state_with_cache_ket(state: Tuple[complex], state_index: int, num_states: int) \
+        -> Tuple[array, array, float]:
+
     state = array(state)
 
     # generate projectors
@@ -133,12 +138,13 @@ def measure_entangled_state_with_cache_ket(state: Tuple[complex], state_index: i
     else:
         state0 = (projector0 @ state) / sqrt(prob_0)
 
-    return (state0, state1, prob_0)
+    return state0, state1, prob_0
 
 
 @lru_cache(maxsize=1000)
-def measure_multiple_with_cache_ket(state: Tuple[complex], num_states: int, length_diff: int) -> Tuple[
-        Tuple[Tuple[complex]], Tuple[float]]:
+def measure_multiple_with_cache_ket(state: Tuple[complex], num_states: int, length_diff: int) \
+        -> Tuple[List[array], List[float]]:
+
     state = array(state)
     basis_count = 2 ** num_states
 
@@ -163,11 +169,12 @@ def measure_multiple_with_cache_ket(state: Tuple[complex], num_states: int, leng
             new_state = tuple(new_state)
             return_states[i] = new_state
 
-    return (tuple(return_states), tuple(probabilities))
+    return return_states, probabilities
 
 
 @lru_cache(maxsize=1000)
 def measure_state_with_cache_density(state: Tuple[Tuple[complex, complex]]) -> float:
+
     state = array(state)
     M0 = array([[1, 0], [0, 0]], dtype=complex)
 
@@ -177,8 +184,9 @@ def measure_state_with_cache_density(state: Tuple[Tuple[complex, complex]]) -> f
 
 
 @lru_cache(maxsize=1000)
-def measure_entangled_state_with_cache_density(state: Tuple[Tuple[complex]], state_index: int, num_states: int) -> Tuple[
-        Tuple[complex], Tuple[complex], float]:
+def measure_entangled_state_with_cache_density(state: Tuple[Tuple[complex]], state_index: int, num_states: int) \
+        -> Tuple[array, array, float]:
+
     state = array(state)
 
     # generate projectors
@@ -205,12 +213,13 @@ def measure_entangled_state_with_cache_density(state: Tuple[Tuple[complex]], sta
     else:
         state0 = (projector0 @ state @ projector0) / prob_0
 
-    return (state0, state1, prob_0)
+    return state0, state1, prob_0
 
 
 @lru_cache(maxsize=1000)
-def measure_multiple_with_cache_density(state: Tuple[Tuple[complex]], num_states: int, length_diff: int) -> Tuple[
-        Tuple[Tuple[complex]], Tuple[float]]:
+def measure_multiple_with_cache_density(state: Tuple[Tuple[complex]], num_states: int, length_diff: int) \
+        -> Tuple[List[array], List[float]]:
+
     state = array(state)
     basis_count = 2 ** num_states
 
@@ -235,5 +244,173 @@ def measure_multiple_with_cache_density(state: Tuple[Tuple[complex]], num_states
             new_state = tuple(new_state)
             return_states[i] = new_state
 
-    return (tuple(return_states), tuple(probabilities))
+    return return_states, probabilities
+
+
+@lru_cache(maxsize=1000)
+def measure_state_with_cache_fock_density(state: Tuple[Tuple[complex]], povms: Tuple[Tuple[Tuple[complex]]]) \
+        -> Tuple[List[array], List[float]]:
+    state = array(state)
+    povms = [array(povm) for povm in povms]
+
+    # probabilities of getting different outcomes according to POVM operators
+    prob_list = [trace(state @ povm).real for povm in povms]
+    state_list = []
+
+    # get output states
+    for i in range(len(prob_list)):
+        if prob_list[i] <= 0:
+            state_post_meas = None
+        else:
+            measure_op = sqrtm(povms[i])
+            state_post_meas = (measure_op @ state @ measure_op) / prob_list[i]
+
+        state_list.append(state_post_meas)
+
+    return state_list, prob_list
+
+
+@lru_cache(maxsize=1000)
+def measure_entangled_state_with_cache_fock_density(state: Tuple[Tuple[complex]], system_index: int, num_systems: int,
+                                                    povms: Tuple[Tuple[Tuple[complex]]], truncation: int = 1) \
+        -> Tuple[List[array], List[float]]:
+
+    """Measure one subsystem of a larger composite system.
+
+    The measurement SHOULD NOT be entangling measurement, and thus POVM operators should be precisely consisted of
+    operators on the subsystem's Hilbert space alone.
+
+    Args:
+        state (Tuple[Tuple[complex]]): state to measure
+        system_index (int): index of measured subsystem within state.
+        num_systems (int): number of total systems in the state.
+        povms (Tuple[Tuple[Tuple[complex]]]): tuple listing all POVM operators to use for measurement
+        truncation (int): fock space truncation, 1 for qubit system (default 1).
+
+    Returns:
+        Tuple[List[array], List[float]]: tuple with two sub-lists.
+            The first lists each output state, corresponding with the measurement of each POVM.
+            The second lists the probability for each measurement.
+    """
+
+    state = array(state)
+    povms = [array(povm) for povm in povms]
+
+    # generate POVM operators on total Hilbert space
+    povm_list = []
+    left_dim = (truncation + 1) ** system_index
+    right_dim = (truncation + 1) ** (num_systems - system_index - 1)
+    for povm in povms:
+        povm_tot = kron(kron(identity(left_dim), povm), identity(right_dim))
+        povm_list.append(povm_tot)
+
+    # list of probabilities of getting different outcomes from POVM
+    prob_list = [trace(state @ povm).real for povm in povm_list]
+    state_list = []
+
+    for i in range(len(prob_list)):
+        if prob_list[i] <= 0:
+            state_post_meas = None
+        else:
+            measure_op = sqrtm(povm_list[i])
+            state_post_meas = (measure_op @ state @ measure_op) / prob_list[i]
+
+        state_list.append(state_post_meas)
+
+    return state_list, prob_list
+
+
+@lru_cache(maxsize=1000)
+def measure_multiple_with_cache_fock_density(state: Tuple[Tuple[complex]], indices: Tuple[int], num_systems: int,
+                                             povms: Tuple[Tuple[Tuple[complex]]], truncation: int = 1) \
+        -> Tuple[List[array], List[float]]:
+
+    """Measure multiple subsystems of a larger composite system.
+
+    Should be called by Quantum Managers.
+    This function will facilitate entangling measurement, e.g. BSM with two photon detectors behind a beamsplitter.
+    Such measurement operators are consisted of mixed operators on different subsystems' Hilbert spaces.
+    For current implementation, the Hilbert space on which those mixed operators act on is separated from the rest of
+    the total Hilbert space, and we assume that the involved subsystems have already been moved close in terms of keys.
+    i.e., elements in `indices` argument are no less than 0 and no greater than `num_systems`
+    (relative indices w.r.t. the measured state), and the elements MUST BE consecutive.
+
+    E.g. For a total system consisted of 4 subsystems (0, 1, 2, 3) each with dimension d,
+    if the entangling measurement happens on (1, 2), then measurement operators on total space will be constructed as
+        O_tot = kron(kron(identity(d), O), identity(d)),
+    where the measurement operator on (1, 2) subspace needs to be generated beforehand to feed in the function.
+
+    Args:
+        state (Tuple[Tuple[complex]]): state to measure.
+        indices (Tuple[int]): indices within combined state to measure.
+        num_systems (int): number of total systems in the state.
+        povms (Tuple[Tuple[Tuple[complex]]]): tuple listing all POVM operators to use for measurement.
+        truncation (int): fock space truncation, 1 for qubit system (default 1).
+
+    Returns:
+        Tuple[List[array], List[float]]: tuple with two sub-lists.
+            The first lists each output state, corresponding with the measurement of each POVM.
+            The second lists the probability for each measurement.
+    """
+
+    state = array(state)
+    povms = [array(povm) for povm in povms]
+
+    # judge if elements in `indices` are consecutive
+    init_meas_sys_idx = min(indices)
+    fin_meas_sys_idx = max(indices)
+    num = len(indices)
+    if (fin_meas_sys_idx - init_meas_sys_idx + 1 != num) or (list(indices) != sorted(indices)):
+        raise ValueError("Indices should be consecutive; got {}".format(indices))
+
+    povm_list = []
+    left_dim = (truncation + 1) ** init_meas_sys_idx
+    right_dim = (truncation + 1) ** (num_systems - fin_meas_sys_idx - 1)
+    for povm in povms:
+        povm_tot = kron(kron(identity(left_dim), povm), identity(right_dim))
+        povm_list.append(povm_tot)
+
+    # list of probabilities of getting different outcomes from POVM
+    prob_list = [trace(state @ povm).real for povm in povm_list]
+    state_list = []
+
+    for i in range(len(prob_list)):
+        if prob_list[i] <= 0:
+            state_post_meas = None
+        else:
+            measure_op = sqrtm(povm_list[i])
+            state_post_meas = (measure_op @ state @ measure_op) / prob_list[i]
+
+        state_list.append(state_post_meas)
+
+    # return post-measurement states and measurement outcome probabilities in the order of fed-in POVM operators
+    return state_list, prob_list
+
+
+@lru_cache(maxsize=1000)
+def density_partial_trace(state: Tuple[Tuple[complex]], indices: Tuple[int], num_systems: int, truncation: int = 1) \
+        -> array:
+
+    """Traces out subsystems systems at given indices.
+
+    Args:
+        state (Tuple[Tuple[complex]]: input state.
+        indices (Tuple[int]): indices of subsystems to trace out of state.
+            should be sorted in increasing order.
+        num_systems (int): number of total subsystems in the state.
+        truncation (int): fock space truncation, 1 for qubit system (default 1).
+
+    Returns:
+        array: output state with reduced number of subsystems `num_systems - len(indices)`.
+    """
+
+    temp = array(state)
+
+    for i, idx in enumerate(indices):
+        offset = num_systems - i
+        temp = temp.reshape((truncation+1,) * offset * 2)
+        temp = trace(temp, axis1=(idx-i), axis2=(offset+idx-i))
 
+    output_dim = (truncation + 1) ** (num_systems - len(indices))
+    output_state = temp.reshape((output_dim, output_dim))
+    return output_state
diff --git a/src/kernel/timeline.py b/src/kernel/timeline.py
index 989c0344..3cea4fbf 100644
--- a/src/kernel/timeline.py
+++ b/src/kernel/timeline.py
@@ -21,8 +21,10 @@
 from ..utils import log
 from .quantum_manager import (QuantumManagerKet,
                               QuantumManagerDensity,
+                              QuantumManagerDensityFock,
                               KET_STATE_FORMALISM,
-                              DENSITY_MATRIX_FORMALISM)
+                              DENSITY_MATRIX_FORMALISM,
+                              FOCK_DENSITY_MATRIX_FORMALISM)
 
 CARRIAGE_RETURN = '\r'
 SLEEP_SECONDS = 3
@@ -58,11 +60,13 @@ class Timeline:
         quantum_manager (QuantumManager): quantum state manager.
     """
 
-    def __init__(self, stop_time=inf, formalism=KET_STATE_FORMALISM):
+    def __init__(self, stop_time=inf, formalism=KET_STATE_FORMALISM, truncation=1):
         """Constructor for timeline.
 
         Args:
             stop_time (int): stop time (in ps) of simulation (default inf).
+            formalism (str): formalism of quantum state representation.
+            truncation (int): truncation of Hilbert space (currently only for Fock representation).
         """
         self.events: EventList = EventList()
         self.entities: Dict[str, "Entity"] = {}
@@ -77,6 +81,8 @@ def __init__(self, stop_time=inf, formalism=KET_STATE_FORMALISM):
             self.quantum_manager = QuantumManagerKet()
         elif formalism == DENSITY_MATRIX_FORMALISM:
             self.quantum_manager = QuantumManagerDensity()
+        elif formalism == FOCK_DENSITY_MATRIX_FORMALISM:
+            self.quantum_manager = QuantumManagerDensityFock(truncation=truncation)
         else:
             raise ValueError(f"Invalid formalism {formalism}")
 
@@ -188,10 +194,10 @@ def print_time(self):
             stdout.flush()
             sleep(SLEEP_SECONDS)
 
-    def ns_to_human_time(self, nanoseconds: int) -> str:
+    def ns_to_human_time(self, nanoseconds: float) -> str:
         milliseconds = nanoseconds / NANOSECONDS_PER_MILLISECOND
         return str(timedelta(milliseconds=milliseconds))
 
     @staticmethod
-    def convert_to_nanoseconds(picoseconds: int) -> int:
+    def convert_to_nanoseconds(picoseconds: int) -> float:
         return picoseconds / PICOSECONDS_PER_NANOSECOND
diff --git a/src/topology/node.py b/src/topology/node.py
index 59faf96a..5aa72142 100644
--- a/src/topology/node.py
+++ b/src/topology/node.py
@@ -170,7 +170,7 @@ def receive_qubit(self, src: str, qubit) -> None:
 
         self.components[self.first_component_name].get(qubit)
 
-    def get_components_by_type(self, component_type: str) -> [Entity]:
+    def get_components_by_type(self, component_type: str) -> List[Entity]:
         return [comp for comp in self.components.values() if type(comp).__name__ == component_type]
 
     def change_timeline(self, timeline: "Timeline"):
diff --git a/src/utils/encoding.py b/src/utils/encoding.py
index cadafd4a..edd33661 100644
--- a/src/utils/encoding.py
+++ b/src/utils/encoding.py
@@ -37,3 +37,8 @@
     {"name": "absorptive",
      "bases": [((complex(1), complex(0)), (complex(0), complex(1))), None]
      }
+
+fock = \
+    {"name": "fock",
+     "bases": None
+     }


Test patch:
diff --git a/example/photon_test.py b/example/photon_test.py
index 22ddef48..9b6774fc 100644
--- a/example/photon_test.py
+++ b/example/photon_test.py
@@ -1,7 +1,7 @@
 import time
 
 from numpy import sin, cos
-from numpy.random import random_sample
+from numpy.random import random_sample, default_rng
 from sequence.components.circuit import Circuit
 from sequence.kernel.quantum_manager import QuantumManagerKet
 from sequence.utils.quantum_state import QuantumState
@@ -10,6 +10,7 @@
 SAMPLE_SIZE = 10000
 
 basis = ((complex(1), complex(0)), (complex(0), complex(1)))
+rng = default_rng()
 qm = QuantumManagerKet()
 meas_circ = Circuit(1)
 meas_circ.measure(0)
@@ -33,7 +34,7 @@ def __init__(self, state):
         self.qs.state = state
 
     def measure(self):
-        return self.qs.measure(basis)
+        return self.qs.measure(basis, rng)
 
 
 class ParticleNew:
@@ -41,7 +42,7 @@ def __init__(self, state):
         self.qs_key = qm.new(state)
 
     def measure(self):
-        return qm.run_circuit(meas_circ, [self.qs_key])
+        return qm.run_circuit(meas_circ, [self.qs_key], rng.random())
 
 
 if __name__ == "__main__":
diff --git a/parallel/tests_parallel/kernel/test_p_timeline.py b/parallel/tests_parallel/kernel/test_p_timeline.py
index ce921b09..f772aa07 100644
--- a/parallel/tests_parallel/kernel/test_p_timeline.py
+++ b/parallel/tests_parallel/kernel/test_p_timeline.py
@@ -6,9 +6,6 @@
 
 from psequence.p_timeline import ParallelTimeline
 
-rank = MPI.COMM_WORLD.Get_rank()
-size = MPI.COMM_WORLD.Get_size()
-
 
 class FakeEntity(Entity):
     def __init__(self, name, tl):
@@ -22,7 +19,7 @@ def add(self):
         self.counter += 1
 
 
-def build_env(lookahead):
+def build_env(lookahead, rank, size):
     tl = ParallelTimeline(lookahead)
 
     entity = FakeEntity(str(rank), tl)
@@ -34,15 +31,22 @@ def build_env(lookahead):
 
 @pytest.mark.mpi
 def test_p_timeline_schedule_local_events():
-    tl, entity = build_env(10)
+    rank = MPI.COMM_WORLD.Get_rank()
+    size = MPI.COMM_WORLD.Get_size()
+
+    tl, entity = build_env(10, rank, size)
     assert len(tl.events) == 0
     event = Event(20, Process(entity, "add", []))
     tl.schedule(event)
     assert len(tl.events) == 1
 
 
+@pytest.mark.mpi
 def test_p_timeline_schedule_remote_events():
-    tl, entity = build_env(10)
+    rank = MPI.COMM_WORLD.Get_rank()
+    size = MPI.COMM_WORLD.Get_size()
+
+    tl, entity = build_env(10, rank, size)
     assert len(tl.events) == 0
     foreign_rank = (int(entity.name) + 1) % size
     foreign_entity_name = str(foreign_rank)
@@ -52,8 +56,12 @@ def test_p_timeline_schedule_remote_events():
     assert len(tl.event_buffer[foreign_rank]) == 1
 
 
+@pytest.mark.mpi
 def test_p_timeline_run():
-    tl, entity = build_env(10)
+    rank = MPI.COMM_WORLD.Get_rank()
+    size = MPI.COMM_WORLD.Get_size()
+
+    tl, entity = build_env(10, rank, size)
     assert len(tl.events) == 0
     foreign_rank = (int(entity.name) + 1) % size
     foreign_entity_name = str(foreign_rank)
diff --git a/tests/components/test_detector.py b/tests/components/test_detector.py
index 5b03ed5f..faf11246 100644
--- a/tests/components/test_detector.py
+++ b/tests/components/test_detector.py
@@ -4,7 +4,8 @@
 from sequence.components.detector import *
 from sequence.components.photon import Photon
 from sequence.kernel.timeline import Timeline
-from sequence.utils.encoding import polarization, time_bin, absorptive
+from sequence.kernel.quantum_manager import FOCK_DENSITY_MATRIX_FORMALISM
+from sequence.utils.encoding import polarization, time_bin, absorptive, fock
 
 SEED = 0
 
@@ -44,7 +45,7 @@ def get_generator(self):
 def test_Detector_init():
     detector, parent, tl = create_detector(dark_count=10)
     tl.init()
-    assert len(tl.events) > 0
+    assert len(tl.events) == 2
 
 
 def test_Detector_get():
@@ -203,7 +204,7 @@ def test_QSDetectorFockDirect():
     period = (1e12 / COUNT_RATE) + 1
     src_list = ["a", "b"]
 
-    tl = Timeline()
+    tl = Timeline(formalism=FOCK_DENSITY_MATRIX_FORMALISM)
 
     qsd = QSDetectorFockDirect("qsd", tl, src_list)
     [qsd.update_detector_params(i, "efficiency", 1) for i in range(2)]
@@ -212,7 +213,7 @@ def test_QSDetectorFockDirect():
     tl.init()
 
     for _ in range(1000):
-        photon = Photon("", tl, encoding_type=absorptive, use_qm=True)
+        photon = Photon("", tl, encoding_type=fock, use_qm=True)
         photon.set_state((0, 1))
         qsd.get(photon, src=src_list[0])
         tl.time += period
@@ -222,7 +223,7 @@ def test_QSDetectorFockDirect():
 
     tl.time = 0
     for _ in range(1000):
-        photon = Photon("", tl, encoding_type=absorptive, use_qm=True)
+        photon = Photon("", tl, encoding_type=fock, use_qm=True)
         photon.set_state((0, 1))
         qsd.get(photon, src=src_list[1])
         tl.time += period
@@ -243,10 +244,11 @@ def get_generator(self):
     COUNT_RATE = 80e6
     period = (1e12 / COUNT_RATE) + 1
     psi_minus = [complex(0), complex(sqrt(1 / 2)), -complex(sqrt(1 / 2)), complex(0)]
+    src_list = ["a", "b"]
 
-    tl = Timeline(formalism="density_matrix")
+    tl = Timeline(formalism=FOCK_DENSITY_MATRIX_FORMALISM)
 
-    qsd = QSDetectorFockInterference("qsd", tl)
+    qsd = QSDetectorFockInterference("qsd", tl, src_list)
     random_control = RandomControl(0)
     qsd.owner = random_control
     [qsd.update_detector_params(i, "efficiency", 1) for i in range(2)]
@@ -254,23 +256,25 @@ def get_generator(self):
 
     tl.init()
 
-    # measure unentangled
-    for _ in range(NUM_TRIALS):
-        p0 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        qsd.get(p0)
-        tl.time += period
-
-    times = qsd.get_photon_times()
-    assert abs(len(times[0]) / len(times[1])) - 1 < 0.1
+    # # measure unentangled
+    # for _ in range(NUM_TRIALS):
+    #     p0 = Photon("", tl, encoding_type=fock, use_qm=True)
+    #     qsd.get(p0, src=src_list[0])
+    #     tl.time += period
+    #
+    # times = qsd.get_photon_times()
+    # assert abs(len(times[0]) / len(times[1])) - 1 < 0.1
 
     # measure unentangled, two incident real photons
     tl.time = 0
     clear_qsd_detectors(qsd)
     for _ in range(NUM_TRIALS):
-        p0 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p1 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        qsd.get(p0)
-        qsd.get(p1)
+        p0 = Photon("", tl, encoding_type=fock, use_qm=True)
+        p1 = Photon("", tl, encoding_type=fock, use_qm=True)
+        p0.set_state((complex(0), complex(1)))
+        p1.set_state((complex(0), complex(1)))
+        qsd.get(p0, src=src_list[0])
+        qsd.get(p1, src=src_list[1])
         tl.time += period
 
     times = qsd.get_photon_times()
@@ -280,47 +284,47 @@ def get_generator(self):
     tl.time = 0
     clear_qsd_detectors(qsd)
     for _ in range(NUM_TRIALS):
-        p0 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p1 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p0.is_null = True
-        p0.combine_state(p1)
-        p0.set_state(psi_minus)
-        qsd.get(p0)
-        qsd.get(p1)
+        p0 = Photon("", tl, encoding_type=fock, use_qm=True)
+        p1 = Photon("", tl, encoding_type=fock, use_qm=True)
+        key0 = p0.quantum_state
+        key1 = p1.quantum_state
+        tl.quantum_manager.set([key0, key1], psi_minus)
+        qsd.get(p0, src=src_list[0])
+        qsd.get(p1, src=src_list[1])
         tl.time += period
 
     times = qsd.get_photon_times()
-    assert abs(len(times[0]) / len(times[1])) - 1 < 0.1
+    assert len(times[1]) == NUM_TRIALS
 
     # measure entangled, pi/2 phase
     tl.time = 0
     clear_qsd_detectors(qsd)
     qsd.set_phase(np.pi/2)
     for _ in range(NUM_TRIALS):
-        p0 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p1 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p0.is_null = True
-        p0.combine_state(p1)
-        p0.set_state(psi_minus)
-        qsd.get(p0)
-        qsd.get(p1)
+        p0 = Photon("", tl, encoding_type=fock, use_qm=True)
+        p1 = Photon("", tl, encoding_type=fock, use_qm=True)
+        key0 = p0.quantum_state
+        key1 = p1.quantum_state
+        tl.quantum_manager.set([key0, key1], psi_minus)
+        qsd.get(p0, src=src_list[0])
+        qsd.get(p1, src=src_list[1])
         tl.time += period
 
     times = qsd.get_photon_times()
-    assert len(times[1]) == NUM_TRIALS
+    assert abs(len(times[0]) / len(times[1]) - 1) < 0.1
 
-    # measure entangled, 3pi/2 phase
+    # measure entangled, pi phase
     tl.time = 0
     clear_qsd_detectors(qsd)
-    qsd.set_phase(3*np.pi/2)
+    qsd.set_phase(np.pi)
     for _ in range(NUM_TRIALS):
-        p0 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p1 = Photon("", tl, encoding_type=absorptive, use_qm=True)
-        p0.is_null = True
-        p0.combine_state(p1)
-        p0.set_state(psi_minus)
-        qsd.get(p0)
-        qsd.get(p1)
+        p0 = Photon("", tl, encoding_type=fock, use_qm=True)
+        p1 = Photon("", tl, encoding_type=fock, use_qm=True)
+        key0 = p0.quantum_state
+        key1 = p1.quantum_state
+        tl.quantum_manager.set([key0, key1], psi_minus)
+        qsd.get(p0, src=src_list[0])
+        qsd.get(p1, src=src_list[1])
         tl.time += period
 
     times = qsd.get_photon_times()
diff --git a/tests/components/test_memory.py b/tests/components/test_memory.py
index 505d692c..ac2a909a 100644
--- a/tests/components/test_memory.py
+++ b/tests/components/test_memory.py
@@ -203,7 +203,8 @@ def test_Memory__schedule_expiration():
 def test_Absorptive_prepare():
     PREPARE_TIME = 100
     tl = Timeline()
-    mem = AbsorptiveMemory("mem", tl, 1, 80e6, 1, perfect_efficiency, 100, 0, 500, 0, PREPARE_TIME)
+    mem = AbsorptiveMemory("mem", tl, 80e6, 1, perfect_efficiency, 100, 1550, PREPARE_TIME)
+    # mem = AbsorptiveMemory("mem", tl, 80e6, 1, perfect_efficiency, 100, 0, 500, 0, PREPARE_TIME)
 
     tl.init()
     mem.prepare()
@@ -222,7 +223,7 @@ def test_Absorptive_get_all_bins():
     PERIOD = 1
     MODE_NUM = 100
     tl = Timeline()
-    mem = AbsorptiveMemory("mem", tl, 1, 1e12/PERIOD, 1, perfect_efficiency, MODE_NUM, 0, 500, 0, 100)
+    mem = AbsorptiveMemory("mem", tl, 1e12/PERIOD, 1, perfect_efficiency, MODE_NUM, 500)
     mem.is_prepared = True
     tl.init()
 
@@ -249,7 +250,7 @@ def test_Absorptive_get_all_bins():
 def test_Absorpive_get_skip_bins():
     PERIOD = 1
     tl = Timeline()
-    mem = AbsorptiveMemory("mem", tl, 1, 1e12/PERIOD, 1, perfect_efficiency, 100, 0, 500, 0, 100)
+    mem = AbsorptiveMemory("mem", tl, 1e12/PERIOD, 1, perfect_efficiency, 100, 500)
     mem.is_prepared = True
     tl.init()
 
@@ -271,7 +272,7 @@ def test_Absorptive_retrieve():
     PERIOD = 1
     MODE_NUM = 100
     tl = Timeline()
-    mem = AbsorptiveMemory("mem", tl, 1, 1e12/PERIOD, 1, perfect_efficiency, MODE_NUM, 0, 500, 0, 100)
+    mem = AbsorptiveMemory("mem", tl, 1e12/PERIOD, 1, perfect_efficiency, MODE_NUM, 500)
     parent = DumbParent(mem)
     mem.is_prepared = True
     tl.init()
@@ -296,9 +297,18 @@ def test_Absorptive_retrieve():
 
     # retrieve photons in reverse order
     mem.is_prepared = True
+    mem.is_spinwave = True
     mem.is_reversed = True
     parent.reset()
     tl.time = 0
+    mem.absorb_start_time = 0
+    photons = [None] * MODE_NUM
+    for i in range(MODE_NUM):
+        photon = Photon(str(i), tl, 500)
+        photons[i] = photon
+        mem.stored_photons[i] = {"photon": photon, "time": i}
+        mem.excited_photons.append(photon)
+
     mem.retrieve()
     assert len(tl.events) == MODE_NUM
     tl.run()
@@ -313,7 +323,7 @@ def test_Absorptive_retrieve():
 
 def test_Absorptive_expire():
     tl = Timeline()
-    mem = AbsorptiveMemory("mem", tl, 1, 80e6, 1, perfect_efficiency, 100, 0, 500, 0, 100)
+    mem = AbsorptiveMemory("mem", tl, 80e6, 1, perfect_efficiency, 100, 500)
     parent = DumbParent(mem)
 
     process = Process(mem, "expire", [])
diff --git a/tests/kernel/test_quantum_manager.py b/tests/kernel/test_quantum_manager.py
index 38dac3f7..75c11c5c 100644
--- a/tests/kernel/test_quantum_manager.py
+++ b/tests/kernel/test_quantum_manager.py
@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.linalg import fractional_matrix_power
 import math
 
 from sequence.kernel.quantum_manager import *
@@ -341,3 +342,214 @@ def test_qmanager__measure_density():
 
     assert abs((len(meas_0) / NUM_TESTS) - 0.5) < 0.1
 
+
+def test_qmanager__prepare_state_fock():
+    qm = QuantumManagerDensityFock(truncation=2)
+    state_0 = [1, 0, 0]
+    state_1 = [0, 1, 0]
+    combined = np.kron(state_0, state_1)
+    desired_state = np.outer(combined, combined.conj())
+
+    # test single state
+    key_0 = qm.new()
+    key_1 = qm.new()
+    qm.set([key_0, key_1], combined)
+
+    new_state, all_keys = qm._prepare_state([key_0, key_1])
+    assert np.all(new_state == np.array(desired_state))
+    assert np.all(np.array(all_keys) == np.array([key_0, key_1]))
+
+    # test disjoint state
+    key_0 = qm.new(state_0)
+    key_1 = qm.new(state_1)
+
+    new_state, all_keys = qm._prepare_state([key_0, key_1])
+    assert np.all(new_state == np.array(desired_state))
+    assert np.all(np.array(all_keys) == np.array([key_0, key_1]))
+
+    # test state that needs swapping
+    key_1 = qm.new()
+    key_0 = qm.new()
+    combined_2 = np.kron(state_1, state_0)
+    qm.set([key_1, key_0], combined_2)
+
+    new_state, all_keys = qm._prepare_state([key_0, key_1])
+    assert np.all(new_state == np.array(desired_state))
+    assert np.all(np.array(all_keys) == np.array([key_0, key_1]))
+
+    # test larger state that needs swapping
+    key_0 = qm.new()
+    key_1 = qm.new()
+    key_2 = qm.new()
+    combined_3 = np.kron(np.kron(state_0, state_0), state_1)
+    qm.set([key_0, key_1, key_2], combined_3)
+
+    combined_correct = np.kron(np.kron(state_0, state_1), state_0)
+    desired_state_alt = np.outer(combined_correct, combined_correct.conj())
+    new_state, all_keys = qm._prepare_state([key_0, key_2])
+    assert np.all(new_state == np.array(desired_state_alt))
+    assert np.all(np.array(all_keys) == np.array([key_0, key_2, key_1]))
+
+
+def test_qmanager_build_ladder_fock():
+    TRUNCATION = 2
+
+    qm = QuantumManagerDensityFock(truncation=TRUNCATION)
+    create, destroy = qm.build_ladder()
+
+    # check dims
+    assert len(create.shape) == 2
+    assert len(destroy.shape) == 2
+    assert create.shape[0] == create.shape[1] == TRUNCATION + 1
+    assert destroy.shape[0] == destroy.shape[1] == TRUNCATION + 1
+
+    # check properties
+    state = [0] * (TRUNCATION + 1)
+    state[-1] = 1
+    state = array(state)
+    state_create = create @ state
+    assert all(state_create == 0)
+
+    state = [0] * (TRUNCATION + 1)
+    state[0] = 1
+    state = array(state)
+    state_destroy = destroy @ state
+    assert all(state_destroy == 0)
+
+    for i in range(TRUNCATION):
+        state = [0] * (TRUNCATION + 1)
+        state[i] = 1
+        state = array(state)
+
+        state_create = create @ state
+        assert state_create[i+1] == sqrt(i+1)
+        assert sum(state_create) == sqrt(i+1)
+
+    for i in range(TRUNCATION):
+        state = [0] * (TRUNCATION + 1)
+        state[i+1] = 1
+        state = array(state)
+
+        state_destroy = destroy @ state
+        assert state_destroy[i] == sqrt(i+1)
+        assert sum(state_destroy) == sqrt(i+1)
+
+
+def test_qmanager_apply_operator_fock():
+    TRUNCATION = 2
+
+    qm = QuantumManagerDensityFock(truncation=2)
+    create, destroy = qm.build_ladder()
+    excited = np.zeros((TRUNCATION + 1, TRUNCATION + 1))
+    excited[1, 1] = 1
+
+    # single state
+    key = qm.new()
+    qm.apply_operator(create, [key])
+    new_state = qm.get(key)
+    assert np.all(new_state.state == excited)
+
+    # compound state
+    key1 = qm.new()
+    key2 = qm.new()
+    oper = np.eye((TRUNCATION + 1) ** 2)
+    qm.apply_operator(oper, [key1, key2])
+    qm.apply_operator(create, [key1])
+    new_state = qm.get(key1)
+
+    desired2 = np.zeros((TRUNCATION + 1, TRUNCATION + 1))
+    desired2[0, 0] = 1
+    desired = np.kron(excited, desired2)
+    assert np.all(new_state.state == desired)
+
+    # comound operator
+    key1 = qm.new()
+    key2 = qm.new()
+    oper = np.kron(create, create)
+    qm.apply_operator(oper, [key1, key2])
+    new_state = qm.get(key1)
+
+    desired = np.kron(excited, excited)
+    assert np.all(new_state.state == desired)
+
+
+def test_qmanager_measure_fock():
+    NUM_TESTS = 1000
+    TRUNCATION = 2
+
+    qm = QuantumManagerDensityFock(truncation=2)
+    create, destroy = qm.build_ladder()
+    series_elem_list = [((-1) ** i) * fractional_matrix_power(create, i + 1).dot(
+        fractional_matrix_power(destroy, i + 1)) / math.factorial(i + 1) for i in range(TRUNCATION)]
+    povm1 = sum(series_elem_list)
+    povm0 = eye(TRUNCATION + 1) - povm1
+
+    # single state
+    meas_0 = []
+    meas_1 = []
+    state_single = [math.sqrt(1 / 2), math.sqrt(1 / 2), 0]
+    state = np.outer(state_single, state_single)
+    for _ in range(NUM_TESTS):
+        key = qm.new(state)
+        samp = np.random.random()
+        res = qm.measure([key], [povm0, povm1], samp)
+        if res:
+            meas_1.append(key)
+        else:
+            meas_0.append(key)
+
+    assert abs((len(meas_0) / NUM_TESTS) - 0.5) < 0.1
+    for key in meas_0:
+        assert qm.states[key] is None
+    for key in meas_1:
+        assert qm.states[key] is None
+
+    # single state in multi-qubit system
+    meas_0 = []
+    meas_1 = []
+    for _ in range(NUM_TESTS):
+        key1 = qm.new(state)
+        key2 = qm.new()
+        samp = np.random.random()
+        # compound
+        oper = np.identity((TRUNCATION + 1) ** 2)
+        qm.apply_operator(oper, [key1, key2])
+        res = qm.measure([key1], [povm0, povm1], samp)
+        if res:
+            meas_1.append(key2)
+        else:
+            meas_0.append(key2)
+
+    assert abs((len(meas_0) / NUM_TESTS) - 0.5) < 0.1
+
+    desired = [0] * (TRUNCATION + 1)
+    desired[0] = 1
+    desired_density = np.outer(desired, desired)
+    for key in meas_0:
+        assert np.all(qm.get(key).state == desired_density)
+    for key in meas_1:
+        assert np.all(qm.get(key).state == desired_density)
+
+    # multiple state
+    povm00 = np.kron(povm0, povm0)
+    povm01 = np.kron(povm0, povm1)
+    povm10 = np.kron(povm1, povm0)
+    povm11 = np.kron(povm1, povm1)
+    meas_0 = []
+    meas_2 = []
+    for _ in range(NUM_TESTS):
+        key1 = qm.new(state)
+        key2 = qm.new()
+        samp = np.random.random()
+        # compound
+        oper = np.identity((TRUNCATION + 1) ** 2)
+        qm.apply_operator(oper, [key1, key2])
+        res = qm.measure([key1, key2], [povm00, povm01, povm10, povm11], samp)
+        if res == 2:
+            meas_2.append(key1)
+        elif res == 0:
+            meas_0.append(key1)
+        else:
+            raise Exception()
+
+    assert abs((len(meas_0) / NUM_TESTS) - 0.5) < 0.1