Parallel Simulation of Quantum Networks with Distributed Quantum State Management

2.1 Quantum Information and Quantum Network Protocols

Superposition allows a quantum system to be simultaneously in two different states. These states can be mathematically represented as a probability distribution of the possible measurement outcomes on the quantum system. For example, the measurement of a particle in the superposition state \(| \psi \rangle := \alpha | 0 \rangle + \beta | 1 \rangle\) yields the measurement outcome \(|0 \rangle\) with probability \(|\alpha |^2\) and \(|1 \rangle\) with probability \(|\beta |^2\), where \(\alpha\) and \(\beta\) are two complex numbers. Upon measurement, the quantum state \(| \psi \rangle\) collapses to either \(|0\rangle\) or \(|1\rangle\). Optical quantum networks use photons to carry quantum states.

Multiparticle quantum systems can exhibit entanglement, where the quantum state of each particle in the group cannot be described independently of the other particles. The simplest example is the EPR pair named after Einstein, Podolsky, and Rosen [13]: \(| \phi ^+ \rangle := \frac{1}{\sqrt {2}} (| 00 \rangle + | 11 \rangle)\). Here the measurement outcomes of the two particles are perfectly correlated with a 50% chance of being \(|0\rangle\), \(|0\rangle\) and 50% chance of being \(|1\rangle\), \(|1\rangle\). The measurement of the first particle immediately collapses the state of the other particle to the same state, even if the particles are spatially separated. We note that with the increasing number of particles, the space required to represent quantum states grows exponentially, presenting a significant challenge for computer simulations.

Fig. 1.

View Figure

One can use entanglement to teleport an arbitrary quantum state over a long distance without worrying about state loss due to channel attenuation [4]. The teleportation protocol is depicted in Figure 1(a). Solid lines represent qubits (a two-state quantum system), double lines classical bits, circles quantum states, and boxes quantum operations (quantum gates). Initially the sender and the receiver share an entangled qubit pair (green circles) that is distributed between them using the quantum network as described below. The arbitrary quantum state \(|\psi \rangle\) we want to transmit is prepared in another qubit (blue circle). The sender performs the indicated elementary operations on its two qubits (CNOT and Hadamard gates [35]) and measures them. One of four possible measurement outcomes can be realized: 00, 01, 10, or 11. The measurement outcome is then sent to the receiver by using a classical channel, and the receiver uses it to perform controlled quantum operations on its qubit in a process called Pauli frame correction [26]. After the correction, the receiver’s qubit (red circle) will be in state \(|\psi \rangle\). Teleportation is one example of a quantum network protocol, and it allows quantum state transmission by distributing only EPR pairs and classical information.

Reliable, long-distance entanglement distribution is an important service provided by quantum networks; and many other protocols depend on it, including teleportation. The success probability of entanglement distribution decreases exponentially with distance. Fortunately, the entanglement swapping protocol [56] can be used to “stitch together” multiple shorter-distance entanglements for a single longer-distance entanglement. The swapping protocol is illustrated in Figure 1(b). In the first step, entanglement is generated between the adjacent nodes A-B and B-C. Entangled photons may be stored in quantum memories to ensure that they are available precisely when needed [23]. To create entanglement between nodes A-C, node B performs a Bell state measurement (BSM) on its two qubits. The measurement outcome is sent from node B to either node A or node C to perform Pauli frame correction. After the correction, the qubits on nodes A and C become entangled in the EPR state.

In this work, we simulate entanglement distribution and entanglement swapping in networks with up to 1,024 nodes. Since entanglement link generation and swapping involve the measurement of qubits, the final state of the remaining two qubits remains bounded with dimension four. The space required for state representation thus grows only linearly with the size of our network. Despite this state simplification, such large-scale simulation still involves unique challenges as discussed in Section 3 and remains a computationally challenging problem.

2.2 SeQUeNCe: A Customizable Quantum Network Simulator

SeQUeNCe is an open source quantum network simulator [54] that allows simulations with photon-level details. The SeQUeNCe project website [41] provides the complete source code, documentation, and tutorials that allow easy usability and extendability.² The many use cases of the simulator include understanding the tradeoffs of alternative quantum network architectures, optimizing quantum hardware, and developing a robust control plane. SeQUeNCe uses a modularized design, as shown in Figure 2. This enables researchers to customize the simulated network and reuse existing models in a flexible way. The Application module represents quantum network applications and their service requests. The resource management and network management modules use classical control messages to manage the allocation of network resources. The entanglement management and hardware modules encompass protocols and hardware operations that act on qubits that are located in quantum memories or in photons at telecommunication wavelengths in optical network links. The SeQUeNCe simulator stores these states in the quantum state manager. For sequential simulation, a local QSM is sufficient to trace the quantum states. However, parallel simulation requires the simulator to synchronize multiple QSMs, a requirement that is addressed in this work. This work also significantly redesigns the simulation kernel to allow parallel simulations. Further details about the design and use cases of the sequential version of the SeQUeNCe simulator can be found in our earlier work [52—54].

Fig. 2.

View Figure

2.3 Parallel Simulation of Classical Networks

Researchers have built time-parallel and space-parallel network simulators to improve scalability in terms of both execution speed and network size. Time-parallel simulation methods [1, 24, 51] parallelize a fixed-size problem by partitioning the simulation time axis into intervals and assigning a processor to simulate the system over its assigned time interval. Parallel (and distributed) discrete-event methods use a space-parallel approach that divides a single network simulation topology into distinct pieces and executes them in parallel on multiple interconnected processors. The efforts to parallelize simulations of classical networks started more than 20 years ago with Parallel/Distributed ns [42] and the Georgia Tech Network Simulator [43]. Distributed simulation support was added to ns-3 [38] in 2011 [40].

Fig. 3.

View Figure

Ensuring that events are processed in time-stamped order is a fundamental part of distributed simulation—the so-called causality constraint. Figure 3 illustrates the virtual time for both sequential and parallel simulation of a classical network. For the sequential case (insert on the left), the simulation uses a global virtual clock, while the parallel case (insert on the right) has three processes, each with its own corresponding local virtual time that has to be synchronized with the global virtual time. Almost all parallel network simulators use a conservative synchronization algorithm that takes precautions to avoid violating the causality constraint. The synchronization overhead significantly limits the speedup that can be achieved with parallel and distributed simulations. In fact, synchronization overheads in conservative approaches can cause parallel execution to be slower than sequential execution [25]. In contrast, optimistic synchronization algorithms allow violations to occur but are able to detect and recover from them. Optimistic synchronization has been applied successfully in classical network contexts such as HPC, storage, and HPC interconnects [6, 7].

3.1 Requirement Analysis

The accuracy of the simulation depends on the correct order of event execution and the accurate modeling of the network state. During a parallel simulation, a large network is partitioned into multiple non-overlapping subnetworks. Events in a subnetwork are executed by one assigned process. Event execution on one process may affect the order of events and the state of the network on the other processes. The accuracy of the simulation requires the simulator to synchronize these effects among processes. We summarize three methods that can affect the state of nodes on the other subnetworks.

The first two methods affect states on other processes by generating events in the future. Processes can synchronize their states by exchanging cross-process events, and the synchronization algorithm can successfully avoid out-of-order event execution. The third method affects states on other processes by changing a shared state. To ensure consistency of shared states, we have entangled qubits share a single entanglement state instance instead of synchronizing copies of the state across processes. A read/write lock is applied to avoid conflict when accessing those shared states.

Random number generation is another factor affecting reproducibility. Given the randomness in measuring qubits, the parallel simulator should generate random numbers identical to those of the sequential simulation. Therefore, we assign a random number generator (RNG) to every network node in both the sequential and parallel simulation. Events executed on a node or its components use the RNG on the node to generate a random number. Since components connected to a single network node are never split between processes, identical results will be produced as long as events on the node are executed in the correct sequential order.

The scalability of the simulator depends on the reduction of parallelization overhead. In a parallel quantum network simulator, the overhead results from synchronizing processes, exchanging data and accessing the shared quantum states. Our design aims to reduce the execution time and complexity of these functions to improve the simulation scalability.

The usability of the parallel simulator depends on its ease of use. Models developed for parallel simulation ought to be compatible with the sequential simulation (one can configure the parallel simulator to run experiments with one processor). Therefore, our design should allow efficient and error-free simulation design; in other words, the parallelization should be systematically handled by the simulation kernel, with the process being completely transparent to users.

3.2 Parallel Simulator Architecture

Figure 4 depicts the architecture of our parallel simulator for quantum networks. We parallelize SeQUeNCe by leveraging multiprocessing, a technique in which separate processes handle separate subnetworks. We use MPI [48] for interprocess communication including cross-process events, produced by cross-process quantum and classical channels, for both shared- and distributed-memory systems. These two techniques allow users to run a parallel simulation on an MPI cluster. Although multiprocessing may incur a higher communication cost than multithreading does, HPC cluster optimizations have reduced the MPI communication overhead [10]. Our experiments in Section 6 also demonstrate that MPI communication has a limited impact on the parallelization overhead.

Fig. 4.

View Figure

To track simulated quantum states, we design a QSM that provides an interface for manipulating and reading quantum states. For shared entanglement, we also design a global QSM as a separate server process. During the simulation, all cross-process entanglement states are stored in the global QSM, while non-entangled qubits or those entangled with only local qubits are stored in a local QSM. Entangled qubits hold the key for accessing their shared state only within a local or global QSM. An individual process may then communicate with the global QSM to retrieve or update the state of an entangled qubit, either directly or through a quantum circuit. Quantum circuit operations, represented as unitary matrices, are not generally commutative. These operations therefore cannot be reordered, implying a requirement for proper ordering with respect to the simulation time for global QSM requests. Since quantum circuits may operate only on qubits stored in local hardware, however, circuits on separate processes will not overlap. Thus, the execution order of requests within one synchronization window does not affect the final results. The global QSM may be implemented as a single-process program, a multithreaded program, or a cluster according to the scale of the parallel simulation. Thus, the global QSM can handle requests from different processes more efficiently and reduce potential resource competition. While interprocess events are exchanged by using MPI, the simpler point-to-point communication between the local QSMs and the global QSM utilizes TCP/IP communication. This is to avoid the additional overhead from using MPI for QSM communication. Further, QSM implementation details are described in Section 4.

Fig. 5.

View Figure

To ensure proper execution of events, we adjust the conservative synchronization algorithm [33], not only to exchange cross-process events but also to synchronize states on the global QSM. The existing conservative algorithm has been demonstrated to work for networks without entangled states. Figure 5 depicts the operation of the conservative algorithm. Events on separate processes execute in parallel during the Computing phase. If an event execution schedules a future event to be run on a different process, then the future event is stored in the output queue OutQ. When a process executes all local events up to \(sync\_time\), it enters the Waiting phase until all processes are finished. This \(sync\_time\) is chosen to guarantee that all events in the OutQ will have future timestamps not yet reached by the simulation engine (as discussed below). The processes then exchange events in the OutQ during the Communicating phase. Additionally, a future \(sync\_time\) is calculated by using a lookahead value (lookahead).

We describe our implementation of the conservative algorithm with the addition of a global QSM in Algorithm 1. In this algorithm, we introduce the following data structures:

We also include the key configuration parameters lookahead (as defined above) and \(end\_time\) (defined below).

At each iteration of the algorithm loop, we check to see whether \(local\_time\) (the simulation time of event execution on the process) is less than \(end\_time\) (the simulation time at which event execution should halt). If the condition is met, then the algorithm moves into the Communicating phase; otherwise, the loop exits because the simulation has finished.

In the Communicating phase, MPI processes first calculate the minimum timestamp of events stored locally (\(local\_min\_time\)) by comparing the minimum timestamp of the events in the queue (Queue) and output queue (OutQ). Next, MPI processes move events in the local output queue to the corresponding process and broadcast their \(local\_min\_time\). The function “MPI_EXCHANGE_DATA” returns events from other processes to the input queue (InQ) and a list of \(local\_min\_time\) from other processes in the \(remote\_min\_times\) list. By calculating the minimum value of the \(remote\_min\_times\) list and \(local\_min\_time\), the algorithm determines \(global\_min\_time\), which is the minimum timestamp of all events in the simulation. If all processes have an empty event queue, then \(global\_min\_time\) will be infinity, and execution is halted. Otherwise, we can use \(global\_min\_time\) and the lookahead value to get the time window [34] bounded by a future \(sync\_time\) (in simulation time). If the time window extends beyond \(end\_time\), then we instead set \(sync\_time\) to be \(end\_time\).

As shown in Figure 5(a), the lookahead value is set so that all cross-process events may be scheduled for simulation times only after \(sync\_time\). We subsequently move into the Computing phase, where events in the local Queue are executed and \(local\_time\) is updated accordingly. After all the events before \(sync\_time\) are executed, we update \(local\_time\) to be \(sync\_time\). Additionally, the process sends a message to the global QSM to serve all pending requests. The process is unblocked upon receiving a response from the global QSM to ensure that all earlier requests on the global QSM are completely handled.

The parallel simulation kernel also supports sequential models for compatibility when needed. As shown in Figure 4, models can schedule events and manipulate quantum states using interfaces offered by the kernel. When the simulation is executed sequentially, the event scheduler does not need to communicate with other processes, and the local QSM does not need to connect with the global QSM. Therefore, the same model can be executed by using sequential or parallel simulation with only a simple setup configuration change.

3.3 Parallelization Overhead Sources

The proposed architecture introduces parallelization overhead. This overhead may be measured as an increase in execution time beyond that required for the execution of simulated events. Figure 5(b) shows the state of processes using the conservative synchronization algorithm. Processes have three states: Computing, Communicating, and Waiting. The Computing state denotes that the process is executing events whose timestamps are within the time window. The Communicating state denotes the MPI communication that creates a synchronization barrier among processes. The Waiting state denotes that a process has completed the execution of local events in the time window and is now waiting for all other processes to reach the communication barrier. We identify four sources of overhead. First, the barrier caused by the MPI communication makes one process (e.g., with a light workload) wait for other processes (e.g., with heavy workloads) in the Waiting state. The second source is due to the data exchange among processes. The number of exchanged events affects the communication overhead. The third source is the cost for barrier synchronization (e.g., the calculation of the local minimum timestamp for the next computing phase and the time window until the next synchronization). Large lookahead is essential for minimizing parallel synchronization overhead; this will be discussed in Section 5. Apart from the overhead for the synchronization algorithm, the global QSM introduces the fourth overhead source. In the parallel simulation, the global QSM handles operations on cross-process entangled states. As a result, the time spent on communication with the global QSM increases the overall execution time for manipulating entangled states. A parallel simulation with a larger amount of cross-process entanglement may thus introduce larger parallelization overhead due to higher global QSM communication. In Section 4 we introduce multiple techniques for reducing the overhead from these sources.

4.1 Local Quantum State Manager

The QSM within the simulation kernel is designed to provide methods for manipulating quantum states. Models of hardware and protocols may manipulate the quantum state of qubits with a qubit key that is produced by a universal unique identifier [28]. The QSM maintains the relation between keys and quantum states as a hash table. For entangled states, multiple keys may point to the same quantum state object. Here we introduce three important functions provided by the QSM:

Quantum states (including entanglement) are stored as amplitude vectors (Bra-ket notation [44]) or matrices (density matrix [16]) together with a list that indicates the order of entangled qubits as a list of keys. The get and set functions allow models to read and write the quantum state of one or many qubits, including any entanglement. Compared with the set function, hardware and protocol models mostly use the run function to change the quantum state of qubits. For example, entanglement purification [11] and swapping [56] protocols usually describe their protocols as a quantum circuit. A quantum circuit consists of quantum logic gates and/or measurements. The gates in the circuit may be represented by a unitary matrix. Unlike quantum logic gates, measurement is a stochastic nonreversible operation. The simulator thus needs to use a random number to determine the measured state. To reproduce identical outcomes, the run function takes the random number \(prob\_sample\) from the caller function if the circuit contains a measurement operation. As described in Section 3.1, this random \(prob\_sample\) is produced by the random number generator on the calling node to ensure the consistency of simulations.

In SeQUeNCe, we use the QuTiP [22] Python library to generate the unitary matrix of a given circuit. The output state is produced by matrix multiplication of the unitary matrix of the circuit gates with the input quantum state of the qubit(s), which is a CPU-intensive task. Furthermore, the simulation of a quantum network involves repeated computation of this matrix multiplication. Therefore, we utilize the memoization technique [32, 46] to speed up redundant functions, such as the run function. We implement least recently used caches with fixed sizes to cache the result of these multiplications.

4.2 Global Quantum State Manager

Although unentangled states or entangled states confined to one process may be managed by a single local QSM, cross-process entangled states may need to be accessed by multiple local QSMs. We thus designed a global QSM as a server process to ensure the consistency and accuracy of quantum states within the parallel simulation. While the quantum state of qubits is stored in the centralized global QSM, the local QSMs on separate MPI processes provide the transparent functions to manipulate the quantum state wherever the storage of state is.

Since the global QSM is executed on an independent process, it must receive requests and send responses over the network, in contrast to a local QSM. Multiple MPI processes may communicate with the global QSM simultaneously, a process that can result in resource competition. Such competition may increase the processing time on the global QSM when performing large-scale parallel simulations. The communication time and resource competition constitute the parallelization overhead from the global QSM, as described in Section 3.3. The global QSM may thus become a bottleneck for parallel simulation. To address this, we implement multiple techniques to efficiently reduce the global QSM overhead.

4.2.1 Mitigating Communication Overhead.

Fig. 6.

View Figure

The total communication time depends on the number of transmissions and the latency of each transmission. In this work we reduced the number of transmissions between the local and global QSMs. These transmissions include three types of communication categorized by their goals:

—	Forward function call. The local QSM forwards function calls to the global QSM and waits for a response.
—	Synchronize quantum state. Operations on the local QSM may cause a state inconsistency, which requires the local QSM to update the state in the global QSM.
—	Acknowledge state of global QSM. This is for a local QSM to confirm that the global QSM has completed executing all the forwarded functions.

To reduce the number of transmissions, we first offload most of the computational work to the local QSM. We use the entanglement generation using the meet-in-the-middle protocol [23] as an example to show interactions between the local and global QSM in Figure 6. Qubits \(q_A\) and \(q_B\) stored in the quantum memories are entangled with photons \(q_a\) and \(q_b\) with the shared quantum state \(S_A\) (\(S_B\)). Photons are sent to the intermediate BSM node to generate entanglement between \(q_A\) and \(q_B\). In the parallel simulation, qubits either belong to one process, shown in Figure 6(a), or belong to two processes, shown in Figure 6(b). If the MPI process contains all required qubits, then the local QSM can produce the entangled state between \(q_A\) and \(q_B\) locally without forwarding the function to the global QSM. Otherwise, the local QSM must forward the function to the global QSM since the “local QSM 2” executing events on the BSM node do not hold the state \(S_A\). As a result, a simulation with fewer cross-process entanglement states has fewer communication requests with the global QSM and thus reduces the communication overhead.

The challenge of offloading work to the local QSMs is in ensuring the consistency of quantum states between the local and global QSMs. For example, in Figure 6(b) the quantum state \(S_A\) in the global QSM should be identical to the state \(S_A\) stored in “local QSM 1” before the global QSM executes the request from “local QSM 2.” The timing of synchronizing local states with the global QSM thus becomes an issue, since the parallel simulator must keep this consistency while avoiding redundant synchronization.

Investigation of existing entanglement protocols shows that the synchronization communication is triggered by two conditions: (1) transmitted qubits sharing an entangled state and (2) forwarded function calls requiring a local quantum state. The first condition may separate two entangled qubits from one process to two processes, so we synchronize the entanglement state to the global QSM. The management of this state is moved from the local QSM to the global QSM. To the best of our knowledge, all entanglement generation protocols rely on quantum channels to entangle qubits in two locations. The synchronization triggered by the first condition can avoid inconsistent states between the local and global QSMs. The second condition requires the global QSM to execute operations. The local QSM should push the local quantum state to the global QSM before forwarding the function call. The correctness of entanglement protocols, such as entanglement swapping, can be protected by these triggering conditions.

The synchronization requirement demands the management of qubits at the global QSM. The local QSM forwards the function calls related to either unknown qubits (i.e., qubits produced by other processes) or qubits managed by the global QSM. Once the qubits lose their entangled state, the management of qubits will be transferred back to the local QSM. In the parallel SeQUeNCe simulator, either the measurement on qubits or the set function can return the qubits back to the local QSM by destroying an entangled state.

In addition to offloading work from the global QSM as much as possible, we reduce communication time by grouping multiple requests into one packet. We observe that some requests to the global QSM do not expect responses, such as the run function without measurement and the set function. Therefore, the simulator caches requests with no responses. These requests are transmitted to the global QSM only if a new request needs results from the global QSM or the simulator forcefully transmits all requests in the cache using the notify_global_QSM_about_synchronization() function from Algorithm 1. This mechanism leads to efficient state synchronization since the synchronization does not require a response from the global QSM.

We implement the communication between the local and global QSM using the TCP/IP protocol. Function calls are serialized in JSON format. We execute the global QSM process on the same cluster where MPI processes execute the events in parallel run, to reduce the communication latency.

6.1 Chain of Quantum Routers

Fig. 8.

View Figure

The first kind of network has a quantum router chain that serves only one flow that distributes entanglement between the two end routers. The quantum network uses a chain of routers to mitigate the exponentially increased loss rate of photons caused by the increased distance [5]. Figure 8(a) illustrates the topology of a quantum network with five routers. In this type of network, adjacent nodes first use the meet-in-the-middle generation scheme [23] to distribute entanglement over a short distance. Then, the entanglement swapping [56] and purification protocols [11] are used to distribute entanglement over a longer distance with high fidelity. For the parallel simulation experiments, we evenly divide the number of routers by the number of processes (\(n = 1024/p\)) and assign them sequentially to each process (i.e., routers 1 to n will be assigned process 1, routers n+1 to 2n to process 2, etc.). For simulations with the increased lookahead value introduced in Section 5, the BSM nodes that connect two different processes are placed in a separate process where the simulation time is allowed to lag.

All our linear topology simulations use a 1 km distance between all adjacent routers. The BSM node is in the middle of each link; in other words, the length of the optical fiber is 0.5 km. All routers contain 100 quantum memories.

We first evaluate the strong scaling of the simulator in a quantum network with 1,024 routers and a total of 102,400 quantum memories. Strong scaling requires simulations with a fixed total workload and an increasing number of processes. We simulate the network with 1 to 128 processes. For this and future simulations using 1 process, the sequential version of SeQUeNCe was used (rather than the parallel version running with 1 MPI process). Since one computing node contains only 36 cores, parallel simulations with 64 and 128 processes are executed on 2 and 4 computing nodes, respectively. Figure 8(b) illustrates the average execution time of each simulation process for 10 separate simulations and the speedup ratio achieved by the parallel simulation. Since the standard deviation of the measured execution time ranges from 7.4 to 189.6 s for strong scaling (with 32 and 1 processes, respectively) and from 0.8 to 13.48 s for weak scaling (with 16 and 128 processes, respectively), we do not include error bars in the plot. We observe that the execution time decreases with the increasing number of processes. The parallel simulation can speed up the simulation by a factor of up to 25. When we double the number of processes, the computation time is halved (indicated by the blue part). The synchronization time (indicated by the red part) is roughly constant, which leads to a growing portion of the total execution time with the growing number of processes (from 0.07% to 0.67% of all execution time) as the computation time drops. We also notice that the time for the socket communication between the local and global QSMs is trivial in this case. The time of socket communication is increased from 2.0 s to 18.6 s as the number of processes increases from 2 to 128.

We also evaluate the performance of the simulator with weak scaling using a scaled workload. We fix the number of routers on one process to 8 and increase the size of the network with the number of processes. Figure 8(c) illustrates the execution time and efficiency of the simulations. We observe that the execution time increases as the model gets larger and runs on more processes. The majority of overhead results from the synchronization algorithm; little overhead is generated from the socket communication as local QSMs process more than 90% of requests without socket communication. We explore the cause of this increased synchronization overhead below. Note that we scale the workload by increasing the size of the network; therefore, the workload does not exactly increase linearly and thus causes the increased time (from 40.9 to 60.7 s) of event execution. In addition, the parallelization efficiency decreases with the number of processes resulting from the increased overhead. The efficiency of the 2-process simulation is 0.96, whereas the efficiency of the 128-process simulation is 0.19.

Both strong-scaling and weak-scaling evaluations show an increasing overhead. We observe that the performance in both cases drops significantly when the number of processes increases from 32 to 64. There are two reasons for this drop. First, the performance of the multithreaded global QSM drops when the number of processes is larger than 32. Since we execute the global QSM on one computer that serves all MPI processes, the capability of the global QSM is limited by the number of processors and is affected by the MPI processes on the same machine. When the number of processes is less than or equal to 36, the number of processors on one Broadwell node, the global QSM threads do not compete with each other. However, when the number of processes is greater than 36, threads of global QSM compete for computation resources and thus show worse performance. We observe that the global QSM achieves a similar processing speed – defined as the time to complete a communication request with a response from the global QSM – for parallel simulation with up to 32 processes. The processing time of the global QSM is doubled for simulations with 64 processes and tripled for 128 processes. The degraded performance is also demonstrated with the growing socket communication time as the number of processes increases from 32 to 128 (see Figure 8(c)).

The second overhead is due to the unbalanced workload. The red bars in Figure 8 indicate the synchronization overhead including the unbalanced workload and the transmission of data. Therefore, we design further experiments to distinguish and quantify which one actually results in the increased synchronization time. The experiments duplicate the cross-process data before exchanging events. When the processes receive data from other processes, the duplicated data is discarded. These experiments increase the amount of cross-process transmission during synchronization without affecting the distribution of the workload. Figure 8(d) shows the execution time of simulations that duplicate data from 0 (i.e., no duplication) to 8 times. The four simulations show a similar time in terms of not only event execution and socket communication but also synchronization. Therefore, we conclude that the unbalanced workload during the time window causes a large overhead even if the aggregated workload is balanced. Since the workload balance depends on the simulation setup, addressing this kind of overhead is out of the scope of the design of a parallel simulator.

We also evaluate the lookahead value improvement technique with both strong scaling and weak scaling. Figure 8(e) and (f) illustrate the performance of the simulator with and without this technique, respectively. For parallel simulation, the left bar represents the larger lookahead value (\(T_{cc}/2\)), and the right bar represents the smaller lookahead value (\(T_{qc}\)). As a result, the lookahead value is increased from 2.5 \({\mu }\)s to 150 \({\mu }\)s. This proposed technique reduces the simulation synchronization overhead in all cases and thus improves the scalability of the simulator. We also verify that the event traces, including the number of executed simulation events, timestamp, object and function names, and arguments of the function, between the parallel and sequential simulations are identical in all cases.⁴

6.2 Quantum Network with Caveman Graph Topology

The second network uses the caveman graph topology [49] to form the quantum network. The caveman graph \(Caveman(n,k)\) is formed from a set of n isolated k-node caves, as a clique, by removing one edge from each cave and using it to connect a neighboring cave along a cycle. The community cave structure of the caveman graph allows us to easily group nodes into sets of nodes for parallelization. Nodes in the same cave are assigned to the same process, and neighboring caves are placed together. Compared with networks in Section 6.1, this network serves multiple flows with a smaller number of hops. Figure 9(a) illustrates the network using the caveman graph \(Caveman(4,4)\) as the topology, which denotes that both the size of one cave and the number of caves are 4.

All simulations use the \(Caveman(128,8)\) network topology that serves 1,024 random flows. Every node in the caveman graph represents one quantum router, and every edge in the graph represents optical fibers and the BSM node used for the meet-in-the-middle protocol. The distance between adjacent routers is 10 km. A sufficient number of quantum memories is allocated on each node to support all flows that pass through them. In our setup, 81,500 quantum memories are used in total. At the beginning of the simulation, every router produces a flow that occupies 25 memories in the source and destination routers of the flow and 50 memories in all intermediate routers. The destination routers are selected randomly based on the number of hops in the flow. The number of hops follows an exponential distribution with \(\lambda = 1\). We assume that the number of short-distance flows is greater than the number of long-distance flows in a realistic quantum network.

Fig. 9.

View Figure

Figure 9(b) illustrates the performance of simulations executed on 1 to 128 processes, with execution time shown as an average over five runs. Again, we note that the standard deviation in execution time ranges from only 6.0 to 164.0 s and thus error bars are not plotted. We evaluate the scalability of the parallel simulation by inspecting the execution time and the speedup with strong scaling. We observe that the execution time decreases as the number of processes grows. The parallel simulation reduces the execution time by a factor of up to 20.5. Although the speedup rate of simulating this network is similar to the speedup rate shown in Section 6.1 in scenarios with fewer than 64 processors, the speedup rate decreases around 23% for parallel simulations with 64 and 128 processors. The scalability degradation results from an unbalanced workload. Here we use the number of quantum memories utilized on each process to quantify load imbalance on MPI processes. We calculate the coefficient of variation (CV), which divides the standard deviation by the mean. We observe that the value of CV increases from 0.01 (2 processes) to 0.17 (128 processes). The unbalanced distribution of quantum memories causes an unbalanced workload, which increases the overhead and therefore reduces the scalability of parallelization. However, placing routers in the same cluster within one process significantly reduces the amount of cross-process entanglement. Less than 10% of work is forwarded to the global QSM. Even if the 128-process simulation has the most cross-process entanglement states, the simulation incurs only a 30-s wait time (of a 240-s total runtime) due to responses from the socket communication. The partition scheme used for caveman networks balances the number of nodes in processes and reduces the number of cross-process entanglement states, but the unbalanced distribution of quantum memories still reduces the scalability. Therefore, we use the next simulation to investigate how partition schemes affect the scalability of parallel simulation on a more complicated topology.

6.3 Quantum Network with Autonomous System Network Topology

Fig. 10.

View Figure

We choose the graph resembling the Internet autonomous system (AS) network as the topology of the third simulated network. As with previous simulations, nodes of the graph represent quantum routers, and edges represent optical fibers with BSM nodes used for the meet-in-the-middle protocol. These routers follow the same method introduced in Section 6.2 to produce random flows. We simulate the quantum network with 1,024 routers in the AS network topology. There are 81,300 quantum memories used (again, based on flow usage) to serve 1,024 random flows. We investigate the scalability of the parallel simulator under this scenario. We use Networkx [19] to generate the graph \(AS(1024,0)\) that denotes a random 1,024-node graph using seed 0.

Figure 10(a) illustrates a smaller network \(AS(32,1)\). Unlike the previous two network topologies, this network topology includes many edge routers that are connected by the core routers. These core routers have a higher degree than do the edge routers, which results in an unbalanced distribution of workload among the routers. The unbalanced workload on the routers is a challenging scenario for the network partition algorithm for parallel simulation.

We propose three heuristic solutions for network partitioning: \(Partition 1\), \(Partition 2\), and \(Partition 3\). The proposed solutions use simulated annealing [47] to optimize their energy functions. We define three required components for the simulated annealing: state, neighboring state, and energy function. State defines the mapping relation between nodes and processes. The initial state evenly assigns routers to processes. Neighboring state is the state that could be generated by swapping two routers in different processes. The initial state and neighboring state make the partition scheme always distribute nodes to processes evenly. The simulated annealing method starts from the initial state and randomly jumps between neighboring states. After a few iterations, the state with the lowest energy is the final partition scheme. Given different energy functions, simulated annealing may produce different partition schemes. We design three energy functions using different assumptions.

The energy functions of \(Partition 1\) and \(Partition 2\) are designed to reduce the overhead of the socket communication, while the energy function of \(Partition 3\) is designed to reduce the overhead of unbalanced workloads. We conduct parallel simulations with these three partition solutions on 2 to 128 processes. These simulations help us explore how various network partitions affect the parallelization performance.

Figure 10(b) illustrates the execution time of the simulations. We note again that the maximum standard deviation in simulation time is 220.8 s, so error bars are not plotted. For parallel simulations, the left, middle, and right bars represent partitioning using \(Partition 1\), \(Partition 2\), and \(Partition 3\), respectively. We observe that parallel simulations using \(Partition 3\) consume more time for socket communication than the other two partition methods because of the different objective functions. Figure 10(c) demonstrates that the increased socket communication time results from the additional requests processed by the global QSM. Compared with \(Partition 1\), \(Partition 3\) requires the local QSM to forward 30% more requests to the global QSM. Among the three methods, \(Partition 1\) performs the best when the number of processes is 2, 4, 64, and 128; \(Partition 2\) performs the worst in most cases; \(Partition 3\) performs the best when the number of processes is 8, 16, and 32. When the parallel simulation uses fewer processes, \(Partition 1\) finds a partition with low overhead from the socket communication, and the imbalance of workload is suppressed by the limited number of processes as well as the equal number of nodes in the processes. As the number of processes increases, the overhead from the unbalanced workload dominates the execution time, and the socket communication time is mitigated by the multithreaded global QSM. Thus, \(Partition 3\) shows the best parallel performance, exemplified by its measured memory CV for 32 processes of 0.10 compared to 0.37 and 0.43 for \(Partition 1\) and \(Partition 2\).⁵ When the number of processes is too large, however, equally distributing nodes to processes becomes a bottleneck in \(Partition 3\). For example, with a parallel simulation of 128 processes, we expect every process to control \(81300 / 128 \approx 89\) memories on average, but the heaviest router controls 2,675 memories. Seven routers need to be placed in the same process as the heaviest router, which deteriorates the workload balance. In addition to the decreased overhead on socket communication, the three partition solutions show similar performance when the number of processes exceeds 32.

\(Partition 3\) demonstrates the largest socket communication time. Therefore, those scenarios are good candidates for evaluating the proposed techniques to reduce communication time, as described in Section 4.2. We compare the performance of the parallel simulation with work offloading and requests grouping, respectively. Figure 10(d) illustrates the normalized execution time of the simulation. We normalize the execution time of parallel simulations, with normalized time 1 corresponding to no communication time optimization and lower normalized time being desirable. We observe that both optimization techniques always speed up simulations. Offloading work from the global QSM to the local QSMs has a bigger impact on simulations with fewer processes. Since the simulation using fewer processes leads to less cross-process entanglement, more manipulations of quantum states can be handled locally instead of using the global QSM. Compared with work offloading, grouping multiple requests to one message achieves a more significant performance improvement, especially for simulations with more processes. By grouping requests, parallel simulations can save more than 75% of the communication time. Moreover, the benefit from this optimization is less affected by the increased number of processes. After aggregating both optimization methods, we speed up the simulation using 128 processes by 4.5 times and the simulation using 2 processes by 11.5 times, which significantly improves the scalability of our simulator.