NetSquid, a NETwork Simulator for QUantum Information using Discrete events

In order to bring quantum networks into the real world, we would like to determine the requirements of quantum network protocols including the underlying quantum hardware. Because detailed architecture proposals are generally too complex for mathematical analysis, it is natural to employ numerical simulation. Here we introduce NetSquid, the NETwork Simulator for QUantum Information using Discrete events, a discrete-event based platform for simulating all aspects of quantum networks and modular quantum computing systems, ranging from the physical layer and its control plane up to the application level. We study several use cases to showcase NetSquid’s power, including detailed physical layer simulations of repeater chains based on nitrogen vacancy centres in diamond as well as atomic ensembles. We also study the control plane of a quantum switch beyond its analytically known regime, and showcase NetSquid’s ability to investigate large networks by simulating entanglement distribution over a chain of up to one thousand nodes. Implementing large-scale quantum networks is one of the challenges at the core of quantum communication. Here, the authors present NetSquid, a quantum network simulator that allows studying how such networks can be built, including physical hardware modelling, modularity, scalability, and examples.


I. INTRODUCTION
Quantum communication can be used to connect distant quantum devices into a quantum network.At short distances, networking quantum devices provides a path towards a scalable distributed quantum computer [1].At larger distances, interconnected quantum networks allow for communication tasks between distant users on a quantum internet.Both types of quantum networks have the potential for large societal impact.First, analogous to classical computers, it is likely that any approach for scaling up a quantum computer so that it can solve real world problems impractical to treat on a classical computer, will require the interconnection of different modules [2][3][4].Furthermore, quantum communication networks enable a host of tasks that are impossible using classical communication [5].For both types of networks, many challenges must be overcome before they can fulfil their promise.The exact extent of these challenges remains unknown, and precise requirements to guide the construction of large-scale quantum networks are missing.At the physical layer, many proposals exist for quantum repeaters that can carry qubits over long distances (see e.g.[6][7][8] for an overview).Using analytical methods  and ad-hoc simulations [31][32][33][34][35][36][37][38] rough estimates for the requirements of such hardware proposals have been obtained.Yet, while greatly valuable to set minimal requirements, these studies still provide limited detail.Even for a smallscale quantum network, the intricate interplay between many communicating devices, and the resulting timing dependencies makes a precise analysis challenging.To go beyond, we would like a tool that can incorporate not only a detailed physical modelling, but also account for the implications of timedependent behaviour.Quantum networks cannot be built from quantum hardware alone; in order to scale they need a tightly integrated classical control plane, i.e. protocols responsible for orchestrating quantum network devices to bring entanglement to users.Fundamental differences between quantum and classical information demand an entirely new network stack in order to create entanglement, and run useful applications on future quantum networks [39][40][41][42][43][44].The design of such a control stack is furthermore made challenging by numerous technological limitations of quantum devices.A good example is given by the limited lifetimes of quantum memories, due to which delays in the exchange of classical control messages have a direct impact on the performance of the network.To succeed, we hence need to understand how possible classical control strategies do perform on specific quantum hardware.Finally, to guide overall development, we need to understand the requirements of quantum network applications themselves.Apart from quantum key distribution (QKD) [45][46][47][48][49] and a few select applications [50][51][52][53], little is known about the requirements of quantum applications [5] on imperfect hardware.Analytical tools offer only a limited solution for our needs.Statistical tools (see e.g.[54][55][56][57]) have been used to make predictions about certain aspects of large regular networks using simplified models, but are of limited use for more detailed studies.Information theory [58] can be used to benchmark implementations against the ideal performance.However, it gives no information about how well a specific proposal will perform.As a consequence, numerical methods are of great use to go beyond what is feasible using an analytical study.Ad-hoc simulations of quantum networks have indeed been used to provide further insights on specific aspects of quantum networks (see e.g.[31][32][33][34][35][36][37][38][59][60][61]).However, while greatly informative, setting up ad-hoc simulations for each possible networking scenario to a level of detail that might allow the determination of more precise requirements is cumbersome, and does not straightforwardly lend itself to extensive explorations of new possibilities.We would hence like a simulation platform that satisfies at least the following three features: First, accuracy: the tool should allow modelling the relevant physics.This includes the ability to model timedependent noise and network behaviour.Second, modularity: it should allow stacking protocols and models together in order to construct complicated network simulations out of simple components.This includes the ability to investigate not only the physical layer hardware, but the entirety of the quantum network system including how different control protocols behave on a given hardware setup.Third, scalability: it should allow us to investigate large networks.Evaluating the performance of large classical network systems, including their time-dependent behaviour is the essence of classical network analysis.Yet, even for classical networks, the predictive power of analytical methods is limited due to complex emergent behaviour arising from the interplay between many network devices.Consequently, a crucial tool in the design of such networks are network simulators, which form a tool to test new ideas, and many such simulators exist for purely classical networks [62][63][64].However, such simulators do not allow the simulation of quantum behaviour.In the quantum domain, many simulators are known for the simulation of quantum computers (see e.g.[65]).However, the task of simulating a quantum network differs greatly from simulating the execution of one monolithic quantum system.In the network, many devices are communicating with each other both quantumly and classically, leading to complex stochastic behaviour, and asynchronous and time-dependent events.From the perspective of building a simulation engine, there is also an important difference that allows for gains in the efficiency of the simulation.A simulator for a quantum computation is optimised to track large entangled states.In contrast, in a quantum network the state space grows and shrinks dynamically as qubits get measured or entangled, but for many protocols, at any moment in time the state space describing the quantum state of the network is small.We would thus like a simulator capable of exploiting this advantage.In this paper we introduce the quantum network simulator NetSquid: the NETwork Simulator for QUantum Information using Discrete events.Net-Squid is a software tool (available as a package for Python and previously made freely available online [66]) for accurately simulating quantum networking and modular computing systems that are subject to physical non-idealities.It achieves this by integrating several key technologies: a discrete-event simulation engine, a specialised quantum computing library, a modular framework for modelling quantum hardware devices, and an asynchronous programming framework for describing quantum protocols.We showcase the utility of this tool for a range of applications by studying several use cases: the analysis of a control plane protocol beyond its analytically accessible regime, predicting the performance of protocols on realistic near-term hardware, and benchmarking different quantum devices.These use cases, in combination with a scalability analysis, demonstrate that NetSquid achieves all three features set forth above.Furthermore, they show its potential as a general and versatile design tool for quantum networks, as well as for modular quantum computing architectures.

A. NetSquid in a nutshell
Simulating a quantum network with NetSquid is generally performed in three steps.Firstly, the network is modelled using a modular framework of components and physical models.Next, protocols are assigned to network nodes to describe the intended behaviour.Finally, the simulation is executed for a typically large number of independent runs to collect statistics with which to determine the performance of the network.To explain these steps and the features involved further, we consider a simple use case for illustration.For a more detailed presentation of the available functionality and design of the NetSquid framework see section Design and functionality of NetSquid of the Methods.The scenario we will consider is the analysis of an entanglement distribution protocol over a quantum repeater chain with three nodes.The goal of the analysis is to estimate the average output fidelity of the distributed entangled pairs.The entanglement distribution protocol is depicted in Figure 1(d-e).It works as follows.First, the intermediate node generates two entangled pairs with each of its adjacent neighbours.Entanglement generation is modelled as a stochastic process that succeeds with a certain probability at every attempt.When two pairs are ready at one of the links, the DEJMPS entanglement distillation scheme [67] is run to improve the quality of the entanglement.If it fails, the two links are discarded and the executing nodes restart entanglement generation.Once both distilled states The state of each qubit is updated as it is accessed during the simulation, for instance to apply time-dependent noise from waiting in memory.e) A zoom in of the distillation protocol.The shared quantum states of the qubits are combined in an entangling step, which then shrinks as two of the qubits are measured.The output is randomly sampled, causing the simulation to choose one of two paths by announcing success or failure.f ) A plot illustrating the stochastic paths followed by multiple independent simulation runs over time, labeled by their final end-to-end fidelity F i .The blue dashed line corresponds to the run shown in panel (d).The runs are typically executed in parallel.Their results are statistically analysed to produce performance metrics such as the average outcome fidelity and run duration.
are ready, the intermediate node swaps the entanglement to achieve end-to-end entanglement.We remark that already this simple protocol is rather involved to analyse.
We begin by modelling the network.The basic element of NetSquid's modular framework is the "component".It is capable of describing the physical model composition, quantum and classical communication ports, and, recursively, any subcomponents.All hardware elements, including the network itself, are represented by components.For this example we require three remote nodes linked by two quantum and two classical connections, the setup of which is shown in Figure 1(a).In Figure 1(b,c) the nested structure of these components is highlighted.A se-lection of physical models is used to describe the loss and delay of the fibre optic channels, the decoherence of the quantum memories, and the errors of quantum gates.
Quantum information in NetSquid is represented at the level of qubits, which are treated as objects that dynamically share their quantum states.These internally shared states will automatically merge or "split" -a term we use to mean the separation of a tensor product state into two separately shared sub-states -as qubits entangle or are measured, as illustrated by the distillation protocol in Figure 1 Discrete-event simulation, an established method for simulating classical network systems [68], is a modelling paradigm that progresses time by stepping through a sequence of events -see Figure 2 for a visual explanation.This allows the simulation engine to efficiently handle the control processes and feedback loops characteristic of quantum networking systems, while tracking quantum state decoherence based on the elapsed time between events.A novel requirement for its application to quantum networks is the need to accurately evolve the state of the quantum information present in a network with time.This can be achieved by retroactively updating quantum states when the associated qubits are accessed during an event.While it is possible to efficiently track a density matrix, quantum operations requiring a singular outcome for classical decision making, for instance a quantum measurement, must be probabilistically sampled.A single simulation run thus consists of a sequence of random choices that forms one of many possible paths.In Figure 1 (d) we show such a run for the repeater protocol example, which demonstrates the power of the discrete-event approach for tracking qubit decoherence and handling feedback loops.
The performance metrics of a simulation are determined statistically from many runs.Due to the independence of each run, simulations can be massively parallelised and thereby efficiently executed on computing clusters.For the example at hand we choose as metrics the output fidelity and run duration.In Figure 1 (f) the sampled run from (d), which resulted in perfect fidelity, is plotted in terms of its likelihood and duration together with several other samples, some less successful.By statistically averaging all of the sampled runs the output fidelity and duration can be estimated.
In the following sections, we will outline three use cases of NetSquid: first, a quantum switch, followed by simulations of quantum repeaters based on nitrogen-vacancy technology or atomic-ensemble memories.We will also benchmark NetSquid's scalability in both quantum state size and number of quantum network nodes.Although the use cases each provide relevant insights into the performance of the studied hardware and protocols, we emphasise that one can use NetSquid to simulate arbitrary network topologies.Abstract example of simulating a quantum protocol with discrete events.When setting up the simulation, protocol actions are defined to occur when a specific event occurs, as in the setup of a real system.Instead of performing a continuous time evolution, the simulation advances to the next event, and then automatically executes the actions that should occur when the event takes place.Any action may again define future events to be triggered after a certain (stochastic) amount of time has elapsed.
For concreteness a simplified quantum teleportation example is shown where a qubit, shown as an orange circle with arrow, is teleported between the quantum memories of Alice and Bob.Here, entanglement is produced using an abstract source sending two qubits, shown as blue circles with arrows, to Alice and Bob.Once the qubit has traversed the fibre and reaches Alice's lab, an event is triggered that invokes the simulation of Alice's Bell state measurement (BSM) apparatus.The simulation engine steps from event to events defined by the next action, which generally occur at irregular intervals.This approach allows time-dependent physical non-idealities, such as quantum decoherence, to be accurately tracked.

B. Simulating a quantum network switch beyond its analytically known regime
As a first use case showcasing the power of NetSquid, we study the control plane of a recently introduced quantum switch beyond the regime for which analytical results have been obtained, including its performance under time-dependent memory noise.The switch is a node which is directly connected to each of k users by an optical link.The communications task is distributing Bell pairs and n-partite Greenberger-Horne-Zeilinger (GHZ) states [69] between n ≤ k users.The switch achieves this by connecting Bell pairs which are generated at random intervals on each link.See Figure 3.
Intuitively, the switch can be regarded as a generalisation of a simple repeater performing entanglement swapping with added logic to choose which parties to link.Even with a streamlined physical model, it is already rather challenging to analytically characterise the switch use case [56].
In the following, we recover via simulation a selection of the results from Vardoyan et al. [56], who studied the switch as the central node in a star network, and extend them in two directions.First, we increase the range of parameters for which we can estimate entanglement rates using the same model as used in the work of Vardoyan et al.Second, simulation enables us to investigate more sophisticated models than the exponentially distributed erasure process from their work, in particular we analyse the behaviour of a switch in the presence of memory dephasing noise.
The protocol for generating the target n-partite GHZ states is simple.Entanglement generation is attempted in parallel across all k links.If successful they result in bipartite Bell states that are stored in quantum memories.The switch waits until n Bell pairs have been generated until performing an npartite GHZ measurement, which converts the pairs into a state locally equivalent to a GHZ state.An additional constraint is that the switch has a finite buffer B of number of memories dedicated for each user (see Figure 3).If the number of pairs stored in a link is B and a new pair is generated, the old one is dropped and the new one is stored.
The protocol can be translated to a Markov chain.The state space is represented by a k-length vector where each entry is associated with a link and its value denotes the number of stored links.The switch's mean capacity, i.e. the number of states produced per second, can be derived from the steadystate of the Markov chain [56].
Using NetSquid, it is straightforward to fully reproduce the previous model and study the behaviour of the network without constructing the Markov Chain (details can be found in Supplementary Note 3).In Figure 4(a), we use NetSquid to study the capacity of a switch network serving nine users.The figure shows the capacity (number of produced GHZstates per second), which we investigate for three use cases.First we consider a switch network distributing bipartite entanglement.Second, we consider also a switch-network serving bipartite entanglement but with link generation rates that do not satisfy the stability condition for the Markov Chain if the buffer B is infinitely large, i.e. a regime so far intractable.Third, we consider a switch-network distributing four-partite entanglement where the link generation rates µ differ per user, a regime not studied so far, and compute the capacity.
Beyond rate, it is important to understand the quality of the states produced.Answering this question with Markov chain models seems challenging.In order to analyse entanglement quality, we introduce a more sophisticated decoherence model where the memories suffer from decay over time.In particular, we model decoherence as exponential T 2 noise, which impacts the quality of the state, as expressed in its fidelity with the ideal state.In Figure 4(b), we show the effect of the time-dependent memory noise on the average fidelity.[56].The switch (central node) is connected to a set of users (leaf nodes) via an optical fibre link that distributes perfect Bell pairs at random times, following an exponential distribution with mean rate µ ∝ e −βL , where L denotes the distance of the link and β the attenuation coefficient.Associated with each link there is a buffer that can store B qubits at each side of the link.As soon as n Bell pairs with different leaves are available, the switch performs a measurement in the n-partite Greenberger-Horne-Zeilinger (GHZ) basis, which results in an n-partite GHZ state shared by the leaves.The GHZ-basis measurement consists of: first, controlled-X gates with the same qubit as control; next, a Hadamard (H) gate on the control qubit; finally, measurement of all qubits individually.The figure shows 4 leaf nodes, GHZ size n = 3 and a buffer size B = 2.

C. Sensitivity analysis for the physical modelling of a long range repeater chain
The next use case is the distribution of long-distance entanglement via a chain of quantum repeater nodes [6,9] based on nitrogen-vacancy (NV) centres in diamond [70,71].This example consists of a more detailed physical model and more complicated control plane logic than the quantum switch or the distillation example presented at the start of this section.It is also an example of how NetSquid's modularity supports setting up simulations involving many nodes; in this case the node model and the protocol (which runs locally at a node) only need to be specified once, and can then be assigned to each node in the chain.Furthermore, the use of a discrete-event engine allows the actions of the individual protocols to be simulated asynchronously, in contrast to the typically sequential execution of quantum comput-  Capacity as a function of the buffer size (number of quantum memories that the switch has available per user) for either 2− or 4−qubit Greenberger-Horne-Zeilinger (GHZ)-states.For each scenario, the generation rate µ of pairs varies per user.For the blue scenario (2-partite entanglement, µ = [1.9,1.9, 1.9, 1, 1, 1, 1, 1, 1] MHz), the capacity was determined analytically by Vardoyan et al. using Markov Chain methods [56, Figure 8].Here we extend this to 4-partite entanglement (orange scenario, same µs), for which Vardoyan et al. have found an upper bound (by assuming unbounded buffer and each µ = maximum of original rates = 1.9 MHz) but no exact analytical expression.The green scenario (µ = [15, 1.9, 1.9, 1, does not satisfy the stability condition for the Markov chain for unbounded buffer size (each leaf's rate < half of sum of all rates) so in that case steady-state capacity is not well-defined.We note that regardless of buffer size, the switch has a single link to each user, which is the reason why the capacity does not scale linearly with buffer size.(b) Average fidelity of the produced entanglement on the user nodes (no analytical results known) with unbounded buffer size.The fact that the green curve has lower fidelity than the blue one, while the former has higher rates, can be explained from the fact that the protocol prioritises entanglement which has the longest storage time (see Supplementary Note 3).Each data point represents the average of 40 runs (each 0.1 ms in simulation).Standard deviation is smaller than dot size.
The NV-based quantum processor includes the following three features.First, the nodes have a single communication qubit, i.e. a qubit acting as the optical interface that can be entangled with a remote qubit via photon interference.This seemingly small restriction has important consequences for the communications protocol.In particular, entanglement can not proceed in parallel with both adjacent nodes.As a consequence, operations need to be scheduled in sequence and the state of the communication qubit transferred onto a storage qubit.Second, the qubits in a node are connected with a star topology with the communication qubit located in the centre.Two-qubit gates are only possible between the communication qubit and a storage qubit.Third, communication and storage qubits have unequal coherence times.Furthermore, the storage qubits suffer additional decoherence when the node attempts to generate entanglement.Previous repeater-chain analyses, e.g.[22,23,43], did not take all three into account simultaneously.
Together with the node model, we consider two protocols: swap-asap and nested-with-distill.In swap-asap, as soon as adjacent links are generated the entanglement is swapped.nested-withdistill is a nested protocol [9] with entanglement distillation at every nesting level.For a description of the simulation, including the node model and protocols, see Methods, section Implementing a processing-node repeater chain in NetSquid.
The first question that we investigate is the distance that can be covered by a repeater chain.For this we choose two sets of hardware parameters that we dub near-term and 10× improved (see Supplementary Note 4) and choose two configurations: one without intermediate repeaters and one with three of them.We observe, see Figure 5(a), that the repeater chain performs worse in fidelity than the repeaterless configuration with near-term hardware.For improved hardware, we see two regimes, for short distances the use of repeaters increases rate but lowers fidelity while from 750 km until 1500 km the repeater chain outperforms the no-repeater setup.
The second question that we address is which protocol performs best for a given distance.We consider seven protocols: no repeater, and repeater chains implementing swap-asap or nested-withdistill over 1, 3 or 7 repeaters.The latter is motivated by the fact that the nested-withdistill protocol is defined for 2 n − 1 repeaters (n ≥ 1), and thus 1, 3, and 7 are the first three possible configurations.In Figure 5(b), we sweep over the hardware parameter space for two distances, where we improve all hardware parameters simultaneously and the improvement is quantified by a number we refer to as "improvement factor" (see section How we choose improved hardware parameters of the Methods).For 500 km, we observe that the no-repeater configuration achieves larger or equal fidelity for the entire range studied.However, repeater schemes boost the rate for all parameter values.If we increase the distance to 800 km, then we see that the use of repeaters increases both rate and fidelity for the same range of parameters.If we focus on the repeater scheme, we observe for both distances that for high hardware quality, the nested-with-distill scheme, which includes distillation, is optimal.In contrast, for lower hardware quality, the best-performing scheme that achieves fidelities larger than the classical bound 0.5 is the swap-asap protocol.We note that beyond 700 km the entanglement rate decreases when the hardware is improved.This is due to the presence of dark counts, i.e. false signals that a photon has been detected.At large distances most photons dissipate in the fibre, whereby the majority of detector clicks are dark counts.Because a dark count is mistakenly counted as a successful entanglement generation attempt, improving (i.e.decreasing) the dark count rate in fact results in a lower number of observed detector clicks, from which the (perceived) entanglement rate plotted in Figure 5(a) is calculated.Lastly, in Figure 6, we investigate the sensitivity of the entanglement fidelity for the different hardware parameters.We take as the figure of merit the best fidelity achieved with a swap-asap protocol.
The uniform improvement factor is set to 3, while the following four hardware parameters are varied: a two-qubit gate noise parameter, photon detection probability (excluding transmission), induced storage qubit noise and visibility.We observe that improving the detection probability yields the largest fidelity increase from 2× to 50× improvement, while this increase is smallest for visibility.We also see that improving two-qubit gate noise or induced storage qubit noise on top of an increase in detection probability yields only a small additional fidelity improvement, which however boosts fidelity beyond the classical threshold of 0.5.These observations indicate that detection probability is the most important parameter for realising remote-entanglement generation with the swap-asap scheme, followed by two-qubit gate noise and induced storage qubit noise.

D. Performance comparison between two atomic-ensemble memory types through NetSquid's modular design
Finally, we showcase that NetSquid's modular design greatly reduces the effort of assessing possible hardware development scenarios.We demonstrate the power of this modularity by simulating point-topoint remote-entanglement generation based on either of two types of atomic-ensemble based quantum memories: atomic frequency combs (AFC) [72] and electronically induced transparency (EIT) [73,74] memories.Both simulations are identical except for the choice of a different quantum memory component.
The two types of memories are a promising building block for high-rate remote entanglement generation through quantum repeaters because of their high efficiency (EIT) or their ability for multiplexing (AFC), i.e. to perform many attempts at entanglement generation simultaneously without network components having to wait for the arrival of classical messages that herald successful generation.
The first type of memories, AFCs, are based on a photon-echo process, where an absorbed photon is re-emitted after an engineered duration.In contrast, the second type, EITs, emit the photon after an ondemand interval, due to optical control.In principle the AFC protocol can be extended to offer such ondemand retrieval as well.At this point both technologies are promising candidates and it is not yet clear which outperforms the other and under what circumstances.Atomic-ensemble based repeaters have been analytically and numerically studied before with streamlined physical models [19,75].NetSquid's discreteevent based paradigm allows us to go beyond that by concurrently introducing several non-ideal characteristics.In particular, we include the emission of more than one photon pair, photon distinguishability and time-dependent memory efficiency.Efficiency in this context is the probability that the absorbed photon will be re-emitted.All these characteristics have a significant impact on the performance of the repeater protocol.
In order to compare the two memory types, we simulate many rounds of the BB84 quantum key distribution protocol [76] between two remote nodes, using a single repeater positioned precisely in between them.Entanglement generation is attempted in synchronised rounds over both segments in parallel.At the end of each round, the two end nodes measure in the X-or Z-basis, chosen uniformly at random, and the repeater performs a probabilistic linear-optical Bell-state measurement.Upon a successful outcome, we expect correlation between the measurement outcomes if they were performed in the same basis.As a figure of merit we choose the asymptotic BB84 secret-key rate.The results of our simulations are shown in Figure 7, where the rate at which secret key between the two nodes can be generated is obtained as a function  Fidelity and entanglement distribution rate achieved with near-term and 10× improved hardware (Supplementary Note 4) with the swap-asap protocol.Dashed line represents classical fidelity threshold of 0.5.We observe that for near-term hardware, the use of 3 repeaters yields worse performance in terms of fidelity than the no-repeater setup.For improved hardware we observe (i) that for approx.0 -750 kms, repeaters improve upon rate by orders of magnitude while still producing entanglement (fidelity> 0.5), while (ii) for approx.750 -1500 kms, repeaters outperform in both rate and fidelity.(b-c) Fidelity and rate achieved without and with repeaters (1, 3 or 7 repeaters) as function of a hardware improvement factor (Methods, section How we choose improved hardware parameters) for two typical distances from both distance regime (i) and (ii), for two protocols swap-asap and nested-with-distill.For the repeater case, only the best-performing number-of-repeaters & protocol in terms of achieved fidelity is shown in (b), accompanied by its rate in (c).Each data point represents the average over (a) 200 and (b) 100 runs.Standard deviation is smaller than dot size.
of the distance between the nodes.For the parameters considered (see Supplementary Note 7), we observe that EIT memories outperform AFC memories at short distances.The crossover performance point is reached at ∼ 50 kilometers, beyond which AFC memories outperform EIT memories.
In the use case above, we showcased NetSquid's modularity by only replacing the memory component.We emphasise that this modularity also applies to different parts of the simulation.For example, if the quantum switch should produce a different type of multipartite state than GHZ states, then one only needs to change the circuit at the switch node.A different example is the NV repeater chain, where one could replace the protocol module (currently either swap-asap or nested-with-distill).

E. Fast and scalable quantum network simulation
NetSquid has been designed and optimised to meet several key performance criteria: to be capable of accurate physical modelling, to be scalable to large networks, and to be sufficiently fast to support multi-variate design analyses with adequate statistics.While it is not always possible to jointly satisfy all the criteria for all use cases, NetSquid's design allows the user to prioritise them.We proceed to benchmark NetSquid to demonstrate its capabilities and unique strengths for quantum network simula-tion.

Benchmarking of quantum computation
To accurately model physical non-idealities, it is necessary to choose a representation for quantum states that allows a characterisation of general processes such as amplitude damping, general measurements, or arbitrary rotations.NetSquid provides two representations, or "formalisms", that are capable of universal quantum computation: ket state vectors (KET) and density matrices (DM), both stored using dense arrays.The resource requirements for storage in memory and the computation time associated with applying quantum operations both scale exponentially with the number of qubits.While the density matrix scales less favourably, 2 2n versus 2 n for n qubits, its ability to represent mixed states makes it more versatile for specific applications.Given the exponential scaling, these formalisms are most suitable for simulations in which a typical qubit lifetime involves only a limited number of (entangling) interactions.
When scaling to large network simulations it can happen that hundreds of qubits share the same entangled quantum state.For such use cases, we need a quantum state representation that scales subexponentially in time and space.NetSquid provides two such representations based on the stabiliser state formalism: "stabiliser tableaus" (STAB) and (Methods, section How we choose improved hardware parameters).Then, for each of the four parameters only, we individually decrease their improvement factor to 2, or increase it to 10 or 50.
The figure shows the resulting fidelity (horizontal and vertical grid lines; dashed line indicates maximal fidelity which can be attained classically).
Note that at an improvement factor of 3 (orange line), all 15 parameters are improved by 3 times, resulting in a fidelity of 0.39.In addition, we vary the improvement factor for combinations of two of the four parameters (diagonal lines).The 3× improved parameter values can be found in Supplementary Table II.The other values (at 2/10/50×) are approximately: two-qubit gate fidelity F EC (0.985/0.997/0.9994),detection probability p nofibre det (6.8%/58%/90%), induced storage qubit noise N 1/e (2800/14000/70000), visibility V (95%/99%/99.8%).The fidelities shown are obtained by simulation of the swap-asap protocol (3 repeaters) with a total spanned distance of 500 km.Each data point represents the average of 1000 runs (standard deviation on fidelity < 0.002).
"graph states with local Cliffords" (GSLC) [77,78] that the user can select.Stabiliser states are a subset of quantum states that are closed under the application of Clifford unitaries and single-qubit measurement in the computational basis.In the context of simulations for quantum networks stabiliser states are particularly interesting because many network protocols consist of only Clifford operations and noise can be well approximated by stochastic application of Pauli gates.For a theoretical compar-ison of the STAB and GSLC formalisms see Supplementary Note 1.
The repetitive nature of simulation runs due to the collection of statistics via random sampling allows NetSquid to take advantage of "memoization" for expensive quantum operations, which is a form of caching that stores the outcome of expensive operations and returns them when the same input combinations reoccur to save computation time.Specifically, the action of a quantum operator onto a quantum state for a specific set of qubit indices and other discrete parameters can be efficiently stored, for instance as a sparse matrix.Future matching operator actions can then be reduced to a fast lookup and application, avoiding several expensive computational steps -see the Methods, section Qubits and quantum computation for more details.
In the following we benchmark the performance of the available quantum state formalisms.For this, we first consider the generation of an n qubit entangled GHZ state followed by a measurement of each qubit (see section Benchmarking of the Methods).For a baseline comparison with classical quantum computing simulators we also include the ProjectQ [79] package for Python, which uses a quantum state representation equivalent to our ket vector.We show the average computation time for a single run versus the number of qubits for the different quantum computation libraries in Figure 8(a).The exponential scaling of the universal formalisms in contrast to the stabiliser formalisms is clearly visible, with the density matrix formalism performing noticeably worse.For the ket formalism we also show the effect of memoization, which gives a speed-up roughly between two and five.
Let us next consider a more involved benchmarking use case: the quantum computation involved in simulating a repeater chain i.e.only the manipulation of qubits, postponing all other simulation aspects, such as event processing and component modelling, to the next section.This benchmark involves the following steps: first the N − 1 pairs of qubits along an N node repeater chain are entangled, then each qubit experiences depolarising noise, and finally adjacent qubits on all but the end-nodes do an entanglement swap via a Bell state measurement (BSM).
If the measured qubits are split from their shared quantum states after the BSM, then the size of any state is limited to four qubits.We observe that for the NetSquid formalisms (but not for ProjectQ) keeping qubits "in-place" after each measurement is more performant than "splitting" them below a certain threshold due to the extra overhead of doing the latter.The ket vector formalism is seen to be the most efficient for this benchmarking use case if states are split after measurement.When the measurement operations are performed in-place the GSLC formalism performs  the best beyond 15 qubits.

Benchmarking of event-driven simulations
As explained in the results section, a typical Net-Squid simulation involves repeatedly sampling many independent runs.As such NetSquid is "embarrassingly parallelisable": the reduction in runtime scales linearly with the number of processing cores available, assuming there is sufficient memory available.
Nonetheless, given the computational requirements associated with collecting sufficient statistics and analysing large parameter spaces it remains crucial to optimise the runtime performance per core.
Depending on the size of the network, the detail of the physical modelling, and the duration of the protocols under consideration, the number of events processed for a single simulation run can range anywhere from a few thousand to millions.To efficiently process the dynamic scheduling and handling of events NetSquid uses the discrete-event simulation engine PyDynAA [80] (see section Discrete event simulation of the Methods).NetSquid aims to schedule events as economically as possible, for instance by streamlining the flow of signals and messages between components using inter-connecting ports.
To benchmark the performance of an event-driven simulation run in NetSquid we consider a simple network that extends the single repeater (without distillation) shown in Figure 1 into an N node chain -see Supplementary Note 2 for further details on the simulation setup.For the quantum computation we will use the ket vector formalism based on the benchmarking results from the previous section, and split qubits from their quantum states after measurement to avoid an exponential scaling with the number of nodes.In Figure 9 we show the average computation time for deterministically generating end-to-end entanglement versus the number of nodes in the chain.Also shown is a relative breakdown in terms of the time spent in the NetSquid sub-packages involved, as well as the PyDynAA and NumPy packages.We observe that the biggest contribution to the simulation runtime is the components sub-package, which accounts for 30% of the total at 1000 nodes.The relative time spent in each of the NetSquid subpackages, as well as NumPy and PyDynAA, is seen to remain constant with the number of nodes.The total runtime of each of the NetSquid sub-packages is the sum of many small contributions, with the costliest function for the components sub-package for a 1000 node chain, for example, contributing only 7% to the total.Extending this benchmark simulation with more detailed physical modelling may shift the relative runtime distribution and impact the overall performance.For example, more time may be spent in calls to the "components" and "components.models"subpackages, additional complexity can increase the volume of events processed by the "pydynaa" engine, and extra quantum characteristics can lead to larger quantum states.In case of the latter, however, the effective splitting of quantum states can still allow such networks to scale if independence among physical elements can be preserved.

F. Comparison with other quantum network simulators
Let us compare NetSquid to other existing quantum network simulators.First, SimulaQron [81] and QuNetSim [82] are two simulators that do not aim at realistic physical models of channels and devices, or timing control.Instead, SimulaQron's main purpose is application development.It is meant to be run in a distributed fashion on physically-distinct classical computers.QuNetSim focuses on simplifying the development and implementation of quantum network protocols.
In contrast with SimulaQron and QuNetSim, the simulator SQUANCH [83] allows for quantum network simulation with configurable error models at Other than NetSquid, there now exist three discreteevent quantum simulators: the QuISP [84], qkdX [85] and SeQUeNCe [86] simulators.With these simulators it is possible to accurately characterise complex timing behaviour, however they differ in goals and scope.Similarly to NetSquid, QuISP aims to support the investigation of large networks that consist of too many entangled qubits for full quantum-state tracking.In contrast to NetSquid, which achieves this by managing the size of the state space, and providing the stabiliser representation as one of its quantum state formalisms, QuISP's approach is to track an error model of the qubits in a network instead of their quantum state.qkdX, on the other hand, captures the physics more closely through models of the quantum devices but is restricted to the simulation of quantum key distribution protocols.Lastly, SeQUeNCe, similar to Net-Squid, aims at simulation at the level of hardware, control plane or application.It has a fixed control layer consisting of reprogrammable modules.In contrast, NetSquid's modularity is not tied to a particular network stack design.Furthermore, it is unclear to us how performant SeQUeNCe's quantum simulation engine is: currently, at most a 9-node network has been simulated, whereas NetSquid's flexibility to choose a quantum state representation enables scalability to simulation of networks of up to 1000 nodes.

G. Conclusions
In this work we have presented our design of a modular software framework for simulating scalable quantum networks and accurately modelling the nonidealities of real world physical hardware, providing us with a design tool for future quantum networks.We have showcased its power and also its limitations via example use cases.Let us recap NetSquid's main features.
First, NetSquid allows the modelling of any physical device in the network that can be mapped to qubits.
To demonstrate this we studied two use cases involving nitrogen-vacancy centres in diamond as well as atomic-ensemble based memories.Second, NetSquid is entirely modular, allowing users to set up large scale simulations of complicated networks and to explore variations in the network design; for example, by comparing how different hardware platforms perform in an otherwise identical network layout.Moreover, this modularity makes it possible to explore different control plane protocols for quantum networks in a way that is essentially identical to how such protocols would be executed in the real world.Control programs can be run on any simulated network node, exchanging classical and quantum communication with other nodes as dictated by the protocol.That allows users to investigate the intricate interplay between control plane protocols and the physical devices dictating the performance of the combined quantum network system.
As an example, we studied the control plane of a quantum network switch.NetSquid has also already found use in exploring the interplay between the control plane and the physical layer in [39,87,88].Finally, to allow large scale simulations, the quantum computation library used by NetSquid has been designed to manage the dynamic lifetimes of many qubits across a network.It offers a seamless choice of quantum state representations to support different modelling use cases, allowing both a fully detailed simulation in terms of wave functions or density matrices, or simplified ones using certain stabiliser formalisms.As an example use case, we explored the simulation run-time of a repeater chain with up to one thousand nodes.
In light of the results we have presented, we see a clear application for NetSquid in the broad context of communication networks.It can be used to predict performance with accurate models, to study the stability of large networks, to validate protocol designs, to guide experiment, etc.While we have only touched upon it in our discussion of performance benchmarks, NetSquid would also lend itself well to the study of modular quantum computing architectures, where the timing of control plays a crucial role in studying their scalability.For instance, it might be used to validate the microarchitecture of distributed quantum computers or more generally to simulate different components in modular architectures.

A. Design and functionality of NetSquid
The NetSquid simulator is available as a software package for the Python 3 programming language.
It consists of the sub-packages "qubits", "components", "models", "nodes", "protocols" and "util", which are shown stacked in Figure 10.NetSquid depends on the PyDynAA software library to provide its discrete-event simulation engine [80].Under the hood speed critical routines and classes are written in Cython [89] to give C-like performance, including its interfaces to both PyDynAA and the scientific computation packages NumPy and SciPy.In the following subsections we highlight some of the main design features and functionality of NetSquid; for a more detailed presentation see Supplementary Note 1.The main classes in each (sub-)package are highlighted, and their relationships in terms of inheritance, composition and aggregation are shown.Also shown are the key modules users interact with, which are described in the main text.

Data-Collector
In this paper NetSquid version 0.10 is described.

Discrete event simulation
The PyDynAA package provides a fast, powerful, and lightweight discrete-event simulation engine.It is a C++ port of the core engine layer from the DynAA simulation framework [80], with bindings added for the Python and Cython languages.Dy-nAA defines a concise set of classes and concepts for modelling event-driven simulations.The simulation engine manages a timeline of "events", which can only be manipulated by objects that are subclasses of the "entity" base class.Simulation entities can dynamically schedule events on the timeline and react to events by registering an "event handler" object to wait for event(s) with a specified type, source entity, or identifier to be triggered.
To deal with the timing complexities encountered in NetSquid simulations, an "event expression" class was introduced to PyDynAA to allow entities to also wait on logical combinations of events to occur.Atomic event expressions, which describe regular wait conditions for standard events, can be combined to form composite expressions using logical "and" and "or" operators to any depth.This feature has been used extensively in NetSquid to model both the internal behaviour of hardware components, as well as for programming network protocols.

Qubits and quantum computation
The qubits sub-package of NetSquid defines the "qubit" object that is used to track the flow of quantum information.Qubits internally share quantum state ("QState") objects, which grow and shrink in size as qubits interact or are measured.The "QState" class is an interface that is implemented by a range of different formalisms, as presented in section Benchmarking of quantum computation of the Results and Discussion.Via the qubit-centric API, which provides functions to directly manipulate qubits without knowledge of their shared quantum states, users can program simulations in a formalism agnostic way.Functionality is also provided to automatically convert between quantum states that use different formalisms, and to sample from a distribution of states, which is useful for instance for pure state formalisms.The ket and density matrix formalisms use dense arrays (vectors or matrices, respectively) to represent quantum states.Applying a k qubit operator to an n qubit ket vector state generally involves the computationally expensive task of performing 2 n−k matrix multiplications on 2 k temporary sub-vectors and aggregating the result (only in special cases can this be done in-place) [90,91].The analogous application of an operator to a density matrix is more expensive due to the extra dimension involved.However, as discussed in section Fast and scalable quantum network simulation of the Results and Discussion, the repetitive nature of NetSquid simulations allows us to take advantage of operators frequently being applied to the same qubit indices for states of a given size.For these operators, we compute a 2 n × 2 n dimensional sparse matrix representation of the k qubit operator via tensor products with the identity and memoize this result for the specific indices and size.When the memoization is applicable the com-putational cost of applying a quantum operator can then be reduced to just sparse matrix multiplication onto a dense vector or matrix.Memoization is similarly applicable to general Clifford operators in the stabiliser tableau formalism.To use memoization on operators that depend on a continuous parameter, such as arbitrary rotations, the parameter can be discretised i.e. rounded to some limited precision.

Physical modelling of network components
All physical devices in a quantum network are modelled by a "component" object, and are thereby also all simulation entities, as shown in Figure 10.Components can be composed of subcomponents, which makes setting up networks in NetSquid modular.
The network itself, for instance, can be modelled as a composite component containing "node" and "connection" components; these composite components can in turn contain components such as quantum memories, quantum and classical channels, quantum sources, etc., as illustrated in Figure 1.The physical behaviour of a component is described by composing it of "models", which can specify physical characteristics such as transmission delays or noise such as photon loss or decoherence.Communication between components is facilitated by their "ports", which can be connected together to automatically pass on messages.
NetSquid also allows precise modelling of quantum computation capable devices.For this it provides the "quantum processor" component, a subclass of the quantum memory.This component is capable of executing "quantum programs" i.e. sequences of "instructions" that describe operations such as quantum gates and measurements or physical processes such as photon emission.Quantum programs fully support conditional and iterative statements, as well as parallelisation if the modelled device supports it.When a program is executed its instructions are mapped to the physical instructions on the processor, which model the physical duration and errors associated to carrying out the operation.A physical instruction can be assigned to all memory positions or only to a specific position, as well as directionally between specific memory positions in the case of multi-qubit instructions.

Asynchronous framework for programming protocols
NetSquid provides a "protocol" class to describe the network protocols and classical control plane logic running on a quantum network.Similarly to the component class, a protocol is a simulation entity and can thereby directly interact with the event timeline.Protocols can be nested inside other protocols and may describe both local or remote behaviour across a network.The "node protocol" subclass is specifically restricted to only operating lo-cally on a single node.Inter-protocol communication is possible via a signalling mechanism and a request and response interface defined by the "service protocol" class.Protocols can be programmed using both the standard callback functionality of PyDy-nAA and a tailored asynchronous framework that allows the suspension of a routine conditioned on an "event expression"; for example, to wait for input to arrive on a port, a quantum program to finish, or to pause for a fixed duration.
The "util" sub-package shown in Figure 10 provides a range of utilities for running, recording and interacting with simulations.Functions to control the simulation are defined in the "simtools" module, including functions for inspecting and diagnosing the timeline.A "data collector" class supports the event-driven collection of data during a simulation, which has priority over other event handlers to react to events.The "simstats" module is responsible for collecting a range of statistics during a simulation run, such as the number of events and callbacks processed, the maximum and average size of manipulated quantum states, and a count of all the quantum operations performed.Finally, the "simlog" module allows fine grained logging of the various modules for debugging purposes.

Benchmarking
To perform the benchmarking described in section Fast and scalable quantum network simulation of the Results and Discussion we used computing nodes with two 2.6 GHz Intel Xeon E5-2690 v3 (Haswell) 12 core processors and 64 GB of memory.Because each process only requires a single core, care was taken to ensure sufficient cores and memory were available when running jobs in parallel.The computation time of a process is the arithmetic average of a number of successive iterations; to avoid fluctuations due to interfering CPU processes the reported time is a minimum of five such repeated averages.To perform the simulation profiling the Cython extension modules of both NetSquid and PyDynAA were compiled with profiling on, which adds some runtime overhead.Version 0.10.0 and 0.3.5 of NetSquid and PyDynAA were benchmarked.We benchmarked against ProjectQ version 0.4.2 using its "MainEngine" backend.See Supplementary Note 2 for further details.
Using the same machine, simulations for B. Implementing a processing-node repeater chain in NetSquid Here, we explain the details of the most complex of our three use cases, namely the repeater chain of Nitrogen-Vacancy-based processing nodes from section Sensitivity analysis for the physical modelling of a long range repeater chain of the Results and Discussion (see Supplementary Notes 3 and 7 for details on the other two use cases).We first describe how we modelled the NV hardware, followed by the repeater protocols used.With regard to the physical modelling, let us emphasise that this is well established (see e.g.[92]); the main goal here is to explain how we used this model in a NetSquid implementation.
In our simulations the following NetSquid components model the physical repeater chain: "nodes", each holding a single "quantum processor" modelling the NV centre, and "classical channels" that connect adjacent nodes and are modelled as fibres with a constant transmission time.We choose equal spacing between the nodes.If we were to simulate individual attempts at entanglement generation, we would also need components for transmitting and detecting qubits such as was used in previous Net-Squid simulations of NV centres [39].However, in order to speed up simulations we insert the entangled state between remote NVs using a model.We designed two types of protocols to run on each node of this network that differ in whether they implement a scheme with or without distillation.
In the remainder of this section, we describe the components modelling.More detailed descriptions of the hardware parameters and their values used in our simulation can be found in Supplementary Note 4.

Modelling a nitrogen-vacancy centre in diamond
In NetSquid, the NV centre is modelled by a quantum processor component, which holds a single communication qubit (electronic spin-1 system) and multiple storage qubits ( 13 C nuclear spins).The decay of the state held by a communication qubit or storage qubit is implemented using a noise model, which is based on the relaxation time T 1 and the dephasing time T 2 .If a spin is acted upon after having been idle for time ∆t, then to its state ρ we first apply a quantum channel where and p = 1 − e −∆t/T1 .Subsequently, we apply a dephasing channel where Z = |0 0| − |1 1| and the dephasing probability equals p = 1 2 1 − e −∆t/T2 • e ∆t/(2T1) .
The electron and nuclear spins have different T 1 and T 2 times.We allow the quantum processor to perform the following operations on the electron spin: initialisation (setting the state to |0 ), readout (measurement in the {|0 , |1 } basis) and arbitrary single-qubit rotation.In particular, the latter includes Pauli rotations where θ is the rotation angle, P ∈ {X, Y, Z} and For the nuclear spin, we have only initialisation and rotations R Z (θ) for arbitrary rotation angle θ.In addition, we allow the two-qubit controlled-R X (±θ) gate between an electron (e) and a nuclear (n) spin: We model each noisy operation O noisy as the perfect operation O perfect followed by a noise channel N : If O is a single-qubit rotation, then N is the depolarising channel: is the probability that a measurement outcome 0 (1) is flipped to 1 (0).

Simulation speedup via state insertion
For generating entanglement between the electron spins of two remote NVs, we simulate a scheme based on single-photon detection, following its experimental implementation in [93].NetSquid was used previously to simulate each generation attempt of this scheme, which includes the emission of a single photon by each NV, the transmission of the photons to the midpoint through a noisy and lossy channel, the application of imperfect measurement operators at the midpoint, and the transmission of the measurement outcome back to the two involved nodes [39].For larger internode distances, simulating each attempt requires unfeasibly long simulation times due to the exponential decrease in attempt success rate.To speed up our simulations in the examples studied here, we generate the produced state between adjacent nodes from a model which has shown good agreement with experimental results [93].This procedure includes a random duration and noise induced on the storage qubits, as we describe below.Let us define where p det is the detection probability, p dc the dark count probability, V denotes photon indistinguishability and α is the bright-state parameter (see Supplementary Note 4 for parameter descriptions).We follow the model of the produced entangled state from the experimental work of [93], whose setup consists of a beam splitter with two detectors located between the two adjacent nodes.In their model, the unnormalised state is given by where ± denotes which of the two detectors detected a photon (each occurring with probability 1 2 ).We also follow the model of [93] for double-excitation noise and optical phase uncertainty, by applying a dephasing channel to both qubits with parameter p = p dexc /2, followed by a dephasing channel of one of the qubits, respectively.The success probability of a single attempt is p succ = p 00 + p 01 + p 10 + p 11 .
The time elapsed until the fresh state is put on the electron spins is (k − 1) • ∆t with ∆t := (t emission + L/c), where t emission is the delay until the NV centre emits a photon, L the internode distance and c the speed of light in fibre.Here, k is the number of attempts up to and including successful entanglement generation and is computed by drawing a random sample from the geometric distribution Pr(k) = p succ • (1 − p succ ) k−1 .After the successful generation, we wait for another time ∆t to mimic the photon travel delay and midpoint heralding message delay.Every entanglement generation attempt induces dephasing noise on the storage qubits in the same NV system.We apply the dephasing channel (eq.( 1)) at the end of the successful entanglement generation, where the accumulated dephasing probability is where p single is the single-attempt dephasing probability (see eq. ( 46) in Supplementary Note 4).

How we choose improved hardware parameters
Here, we explain how we choose 'improved' hardware parameters.Let us emphasise that this choice is independent of the setup of our NetSquid simulations and only serves the purpose of showcasing that NetSquid can assess the performance of hardware with a given quality.By 'near-term' hardware, we mean values for the above defined parameters as expected to be achieved in the near future by NV hardware.If we say that an error probability is improved by an improvement factor k, we mean that its corresponding no-error probability equals k √ p ne , where p ne is the no-error probability of the near-term hardware.For example, visibility V is improved as k √ V while the probability of dephasing p of a gate is improved as 1− k √ 1 − p.A factor k = 1 thus corresponds to 'near-term' hardware.By 'uniform hardware improvement by k', we mean that all hardware parameters are improved by a factor k. We do not improve the duration of local operations or the fibre attenuation.The near-term parameter values as well as the individual improvement functions for each parameter can be found in Supplementary Note 4.

NV repeater chain protocols
For the NV repeater chain, we simulated two protocols: swap-asap and nested-with-distill.Both protocols are composed of five building blocks: entgen, store, retrieve, distill and swap.By entgen, we denote the simulation of the entanglement generation protocol based on the description in the previous subsection: two nodes wait until a classical message signals that their respective electron spins hold an entangled pair.In reality, such functionality would be achieved by a link layer protocol [39].store is the mapping of the electron spin state onto a free nuclear spin, and retrieve is the reverse operation.The distill block implements entanglement distillation between two remote NVs for probabilistically improving the quality of entanglement between two nuclear spins (one at each NV), at the cost of reading out entanglement between the two electron spins.It consists of local operations followed by classical communication to determine whether distillation succeeded.The entanglement swap (swap) converts two short-distance entangled qubit pairs A − M and M − B into a single long-distance one A − B, where A, B and M are nodes.It consists of local operations at M , including spin readout, and communicating the measurement outcomes to A and B, followed by A and B updating their knowledge of the precise state A − B they hold in the perfect case.We opt for such tracking as opposed to applying a correction operator to bring A − B back to a canonical state since the correction operator generally cannot be applied to the nuclear spins directly.Details of the tracking are given in Supplementary Note 6.The circuit implementations for the building blocks, "quantum programs" in NetSquid, are given in Supplementary Note 5. Let us explain the swap-asap and nested-withdistill protocols in spirit; the exact protocols run asynchronously on each node and can be found in Supplementary Note 5.In the swap-asap protocol, a repeater node performs entgen with both its neighbours, followed by swap as soon as it holds the two entangled pairs.Next, nested-with-distill is a nested protocol on 2 n + 1 nodes (integer n ≥ 0) with distillation at each nesting level which is based on the BDCZ protocol [9].For nesting level n = 0, there are no repeaters and the two nodes only perform entgen once.For nesting level n > 0, the chain is divided into a left part and a right part of 2 n−1 + 1 nodes, and the middle node (included in both parts) in the chain generates twice an entangled pair with the left end node following the (n − 1)-level protocol; store is applied in between to free the electron spin.Subsequently, distill is performed with the two pairs as input (restart if distillation fails), after which the same procedure is performed on the right.Once the right part has finished, the middle node performs swap to connect the end nodes.If needed, store and retrieve are applied prior to distill and swap in order achieve the desired configuration of qubits in the quantum processor, e.g. for distill to ensure that the two involved NVs hold an electron-electron and nuclearnuclear pair of qubits, instead of electron-nuclear for both entangled pairs.

IV. DATA AVAILABILITY
The data presented in this paper have been made available at https://doi.org/10.34894/URV169 [94].

V. CODE AVAILABILITY
The NetSquid-based simulation code that was used for the simulations in this paper has been made available at https://doi.org/10.34894/DU3FTS [95].
This section supplements the Methods, section Design and functionality of NetSquid, by going into more depth on specific details of NetSquid's design.The version of NetSquid that we consider is 0.10.For up-todate documentation of the latest NetSquid version, including a detailed user tutorial, code examples, and its application programming interface, please visit the NetSquid website: https://netsquid.org[66].

A. Qubits and their quantum state formalisms
The qubits sub-package of NetSquid, shown in Figure 10 (main text), provides a specialised quantum computation library for tracking the lifetimes of many qubits across a quantum network.A class diagram of the main classes present in this sub-package is shown in in Supplementary Figure 11.Rather than assigning a single quantum state for a predefined number of qubits, both the number of qubits and the quantum states describing them are managed dynamically during a simulation run.Every Qubit (Qubit) object references a shared quantum state (QState) object, which varies in size according to the number of qubits sharing it.When two or more qubits interact, for instance via a multi-qubit operation, their respective shared quantum states are merged together.On the other hand, when a qubit is projectively measured or discarded it can be split from the quantum state it's sharing and optionally be assigned a new single-qubit state.The QState class is an interface for shared quantum states that NetSquid implements for four different quantum state formalisms -described in more detail below.To allow simulations to seamlessly switch between formalisms NetSquid offers a formalism agnostic API, which is defined in the qubitapi module.The functions in this API take as their primary input parameters the qubits to manipulate and the operators (Operator) describing a quantum operation to perform, if applicable.The merging and splitting of shared quantum states is handled automatically under the hood, as are conversions between states using different formalisms (where this is possible).This allows users to program in a "qubit-centric" way, by for instance applying local operations to qubits at a network node without knowledge of their positions within a quantum state representation or any entanglement they may have across the network.We proceed to give a high-level description of the available quantum state formalisms.The first two formalisms are ket state vectors (KET) and density matrices (DM), which both enable universal quantum computation.A ket state vector represents a quantum pure state, while a density matrix can represent statistical ensembles of pure states.The stabiliser formalism (STAB) [77,96] and graph states with local Cliffords formalisms (GSLC) [78] can only represent stabiliser states.Stabiliser states form a subset of all quantum states that are closed under the application of: • Clifford gates.Each Clifford gate can be written as circuit consisting of the following three gates only: a Can only represent stabiliser states.The only operators that can operate on these states are Clifford operators.b A stricter upper bound exists, depending on the matrix multiplication implementation.

Stabiliser tableaus (STAB)
In the stabiliser formalism [96], one tracks the generators of the stabiliser group of a state.We briefly explain the concept here; for a more accessible introduction to the topic, we refer to [97].In order to define a stabiliser group, let us give the Pauli group, which consists of strings of Pauli operators with multiplicative phases ±1, ±i: A stabiliser group is a subgroup of the Pauli group which is commutative (i.e.any two elements A and B satisfy A • B = B • A) and moreover does not contain the element −1 In case the stabiliser group contains 2 n elements, there is a unique quantum state |ψ for which each element A from the stabiliser group stabilises |ψ , i.e.A |ψ = |ψ .Not all quantum states have such a corresponding stabiliser group; those that do are called stabiliser states.The intuition behind the stabiliser state formalism is that one tracks how the stabiliser group is altered by Clifford operations and |0 / |1 -basis measurements.Since the stabiliser state belonging to a stabiliser group is unique, one could in principle always convert the group back to any other formalism, such as KET.Concrete examples of stabiliser groups and their corresponding stabiliser states are: • the stabiliser group {1 1 2 , Z}, which corresponds to the state |0 ; • the stabiliser group {1 Rather than tracking the entire 2 n -sized stabiliser group, it suffices to track a generating set, i.e. a set of n Pauli strings whose 2 n product combinations yield precisely the 2 n elements of the stabiliser group.The choice of generators is not unique.For the examples given above, example sets of stabiliser generators are: • for |0 , the stabiliser group is generated by the single element Z, since Z 2 = 1 1 2 • for |00 , the stabiliser group is generated by {Z ⊗ 1 1 2 , 1 1 2 ⊗ Z}, since squaring any of these two yields 1 1 2 ⊗ 1 1 2 , while multiplying them yields Z ⊗ Z; • for the state (|00 + |11 )/ √ 2, one possible set of of generators is {X ⊗ X, Z ⊗ Z}.
In NetSquid we store generators as a stabiliser tableau: where p k , x jk , z jk ∈ {0, 1}, 0 < j, k ≤ n The k-th generator corresponds to the k-th row of this tableau and is given by For updating the stabiliser tableau after the application of a Clifford gate or a |0 / |1 -basis measurement, NetSquid uses the algorithms by [96] and [77].The runtime performance of stabiliser tableau algorithms is a direct function of the number of qubits: linear for applying single-or two-qubit Clifford unitaries, which any Clifford can be composed into, and cubic for single-qubit measurement [96].

Graph states with local Cliffords (GSLC)
The last formalism is GSLC: graph states with local Cliffords [78].Graph states are a subset of all stabiliser states (see [98] for a review) and an n-qubit graph state |ψ can be written as where Z jk indicates a controlled-Z gate |00 00| + |01 01| + |10 10| − |11 11| between qubits j and k, and we have denoted As such, a graph state is completely determined by the set of qubit index pairs (j, k) at which a controlled-Z operation is performed.These indices can be captured in a graph with undirected edges; in eq.(11), the edge set is E. Each stabiliser state can be written as a graph state, followed by the application of single-qubit Clifford operations.Thus, a stabiliser state in the GSLC formalism is represented by a set of edges E and a list of n single qubit Cliffords.There exist 24 single-qubit Cliffords, so the Clifford list only requires O(n) space.For updating the graph and the list of single-qubit Cliffords after the application of a Clifford gate or a |0 / |1 -basis measurement, NetSquid uses the algorithms by [78].The runtime scaling of the graph-state-based formalism depends on the edge degree d of the vertices involved in the operation -constant-time for single-qubit Cliffords, quadratic in d for two-qubit Cliffords and measurement -and thus scales favourably if the graph is sparse.

B. The PyDynAA simulation engine
The discrete-event modelling framework used by NetSquid is provided by the Python package PyDynAA, which is based on the core engine layer of DynAA, a system analysis and design tool [80].This foundation provides a simple yet powerful language for describing large and complex system architectures.To realise PyDynAA, the simulation engine core was written in C++ for increased performance, and bindings to Python were added using Cython.NetSquid takes advantage of the Cython headers exposed by PyDynAA to efficiently integrate the engine into its own compiled C extension libraries.Several of NetSquid's sub-packages depend and build on the classes provided by PyDynAA, as illustrated in Figure 10 (main text).In Supplementary Figure 13 we highlight several of these key classes and how they interact with the simulation timeline in more detail, namely: the simulation engine (SimulationEngine), events (Event and EventType), simulation entities (Entity), and event handlers (EventHandler).We proceed to describe the concepts these classes represent in more detail.Simulation entities represent anything in the simulation world capable of generating or responding to events.They may be dynamically added or removed during a simulation.The Entity superclass provides methods for scheduling events to the timeline at specific instances and waiting for them to trigger.The intended use is that users subclass the Entity class to implement their own entities.The simulation engine efficiently handles the scheduling of events at arbitrary (future) times by storing them in a self-balancing binary search tree.Events may only be scheduled by entities, which ensures that events always have a source entity.If an entity is removed during a simulation, then any future events it had scheduled will no longer trigger.An entity responds to events by registering an event handler object with a callback function.Responses can be associated to a specific type, source, and id (including wildcard combinations).The simulation engine runs by stepping sequentially from event to event in a discrete fashion and checking if any event handlers in its registry match.A hash table together with an efficient hashing algorithm ensure efficient lookups of the event handlers in the registry.
PyDynAA implements an event expression class to allow entities to wait on logical combinations of events.Atomic event expressions, which describe regular wait conditions for standard events, can be combined to form composite expressions using logical and and or operators to any depth.Event expressions enable NetSquid simulations to deal with timing complexities.This feature has been used extensively in NetSquid to model both the internal behaviour of hardware components, as well as for programming network protocols.As example, consider DEJMPS entanglement distillation [67]: two nodes participate in this protocol and a node can only decide whether the distillation succeeded or failed when both its local quantum operations have finished and it has received the measurement outcome from the remote node.Thus, the node waits for the logical and of the receive-event and the event that the local operations have finished.

C. The modular component modelling framework
The physical modelling of network devices is provided by several NetSquid sub-packages: components, models and nodes, which are shown stacked with relation the NetSquid package in Figure 10 (main text).The pivotal base class connecting all them is the component (Component), which is used to model all hardware devices.Specifically, it represents all physical entities in the simulation, and as such sub-classes the entity (Entity),which enables it to interact with the event timeline.In Supplementary Figure 14  can be represented by a single component, which is composed of node and connection sub-components, which in turn are composed of devices such as channels, sources, memories, etc.To streamline and automate the communication between components, including to and from sub-components, components can be linked using ports (Port) that can send, receive and forward both quantum and classical messages (Message).While the component base class defines a modular interface for modelling all kinds of hardware, it doesn't internally implement any event-driven behaviour itself.That behaviour is implemented by a library of base classes that sub-class Component.The right half of Supplementary Figure 14 shows the sub-classing hierarchy of the provided components, ranging from quantum and classical channels, quantum memory and processing devices, sources, detectors, clocks, to nodes, connections, and networks.The quantum processor (QuantumProcessor) is a component from the base class library used for modelling general quantum processing devices.It sub-classes the quantum memory (QuantumMemory) component, from which it inherits a collection of quantum memory positions (MemoryPosition) for tracking the quantum noise of stored qubits.The processor can assign a set of physical instructions to these positions to describe the operations possible for manipulating their stored qubits, such as quantum gates and measurements, or initialisation, absorption, and emission processes.The physical instructions map to general device-independent instructions, for which they specify physical models such as duration and error models specific to the modelled device.This mapping allows users to write quantum programs in terms of device-independent instructions and re-use them across devices.The quantum programs can include classical conditional logic, make use of parallel execution (if supported by the device), and import other programs.

D. Asynchronous programming networks using protocols
While components are entities in the simulation describing physical hardware, protocols -represented by the Protocol base class as shown in Supplementary Figure 15 -are entities that describe the intended virtual behaviour of a simulation.In other words, the protocol base class is used to model the various layers of software running on top of the components at the various nodes and connections of a network.That can include, for instance, any automated control software at the physical or link layers of a quantum network stack, up to higher-level programs written at the application layer.Protocols in NetSquid can be likened to background processes: they can be started, stopped, as well as reset to clear any state.They can also be nested i.e. a protocol can manage the execution of sub-protocols under its control.To communicate changes of state, such as a successful or failed run, protocols can use a signalling mechanism (Signal).
NetSquid defines several sub-classes of the protocol base class that add extra restrictions or functionality.
To restrict the influence of a protocol to only a local set of nodes the local protocol (LocalProtocol) can be used.Similarly, to restrict a protocol to executing on only a single node, which is a typical use case, a node protocol (NodeProtocol) is available.The service protocol (ServiceProtocol) describes a protocol interface in terms of the types of requests and responses they support.Lastly, a data node protocol adds functionality to process data arriving from a port linked to a connection, and the timed node protocol supports carrying out actions at regularly timed intervals.Programming a protocol involves waiting for and responding to events, which is achieved in the simulation engine by defining event handlers that wrap callback functions.As the complexity of a protocol grows, typically the flow and dependencies of the callback calls do too.To make the asynchronous interaction between protocol and component entities easier and more intuitive to program and read, the main execution function of a protocol (the run() method) can be suspended mid-function to wait for certain combinations of events to trigger.This is implemented in Python using the yield statement, which takes as its argument an event expression.Several helper methods have been defined that generate useful event expressions a protocol can await, for instance: await_port_input() to wait for a message to arrive on a port, or await_timer() to have the protocol sleep for some time.
Figure 16: Circuits used to benchmark quantum computation in the Results, section Benchmarking of quantum computation, for n qubits.For panel (b) the CNOT control line crossing the ellipses represents multiple lines for n > 6 qubits, following the pattern of q 2 and q 3 .Similarly, the classical control lines represent an AND of the measurement results for q 3 , q 5 , . . ., q n−1 and q 2 , q 4 , . . ., q n−2 to determine the control of the X and Z gates, respectively.The noise gates denoted by {X,Y,Z} cycle through the Pauli gates (see main text).Note that this circuit always requires an even number of qubits.
To benchmark the runtimes of quantum computation circuits the processes were timed in isolation from any setup code using the Python timeit package.Python garbage collection is disabled during the timing of each process.To avoid fluctuations due to interfering CPU processes the reported time is a minimum of five repetitions.

B. Runtime profiling of a repeater chain simulation
The runtime profiling of NetSquid presented in the Results, section Benchmarking of event-driven simulations, is performed for a simple repeater chain.The network setup of this simulation extends the single repeater presented in Supplementary Figure 1 to a chain of nodes by adding the entangling connection shown between each pair of neighbouring nodes.Direct classical connections are connected between each node and one of the end-nodes, rather than between neighbouring nodes, and are used to transmit the swapping corrections.The chosen configuration for this network does not need to be physically realistic; it suffices for it to be representative of the typical computational complexities.The nodes are placed at 20km intervals and the channels transmit messages at the speed of light in fibre.The entanglement sources, assumed to be perfect, are all synchronised and operate at a frequency of 100 kHz.Physical non-idealities are represented by adding time-dependent depolarising noise to both the quantum channels and quantum memories, as well as dephasing noise to quantum gates.The corresponding depolarising and dephasing times are 0.1 s and 0.04 s, which correspond to the T 1 and T 2 times presented in section Modelling a nitrogen-vacancy centre in diamond of the Methods.In a simulation run entanglement is created once between the end-nodes by performing entanglement swaps along the chain.Protocols are assigned to all but the end-nodes to perform entanglement swaps after each round of entanglement generation, and send their measurement results as corrections to the same end-node.
A protocol running on the end-node collects these corrections, and applies them if needed.The runtime of this simulation is profiled to determine the distribution of time spent in the functions of NetSquid's sub-packages, as well as its dependency packages NumPy and PyDynAA.To perform this profiling the cProfile package is used.The reported runtime for a given number of nodes is the mean of 400 successive simulation runs.
Electron readout (eq.( 4) and sec.• Dark counts: a photon detector falsely registering.The dark count probability follows a Poisson distribution p dc = 1 − e −tw•λ dark where t w = 25 ns [93] is the duration of the time window at the midpoint.We set λ dark = 1 Hz as the dark count rate.
• Imperfect photon indistinguishability.The generation of entanglement at the middle station is based upon the erasure of the which-way information with respect to the path of the photons.Only in case the photons are fully indistinguishable, the which-way information is erased perfectly.The overlap of the photon states is given by the visibility V , which we set to 0.9 [93].
• Double excitation of the electron spin.When triggered to emit a photon by a resonant laser pulse, an NV centre could be excited twice, which results into the emission of two photons.We set its occurrence probability to p dexc = 0.06 [99].
• Photon phase uncertainty.The photons which interfere at the midpoint acquired a phase during transmission and a difference of these phases influences the precise entangled state that is produced [104].Given a standard deviation σ phase = 0.35 rad [99] of the acquired phase, we compute the dephasing probability as [93] c. Nuclear spin dephasing during entanglement generation The initialisation of the electron spin state induces dephasing of the carbon spin states through their hyperfine coupling.Following [105], we model this uncertainty by a dephasing channel for each attempt with dephasing probability The parameter C nucl is the product of the coupling strength between the electron spin and the carbon nuclear spin, and an empirically determined decay constant.Rather than expressing the dephasing probability as function of C nucl , we express the magnitude of nuclear dephasing as N 1/e , the number of electron spin pumping cycles after which the Bloch vector length of a nuclear spin in the state (|0 + |1 )/ √ 2 in the X − Y plane of the Bloch sphere has shrunk to 1/e, when the electron spin state has bright-state parameter α = 0.5 (i.e the electron spin is in the state (|0 + |1 )/ √ 2).Let us compute how p single depends on N 1/e instead of on C nucl .First, we find by direct computation that the equatorial Bloch vector length of a state is shrunk by a factor 1 − 2p after a single application of the single-qubit dephasing channel (eq.( 1)) with probability p (p ≤ 1 2 ).Equating Equating p single from eq. ( 44) with α = 0.5 and p from eq. ( 45), followed by solving for C nucl yields Substituting back into eq.( 44) yields an expression for general α: We set N 1/e = 1400 [106].

d. Local processing parameters
For the dynamics of the electron spin, we use T 1 = 1 hour and T * 2 = 1.46s [107].For the carbon nuclear spin, we take T 1 = 10 hours and T 2 = 1 s (experimentally realised: T 1 = 6m and T 2 ≈ 0.26 − 25s [108]).For the noise of the controlled-R X gate (Methods, section Modelling a nitrogen-vacancy centre in diamond), we set the depolarising probability p = 0.02 (denoted as p EC in Supplementary Table II), since by simulation of the circuit [104, Fig. 2a], we find that this value agrees with the experimentally found effective circuit Figure 17: Quantum circuits used in simulations of the NV repeater chain, acting on an electron (e) and nuclear (n) spin.when its previous action is finished.A simulation run finishes as soon as the two end nodes share a single entangled pair of qubits.
For the swap-asap protocol, the sequence of operations that a node performs depends on its index in the chain (start counting from left to right, nodes have indices 1, 2, 3, . . .).If the index of the node is even, the node sends a request for entgen to its left neighbour, and starts the operation as soon as it has received confirmation.After performing store to free the electron spin, it repeats it for its right neighbour.Oddindexed nodes remain idle until reception of an entgen request, after which they perform store if necessary to free the electron, send a confirmation, sleep for the duration of the message transmission, followed by performing entgen.Once a node has entanglement with both directions, it performs a swap and sends the outcome to the end nodes.
The two end nodes are exceptions to the above.The left end node (i.e. with index 1) behaves like an oddindexed node, but without performing swap.The same holds for the rightmost node (i.e. the node with the largest index), unless its index is even, in which case it initiates and performs entanglement generation with the adjacent node on the left.The nested-with-distill protocol is a variant of the BDCZ protocol [9], adapted to the fact that an NV cannot perform multiple entgen , distill or swap operations in parallel due to its restricted topology (Methods, section Modelling a nitrogen-vacancy centre in diamond).In the adapted version, nodes take the role of initiator of one of the three main actions (entgen, distill, swap) if the action occurs at the highest nesting level that this node belongs to.To be precise, we do the following.In a repeater chain with 2 n + 1 nodes, denote by {0, 1, 2, . . ., 2 n } the indices of the nodes from left to right.A node (not an end node) with index k ∈ {1, 2, 3, . . ., 2 n − 1} initiates an action only if the entanglement that is involved in the task spans precisely f n (k) segments, where End nodes are never initiators.When a node (index k) is triggered, it performs the following checks in order and performs the first action for which the check holds true: Storing and retrieving qubits.Locally mapping the state of a qubit onto a different memory position by the store or retrieve circuits does not alter the correction Pauli corresponding to that qubit.
Entanglement distillation.Suppose that nodes A and B wish to perform the distill protocol, which starts by A and B sharing an electron-electron pair (correction Paulis P e A and P e B at node A and B, respectively) and a nuclear-nuclear pair (P n A and P n B ).In the protocol, first both nodes apply P n • P e to their electron spin qubit.Then, both nodes locally perform the distillation circuit from Supplementary Figure 17(c), followed by sending both the measurement outcome and P n to the other node.The nodes determine whether the distillation succeeded using the condition explained in section 5.In case of failure, the nuclear-nuclear state is discarded.In case of success, one of the nodes in the chain (for example, the one with the lower position index in the chain) sets P n = 1 1, while the other sets Below, in section 6 B of this Supplementary Note, we show that after this procedure, eq. ( 62) still holds.
Entanglement swapping.Suppose that node M wants to execute the swap protocol on shared pairs A − M and M − B, with nodes A and B respectively.We denote M 's correction Paulis as P A M and P B M .First, M performs the Bell-state measurement circuit from Supplementary Figure 17 Both nodes A and B multiply their local Pauli with the Pauli they received from M .The proof that after the swap, eq. ( 62) still holds can be found below in section 6 C of this Supplementary Note.

B. Correctness proof of the correction operator update for distill
Here, we prove that eq. ( 62) holds for the states that are outputted by the protocols for entanglement distillation and swapping explained above.Let us start with entanglement distillation.For this, we denote by 'physical nuclear-nuclear state' the joint state of the nuclear spins of node A and B. By direct computation, one can show the following.We emphasise that using the correction-operator tracking for the store and retrieve operations as described in section 6 A of this Supplementary Note, the physical nuclear-nuclear state between any two nodes does not satisfy eq. ( 62).The reason for this is that the store operation maps the electron spin state to the nuclear spin in a rotated basis, where the rotation operator is a Hadamard gate H (eq. 51).However, the correction operators are not updated when the store is applied (see 'Storing and retrieving qubits' in section 6 A).Consequently, if nodes A and B share the physical nuclear-nuclear state |ψ , then mapping |ψ protocol and atomic frequency comb (AFC) memories as an example for an engineered absorption protocol.For a more detailed comparison of these two types of memories see [112].AFC memories are very promising due to their great potential for multiplexing, while EIT memories promise superior efficiency with limited multiplexing.Here we will give a brief overview of both technologies.
a. AFC AFC memories require an inhomogeneously broadened material (e.g.rare-earth-doped crystals) with equidistant peaks in occupation density created by spectral hole burning [72].When a photon is absorbed by the memory, the system will be in a collective delocalised state (Dicke state [117]) which then rapidly dephases with time.After a fixed time the Dicke state rephases and the absorbed photon is re-emitted.The retrieval efficiency decreases exponentially with storage time.Due to the large linewidth of the inhomogeneously broadened state, AFC memories have great potential for massive multiplexing as the number of modes is not limited by the optical depth of the material.Currently, the efficiencies that have been experimentally achieved are up to 58% [118] reported for a cavity-enhanced AFC protocol and around 8.5% [119,120] for most multiplexed memories.Memory lifetimes of up to 0.53 s [121] have been reported for a dynamically decoupled AFC protocol (however with an efficiency of just 0.5%).In our simulation we set a maximum memory efficiency of 45% (currently optimistic value for multimode memories, but still below the highest demonstrated single-mode storage [118]) with 1000 modes, consisting of 50 spectral and 20 temporal modes (15 x 9 already demonstrated [119], over 1000 total modes also demonstrated [122]), and a lifetime of 1 ms (0.542 ms already achieved for multimode on-demand AFC memory [120]).

b. EIT
EIT [123] is an effect where an opaque sample of atomic-ensembles becomes transparent due to interaction with electromagnetic fields.This creates a steep dispersion which in turn affects the group velocity of incoming light.The effect can be used to slow or completely stop and later re-emit incoming photonic states.Due to the use of a strong resonant control laser to store the incoming light, these kind of memories are subject to background noise in the form of additional photons introduced by the control beam.These additional photons need to be compensated by filtering.Another disadvantage is that multiplexing in EIT is limited to using only the spatial degree of freedom, thus making it more difficult to realise a highly multiplexed protocol.The great advantage of these memories is the high efficiency of 85% [124] at a memory lifetime of 15 µs.In our simulation we set a maximum efficiency of 90% (within reach of [124]) with 1000 spatial modes (experimentally 2 spatial modes have been demonstrated [125]), and a lifetime of 100 µs ( [124] gives 200µs as a theoretical limit) .For a review of EIT see [73].

C. Parameters
For the comparison we used the optimistic elementary link parameters taken from [19](Fig 10 (g)).We then added optimistic values for the the memory parameters (see above) and the additional noise parameters we use, such as the photon visibility.In line with the analysis of [19], we model the source as an almost perfect single-pair source with a small probability p(2) of emitting two pairs.We first list the common parameters of both simulations before we compare the parameters for the memory components in Supplementary Table III.It is worth pointing out that using the same number of modes for both technologies slightly biases this comparison as the high multiplexing potential is the main advantage of AFC memories.

Figure 1 :
Figure 1: Illustrative example of a NetSquid use case.Each sub-figure explains part of the modelling and simulation process.For greater clarity the figures are not based on real simulation data.The scenario shown is a quantum repeater utilising entanglement distillation (see main text).a) The setup of a quantum network using node and connection components.b) A zoom in showing the subcomponents of the entangling connection component.The quantum channels are characterised using fibre delay and loss models.The quantum source samples from an entangled bipartite state sampler when externally triggered by the classical channel.c) A zoom in of the quantum memory positions within a quantum processor illustrating their physical gate topology.The physical single-qubit instructions possible on each memory in this example are the Pauli (X, Y , Z), Hadamard (H), and X-rotation (R X ) gates, and measurement.The blue-dashed arrows show the positions and control direction (where applicable) for which the two-qubit instructions controlled-X (CNOT) and swap are possible.Noise and error models for the memories and gates are also assigned.d) Illustration of a single simulation run.Time progresses by discretely stepping from event to event, with new events generated as the simulation proceeds.Qubits are represented by circles, which are numbered according to the order they were generated.A star shows the moment of generation.The curved lines between qubits denote their entanglement with the colour indicating fidelity.The state of each qubit is updated as it is accessed during the simulation, for instance to apply time-dependent noise from waiting in memory.e) A zoom in of the distillation protocol.The shared quantum states of the qubits are combined in an entangling step, which then shrinks as two of the qubits are measured.The output is randomly sampled, causing the simulation to choose one of two paths by announcing success or failure.f ) A plot illustrating the stochastic paths followed by multiple independent simulation runs over time, labeled by their final end-to-end fidelity F i .The blue dashed line corresponds to the run shown in panel (d).The runs are typically executed in parallel.Their results are statistically analysed to produce performance metrics such as the average outcome fidelity and run duration.

Figure 2 :
Figure 2: Abstract example of simulating aquantum protocol with discrete events.When setting up the simulation, protocol actions are defined to occur when a specific event occurs, as in the setup of a real system.Instead of performing a continuous time evolution, the simulation advances to the next event, and then automatically executes the actions that should occur when the event takes place.Any action may again define future events to be triggered after a certain (stochastic) amount of time has elapsed.For concreteness a simplified quantum teleportation example is shown where a qubit, shown as an orange circle with arrow, is teleported between the quantum memories of Alice and Bob.Here, entanglement is produced using an abstract source sending two qubits, shown as blue circles with arrows, to Alice and Bob.Once the qubit has traversed the fibre and reaches Alice's lab, an event is triggered that invokes the simulation of Alice's Bell state measurement (BSM) apparatus.The simulation engine steps from event to events defined by the next action, which generally occur at irregular intervals.This approach allows time-dependent physical non-idealities, such as quantum decoherence, to be accurately tracked.

Figure 3 :
Figure3: A quantum switch in a star-shaped network topology as studied by Vardoyan et al.[56].The switch (central node) is connected to a set of users (leaf nodes) via an optical fibre link that distributes perfect Bell pairs at random times, following an exponential distribution with mean rate µ ∝ e −βL , where L denotes the distance of the link and β the attenuation coefficient.Associated with each link there is a buffer that can store B qubits at each side of the link.As soon as n Bell pairs with different leaves are available, the switch performs a measurement in the n-partite Greenberger-Horne-Zeilinger (GHZ) basis, which results in an n-partite GHZ state shared by the leaves.The GHZ-basis measurement consists of: first, controlled-X gates with the same qubit as control; next, a Hadamard (H) gate on the control qubit; finally, measurement of all qubits individually.The figure shows 4 leaf nodes, GHZ size n = 3 and a buffer size B = 2.

Figure 4 :
Figure4: Performance analysis of the quantum switch with 9 users using NetSquid.(a) Capacity as a function of the buffer size (number of quantum memories that the switch has available per user) for either 2− or 4−qubit Greenberger-Horne-Zeilinger (GHZ)-states.For each scenario, the generation rate µ of pairs varies per user.For the blue scenario (2-partite entanglement, µ = [1.9,1.9, 1.9, 1, 1, 1, 1, 1, 1] MHz), the capacity was determined analytically by Vardoyan et al. using Markov Chain methods [56, Figure8].Here we extend this to 4-partite entanglement (orange scenario, same µs), for whichVardoyan et   al. have found an upper bound (by assuming unbounded buffer and each µ = maximum of original rates = 1.9 MHz) but no exact analytical expression.The green scenario (µ = [15, 1.9, 1.9, 1,1, 1, 1, 1, 1] MHz)does not satisfy the stability condition for the Markov chain for unbounded buffer size (each leaf's rate < half of sum of all rates) so in that case steady-state capacity is not well-defined.We note that regardless of buffer size, the switch has a single link to each user, which is the reason why the capacity does not scale linearly with buffer size.(b) Average fidelity of the produced entanglement on the user nodes (no analytical results known) with unbounded buffer size.The fact that the green curve has lower fidelity than the blue one, while the former has higher rates, can be explained from the fact that the protocol prioritises entanglement which has the longest storage time (see Supplementary Note 3).Each data point represents the average of 40 runs (each 0.1 ms in simulation).Standard deviation is smaller than dot size.

Figure 5 :
Figure 5: Performance of repeaters based on nitrogen-vacancy (NV) centres in diamond.(a)Fidelity and entanglement distribution rate achieved with near-term and 10× improved hardware (Supplementary Note 4) with the swap-asap protocol.Dashed line represents classical fidelity threshold of 0.5.We observe that for near-term hardware, the use of 3 repeaters yields worse performance in terms of fidelity than the no-repeater setup.For improved hardware we observe (i) that for approx.0 -750 kms, repeaters improve upon rate by orders of magnitude while still producing entanglement (fidelity> 0.5), while (ii) for approx.750 -1500 kms, repeaters outperform in both rate and fidelity.(b-c) Fidelity and rate achieved without and with repeaters (1, 3 or 7 repeaters) as function of a hardware improvement factor (Methods, section How we choose improved hardware parameters) for two typical distances from both distance regime (i) and (ii), for two protocols swap-asap and nested-with-distill.For the repeater case, only the best-performing number-of-repeaters & protocol in terms of achieved fidelity is shown in (b), accompanied by its rate in (c).Each data point represents the average over (a) 200 and (b) 100 runs.Standard deviation is smaller than dot size.

Figure 6 :
Figure 6: Sensitivity of fidelity in various hardware parameters for nitrogen-vacancy (NV) repeater chains.The NV hardware model consists of 15 parameters and from those we focus on four parameters in this figure: (A) two-qubit gate fidelity, (B) detection probability, (C) induced storage qubit noise and (D) visibility.We start by improving all 15 parameters, including the four designated ones, using an improvement factor of 3 (Methods, section How we choose improved hardware parameters).Then, for each of the four parameters only, we individually decrease their improvement factor to 2, or increase it to 10 or 50.The figure shows the resulting fidelity (horizontal and vertical grid lines; dashed line indicates maximal fidelity which can be attained classically).Note that at an improvement factor of 3 (orange line), all 15 parameters are improved by 3 times, resulting in a fidelity of 0.39.In addition, we vary the improvement factor for combinations of two of the four parameters (diagonal lines).The 3× improved parameter values can be found in Supplementary TableII.The other values (at 2/10/50×) are approximately: two-qubit gate fidelity F EC (0.985/0.997/0.9994),detection probability p nofibre The average computation time for a single run versus the number of qubits in the chain are shown for the different quantum computation libraries in Figure 8(b), where we have again included ProjectQ.

Figure 7 :
Figure 7: Performance comparison of a single quantum repeater with atomic frequency comb (AFC) or electronically induced transparency (EIT) quantum memories.Shown are: (a) thesecret key rate in secret bits per entanglement generation attempt, (b) the quantum bit error rate (QBER) in the X and Z bases (c) the average number of attempts necessary for one successful end-to-end entanglement generation.Each data point is obtained using 10.000 (EIT) or 30.000 (AFC) successful attempts at generating entanglement between the end nodes.Solid lines are fits.Note that for the secret key plot we use logarithmic scale with added 0 at the origin of the axes.Error bars denote standard deviation and are symmetrical.

Figure 8 :
Figure 8: Runtime comparison of NetSquid's quantum state formalisms.Runtime comparisons of the available quantum state formalisms in NetSquid as well as ProjectQ ket vector for two benchmark use cases.The KET, DM, STAB and GSLC formalisms refer to the use of ket vectors, density matrices, stabiliser tableaus and graph states with local Cliffords, respectively.(a) Generating a Greenberger-Horne-Zeilinger (GHZ) state.Qubits are split off from the shared quantum state after a measurement.For the KET formalism the effect of turning off memoization (dotted line) is also shown.(b) Quantum computation involved in a repeater chain.Each formalism is shown with qubits split (dotted lines) versus being kept in-place (solid lines) after measurement.

Figure 10 :
Figure 10: Overview of NetSquid's software architecture.The sub-packages that make up the NetSquid package are shown stacked in relation to each other and the PyDynAA package dependency.The main classes in each (sub-)package are highlighted, and their relationships in terms of inheritance, composition and aggregation are shown.Also shown are the key modules users interact with, which are described in the main text.In this paper NetSquid version 0.10 is described.

Figure 5
(bc) were run, which took almost 260 core hours wallclock time in total, while simulations for Figure 7 took roughly 625 core hours.For Figure 4 (≈10 hours in total), Figure 5(a) (≈90 minutes) and Figure 6 (≈30 minutes), a single core Intel Xeon Gold 6230 processor (3.9GHz) with 192 GB RAM was used.
) with parameter p = 4(1 − F )/3 with F the fidelity of the operation.If O is single-qubit initialisation, N = N depol p with parameter p = 2(1 − F ).The noise map of the controlled-R X gate is an identical single-qubit depolarising channel on both involved qubits, i.e.N = N depol p ⊗ N depol p .Finally, we model electron spin readout by a POVM measurement with the Kraus operators

Figure 13 :
Figure 13: Design overview of the PyDynAA package.Schematic overview of key classes defined by the PyDynAA package, the discrete-event simulation engine used by NetSquid.Also shown is the relation of each class to the simulation timeline.Events are scheduled onto the simulation timeline by Entity objects.Entities wait for events to trigger by registering EventHandlers, which respond to an event by passing it as input to a specified callback function.The events to wait for can be specified by their type, id, and source entity.Ellipses indicate that not all of a class's public variables and methods are listed.Omitted from this class diagram is the EventExpression class -see the text for more details.

Figure 15 :
Figure 15: Design overview of protocols in NetSquid.Class diagram of the Protocol class and its subclasses.Ellipses indicate that not all of a class's public variables and methods are listed.

1 .
If it shares entangled pairs with nodes k − f n (k) and k + f n (k), and both are the immediate result of successful distillation: perform swap and send the measurement outcomes to the involved nodes 2. If it holds two entangled pairs with node k − f n (k) and neither pair is the result from successful entanglement distillation: send a request to distill to the node, wait for confirmation, followed by performing distill clicked (Methods, section Simulation speedup via state insertion).If the +-detector clicked, then the state that A and B hold is the desired state |φ[1, −1] , so we set P A = P B = 1 1.If the other detector clicked, then the produced state is |φ[−1, −1] .Therefore, one of the nodes (for example, the one with the higher position index in the chain) sets the correction operator to Z, whereas the other sets it to 1 1, since (1 1 ⊗ Z) |φ[−1, −1] = |φ[1, −1] .
(d).Let us denote the individual measurement outcomes of the circuit as m earlier and m later (both take values from {1, −1}), which correspond to the measured Bell state |φ[a, b] with a = −1 • m earlier • m later and b = m later .Then, M sends the Pauli 1 1 to A, while to B it sends P A M • P B M • Q, where Q is given by

Proposition 1 .
Suppose that nodes A and B share the state |φ[a, b] on the electrons and the physical nuclear-nuclear state |φ[c, d] , where a, b, c, d ∈ {1, −1}.When both nodes execute the distillation circuit from Supplementary Figure 17(c), the resulting state on the carbon nuclear spins is |φ[c, −a • c • d] and the measurement outcome m 1 ∈ {1, −1} on one side is uniformly random, while the outcome of the other node equals m 2 = m 1 • b • c.
Runtime profile of a repeater chain simulation using Netsquid.Runtime profile for a repeater chain simulation with a varying number of nodes in the chain.The maximum quantum state size is four qubits.The total time spent in the functions of each NetSquid subpackage and its main package dependencies (in italics) is shown.The dark hatched bands show the largest contribution from a single function in each NetSquid sub-package, as well as in NumPy and uncategorised (other ) functions.The sub-packages are stacked in the same order as they are listed in the legend.
the physical layer.However, SQUANCH, similar to SimulaQron and QuNetSim, does not use a simulation engine that can accurately track time.Accurate tracking is crucial for e.g.studying time-dependent noise such as memory decoherence.
Design overview of NetSquid's qubits sub-package.The main classes and module of the netsquid.qubitssub-package.Qubit objects can be manipulated, as described for instance by Operator objects, using the functions of the qubitapi module.Under the hood the qubits share a specific sub-class of the QState interface.Ellipses indicate that not all of a class's public variables and methods are listed.
1n z 11 . . .z 1n p 1 n1 . . .x nn z n1 . . .z nn p n x we show a class diagram of the component class and its relationships to other classes from these sub-packages.The modularity of NetSquid's modelling framework is achieved by the composition of components in terms of properties, models, communication ports and subcomponents.A component's properties are values that physically characterise it, such as the length of a channel or the frequency of a source.A special constrained map (ConstrainedMap) container is used to store the properties (as well as the other composed objects) to give control of the expected types and immutability of properties during a simulation.Models (Model) are used to describe the physical behaviour of a component, such as the transmission delay of a channel, or the quantum decoherence of a qubit in memory.Model objects are essentially elaborate functions and generally do not store any state; when a model is called it is passed its component's properties, in addition to any modelling specific input, such as, in the case of a quantum noise model, the qubit to apply noise and the time the qubit has been waiting on a memory.Components can be composed of other subcomponents, which allows networks to be pieced together in a very modular fashion.For instance, a complete network Design overview of components in NetSquid.Class diagram for the Component class, a simulation entity that is used to model all network hardware devices, including composite components such as nodes, connections and the network itself.A component is shown to be composed of properties, ports, models and subcomponents.Ellipses indicate that not all of a class's public variables and methods are listed.

Table II :
Physical parameters dealing with elementary link generation, memory coherence times and duration and fidelities (F ) of the gates.Depicted are both parameters of the dataset 'near-term' and two examples of improved parameter sets (see Methods, section How we choose improved {X, Y } and P n B ∈ {1 1, Z} Y if P n A ∈ {1 1, Z} and P n B ∈ {X, Y } 1 1 otherwise.